How to do matrix multiplication and define likelihood by myself in PyMC


I have a model with 8 parameters and I can use this model to generate 5 distribution function. With 5 distribution function from data, I want to use pymc to constrain these 8 parameters. But this model is a little complicated because it includes matrix multiplication. And I faced some problems. Here is what I tried:

Firstly I draw prior for 8 parameters from uniform distribution.

x = pm.Uniform(“x”, lower=range[0], upper=range[1], shape=(1, 8))

Then I need to do some matrix multiplication, Y1)+Y2

where the shape for Y1 is (8,17), for Y2 is (1,17) and Y1, Y2 are numpy array.

I successfully done this and get a pytensor variable and its shape is (1,17)

Next I need to split matrix1 into two variables and their shapes are (1,2),(1,4),(1,4),(1,5),(1,2).

It seems like I cannot use method in numpy like matrix1[:2], matrix1[2:6], matrix1[6:10], matrix1[10:15], matrix1[15:17].

By using these 5 matrix, I can do another matrix multiplication independently and this time will give me 5 distribution function.

Recall that I have 5 distribution function from data. Now I need to define likelihood. I don’t think I should use pm.Normal(). In fact I want to define it as sum(-(Y_Data-Y)^2/Y_Data^2). But I don’t know what should I do. Maybe I need to use pm.Deterministic().

I appreciate if any one could help me

Hi Jinning,

You can indeed slice/index pytensor tensors just like numpy. Are you getting an error when you try to do this? Could you share the exact code you’ve tried?

You can add a arbitrary likelihood term to your model by using pm.Potential. The only downside is that pm.Potential terms are not taken into account when sampling posterior predictive distributions. Deterministic is just for signaling to PyMC that you would like an arbitrary computation to be sampled and stored in your idata. Using it won’t affect the log likelihood per se.

Hi Jesse,

I made a mistake when slice pytensor. Now it works. And I follow your suggestion using pm.Potential and everything is solved now.
Really thank you!

1 Like

Potential isn’t perfect (for the reasons noted). Looking at your likelihood, you could use a Normal with sigma = Y_data?