Modeling of multiple regression model with array form

Questions v5
, , ,
1

Hello.

I am a beginner in Bayesian Inference and PYMC (Also Python).
I have studied multi-objective predictive control using machine learning and deep learning.

I wonder whether users can develop the regression (or GP model) with a selective input structure as ANN is available.

What I mean is shown in the below figure:

image
image420Ă—147 3.03 KB

(All Beta and Alpha follow a normal distribution)

So, I am trying to implement that with this code:
But it really takes much duration (I think that this code cannot work. Hence I can’t see even errors)

# multiple linear regression

    x = np.random.normal(10,np.sqrt(2),size=(6,10))
    y = np.random.normal(30,np.sqrt(15),size=(6,10))

    with pm.Model() as mod:
        Temp_array = np.matrix([[1,1,0,1,0,1],
                              [0,1,1,0,1,1],
                              [1,1,1,1,0,0],
                              [0,0,1,1,1,1],
                              [1,0,1,0,1,0],
                              [0,1,1,1,1,0]])
    
        Beta1 = pm.Normal('Beta1', mu=5,sd=1)
        Beta2 = pm.Normal('Beta2', mu=10,sd=3)
        Beta3 = pm.Normal('Beta3', mu=4,sd=1.5)
        Beta4 = pm.Normal('Beta4', mu=7,sd=6)
        Beta5 = pm.Normal('Beta5', mu=9,sd=2)
        Beta6 = pm.Normal('Beta6', mu=11,sd=8)
    
        Alpha1 = pm.Normal('Alpha1', mu=5,sd=1)
        Alpha2 = pm.Normal('Alpha2', mu=10,sd=3)
        Alpha3 = pm.Normal('Alpha3', mu=4,sd=1.5)
        Alpha4 = pm.Normal('Alpha4', mu=7,sd=6)
        Alpha5 = pm.Normal('Alpha5', mu=9,sd=2)
        Alpha6 = pm.Normal('Alpha6', mu=11,sd=8)
    
        Beta_Vec = np.matrix([[Beta1],[Beta2],[Beta3],[Beta4],[Beta5],[Beta6]])
        Alpha_Vec = np.matrix([[Alpha1],[Alpha2],[Alpha3],[Alpha4],[Alpha5],[Alpha6]])
    
        epsilon = pm.Normal('epsilon',mu=0,sd=np.sqrt(5), shape=(6,1))
    
        Beta_Mat = Temp_array*Beta_Vec
    
        Slope = Beta_Mat.transpose() * x + Alpha_Vec
    
        Selective_Regressor = pm.Normal('Selective_Regressor',mu=Slope,observed=y)
        tr = pm.sample(10)

Is there anyone who can help me… ?

I would very much appreciate it if you help me.

1 Like
2

I got error like this:
ValueError: setting an array element with a sequence.

3

Your specific error is because you are mixing numpy and pytensor. PyMC variables are symbolic pytensor objects, not floating point numbers, so numpy has no idea what do with them if you try to put them in e.g. an np.matrix (aside: np.matrix is slated for depreciation, you should use np.array everywhere).

You can give a size parameter to PyMC distributions as a short-hand for sampling a group of variables together. Alternatively, you can pass a sequence to the parameters, and PyMC will automatically infer that you mean you want a vector of RVs of the same size as the associated parameters.

Your code can be written:

mus = np.array([5,10,4,7,9,11])
sigmas = np.array([1,3,1.5,7,2,8])

Beta_Vec = pm.Normal('beta', mu=mus, sigma=sigmas)
Alpha_Vec = pm.Normal('alpha', mu=mus, sigma=sigmas)

You don’t seem to be using epsilon anywhere in your model, what is that? If they are regression errors, that’s covered by the pm.Normal in the observed node. You might want to check out some of the regression examples in the example gallery to better get your bearings.

4

Thank you Jesse! That was just test ^^;; So it was only for checking.

5

Could you give me one more idea?

If I want to make hierarchical linear regression with that forms (What I mean is hyper-prior),
should I use two “with statement” like this?

(1) Sampling hyper-parameters, using pymc
(2) Intermediately, calculate the statistic values from (1)
(3) Getting processed hyper-parameters from (2), and sampling parameters for regression model

^^;;; Sorry for bothering you…

6

It’s not bother, you’re helping me avoid work :wink:

Can you write out what you have in mind in math or pseudo-code? I’m not sure I completely understand what you want to do.

I can say that PyMC uses a symbolic computational backend, so you can do any kind of aggregation, combination, or computation you like with random variables. There’s no need to break up the model context, although you’re allowed to if you wish! It doesn’t do anything. It’s always just quietly building a symbolic graph of your model, which isn’t finalized and complied until you hit the pm.sample button (it has a name but it’s trademarked).

7

Really Thanks ^^.

By the way, last question, I will search the solution by myself. I don’t have a specific plan for implementing this.

