# Modeling of multiple regression model with array form

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:

(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

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

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.

``````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.

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

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â€¦

Itâ€™s not bother, youâ€™re helping me avoid work

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).

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)]

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

Dear Jesse.

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 ?..

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

How kindly!!

Sincerely thank you.