However, I get still errors of the first question…
Could you help me again? ㅠㅠ…
(I am actually studying the workbook with PyMC3. So I don’t get used to pytensor…)

# multiple linear regression - Selective
    x = np.random.normal(10,np.sqrt(2),size=(6,10))
    y = np.random.normal(30,np.sqrt(15),size=(6,10))

    with pm.Model() as mod:
        Mu1_array = np.array([[5,10,0,4,0,9],
                              [0,6,2,0,8,11],
                              [13,12,11,10,0,0],
                              [0,0,5,3,7,8],
                              [10,0,11,0,15,0],
                              [0,7,8,6,4,0]])
        Mu2_array = np.array([[5,10,0,4,0,9],
                              [0,6,2,0,8,11],
                              [13,12,11,10,0,0],
                              [0,0,5,3,7,8],
                              [10,0,11,0,15,0],
                              [0,7,8,6,4,0]])
        sd1_array = np.array([[.5,.10,np.inf,.4,np.inf,.9],
                              [np.inf,.6,2,np.inf,.8,.11],
                              [1.3,1.2,1.1,np.inf,np.inf,np.inf],
                              [np.inf,np.inf,.5,.3,7,.8],
                              [np.inf,np.inf,1.1,np.inf,1.5,np.inf],
                              [np.inf,7,8,6,4,np.inf]])
        sd2_array = np.array([[.5,.10,np.inf,.4,np.inf,.9],
                              [np.inf,.6,2,np.inf,.8,.11],
                              [1.3,1.2,1.1,np.inf,np.inf,np.inf],
                              [np.inf,np.inf,.5,.3,7,.8],
                              [np.inf,np.inf,1.1,np.inf,1.5,np.inf],
                              [np.inf,7,8,6,4,np.inf]])
        
        Selective_Regressor = np.ones([6,1])
    
        for i in np.arange(0,5):
            Beta = pm.Normal('Beta', mu=Mu1_array[i,:],sigma=sd1_array[I,:])
            Alpha = pm.Normal('Alpha', mu=Mu2_array[i,:],sigma=sd2_array[I,:])
            epsilon = pm.Normal('epsilon',mu=0,sigma=np.sqrt(5), shape=(6,1))
            Slope = Beta * x + Alpha

            Selective_Regressor[i,0] = pm.Normal(f'Selective_Regressor_{i}',mu=Slope,observed=y)
            tr = pm.sample(10)

The error is : Incompatible Elemwise input shapes [(1, 6), (6, 10)]

8

What line is the error on? What are the shapes of the objects? (1,6) and (6, 10) are big hints – my guess is that (1,6) is Beta and (6,10) is x, so there’s something wrong with Slope = Beta * x + Alpha. Should it be @ instead of *?

You can check all the shapes of your objects by removing the pm.sample, then using pm.model_to_graphviz(mod). This will show you all the random variables in the model, their shapes, and how they connect. It’s a really great tool for quickly debugging models.

A more direct way to check shapes is to use .shape.eval(). For example, print(Beta.shape.eval()) will show you the shape of Beta. .eval() tells pytensor to actually compute the symbolic graph. If you just do print(Beta.shape) you’ll get a symbolic tensor representing the shape of Beta (that’s worth doing once too so you can see how things work).

1 Like
9

Dear Jesse.

I rewrite code referring your answer.

x = np.array([np.random.normal(10,10,size=(1000)),np.random.normal(5,7,size=(1000)),
                  np.random.normal(6,2,size=(1000)),np.random.normal(11,5,size=(1000)),
                 np.random.normal(9,3,size=(1000)),np.random.normal(4,1,size=(1000))])
    y = np.array([np.linspace(20,30,1000),np.linspace(10,40,1000),
                 np.linspace(50,20,1000)])
    
    with pm.Model() as Structured_Model:
        Sigma = np.array([pm.HalfCauchy('Sigma1', beta=1),
                         pm.HalfCauchy('Sigma2', beta=6),
                         pm.HalfCauchy('Sigma3', beta=4)])
        Intercept = np.array([pm.Normal("Intercept1", mu=0, sigma=9),
                              pm.Normal("Intercept2", mu=0, sigma=4),
                              pm.Normal("Intercept3", mu=0, sigma=5)])
        Slope = np.array([[pm.Normal('beta11',10,3),pm.Normal('beta12',1,4),
                           pm.Uniform('beta13',[-np.inf,np.inf]),pm.Normal('beta14',7,1),
                           pm.Normal('beta15',2,0.25),pm.Uniform('beta16',[-np.inf,np.inf])],
                          [pm.Uniform('beta21',[-np.inf, np.inf]),pm.Uniform('beta22',[-np.inf, np.inf]),
                           pm.Normal('beta23',6,2),pm.Normal('beta24',8,np.sqrt(3)),
                           pm.Normal('beta25',3,1),pm.Normal('beta26',5,10)],
                         [pm.Normal('beta31',9,1),pm.Uniform('beta32',[-np.inf, np.inf]),
                          pm.Normal('beta33',12,0.9),pm.Uniform('beta34',[-np.inf, np.inf]),
                          pm.Normal('beta35',12,0.9),pm.Uniform('beta36',[-np.inf, np.inf])]])
        Slope_Temp = Slope @ x
        likelihood = pm.Normal("y", mu=Intercept + np.sum(Slope_Temp,axis=1), sigma=Sigma, observed=y)
        idata = sample(3000)

        TypeError                                 Traceback (most recent call last)
        Cell In[98], line 24
             14 Slope = np.array([[pm.Normal('beta11',10,3),pm.Normal('beta12',1,4),
             15                    pm.Uniform('beta13',[-np.inf,np.inf]),pm.Normal('beta14',7,1),
             16                    pm.Normal('beta15',2,0.25),pm.Uniform('beta16',[-np.inf,np.inf])],
           (...)
             21                   pm.Normal('beta33',12,0.9),pm.Uniform('beta34',[-np.inf, np.inf]),
             22                   pm.Normal('beta35',12,0.9),pm.Uniform('beta36',[-np.inf, np.inf])]])
             23 Slope_Temp = Slope @ x
        ---> 24 likelihood = pm.Normal("y", mu=Intercept + np.sum(Slope_Temp,axis=1), sigma=Sigma, observed=y)
             25 idata = sample(3000)

            TypeError: Variables do not support boolean operations.

Could you recommend examples for structured multiple regression ?..

10

I don’t really know what “structured multiple regression” means to be honest, but from your code it seems you want to do several regression simultaneously, which is no problem. Here’s a similar example, where a user wanted to run a collection of independent regressions all in one go. In your case, you want to run 3 regressions together on the same MCMC chain, which is no problem. There are just some rules to respect:

  1. Don’t mix numpy and pytensor. As I already stated earlier in the thread, you cannot use numpy functions on PyMC random variables, because they are symbolic pytensor variables, and numpy does not know how to handle these. If you want to group together variables, you should generate them together using the size keyword, pass sequences to the distribution prarameters (from which PyMC will use to infer the size), or use pm.math.stack/pm.math.concatenate (as the case may be).
  2. The batch dimension always needs to come first. There’s a demonstration of this here. It will be important for you to familiarize yourself with how dimensions work to do this kind of problem.

Your specific error comes from breaking rule (1) above (mixing pytensor and numpy). When you do this you will get nonsense error like this one, that result from unexpected interactions between the two packages. Never mix them. But there are other problems in your code:

  1. y has shape (3, 1000). Remember, the batch dimension needs to come first, so we will need to transpose it.
  2. You cannot specify pm.Uniform over \mathbb R. That’s called pm.Flat.
  3. You are taking the sum of Slope_Temp over the batch dimension, so the quantity Intercept + Slope_Temp.sum(axis=1) will be only 3 numbers, because Intercept is shape (3,), x is shape (6, 1000), and Slope is shape (3, 6), so (3,) + ((3,6) @ (6, 1000)).sum(axis=1) = (3,) + (3, 1000).sum(axis=1) = (3,) + (3,) = (3,). meaning that all your regressions will be flat lines. I doubt this is what you intended.

Here is how I would start to re-write your code:

with pm.Model() as Structured_Model:
    Sigma = pm.HalfCauchy('Sigma', beta=[1, 6, 4])
    Intercept = pm.Normal("Intercept", mu=0, sigma=[9,4,5])
    
    Slope = pm.math.stack([
        pm.math.stack([pm.Normal('beta11',10,3),
                       pm.Normal('beta12',1,4), 
                       pm.Flat('beta13'),
                       pm.Normal('beta14',7,1),
                       pm.Normal('beta15',2,0.25),
                       pm.Flat('beta16')]),
        pm.math.stack([pm.Flat('beta21'),
                       pm.Flat('beta22'),
                       pm.Normal('beta23',6,2),
                       pm.Normal('beta24',8,np.sqrt(3)),
                       pm.Normal('beta25',3,1),
                       pm.Normal('beta26',5,10)]),
        pm.math.stack([pm.Normal('beta31',9,1),
                       pm.Flat('beta32'),
                       pm.Normal('beta33',12,0.9),
                       pm.Flat('beta34'),
                       pm.Normal('beta35',12,0.9),
                       pm.Flat('beta36')])])

From here, fix the remaining shape and code errors using the above rules.

1 Like
11

How kindly!!

Sincerely thank you.

12

Hello,I am working on a similar problem and what I want to know is how to specify the likelihood probability term, likelihood = pm.Normal("y", mu=Intercept + np.sum(Slope_Temp,axis=1), sigma=Sigma, observed=y), is this the correct ?

13

I discussed the shape logic in point (3) in the above post. It was wrong in the original poster’s model. I suggest you check the shapes of your specific problem and make sure they come out making sense.

1 Like