# pm.Potential and confliction between Theano and Numpy

**URL:** <https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499>\
**Category:** v5\
**Tags:** theano, gaussian\_process\
**Created:** [February 28, 2023, 3:47pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499 "2023-02-28T15:47:46Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![Daniel](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/daniel/32/6184_2.png) [@Daniel](https://discourse.pymc.io/u/Daniel)\
**Post date:** [February 28, 2023, 3:47pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/1 "2023-02-28T15:47:46Z")

</div>

Dear all, I am using the self-defined likelihood function together with pm. Potential.  
However, I met this problem: setting an array element with a sequence. I attached the code here.  
Hereafter, Data is a 1D numpy array, meaning measured data;  
Datac is a constant vector with the same shape of Data, e.g. Datac=np.random.rand(23).

The model has three parameters to be determined, theta=[mu,sig and delta].  
The prior distribution is just uniform for all components of theta.

I suppose it is because I am not using Theano, but how could this problem be solved without theano?

```python
Datac=np.random.rand(23)
import numpy as np
import pymc as pm
import scipy.stats as st

def likli(mu,sig,delta,Data):
    D_len=Data.shape[0]
    l=np.ones(D_len)
    R=np.ones((D_len,D_len))
    for i in range(D_len): 
        for j in range(D_len):
            tau=np.abs(Datac[i]-Datac[j])
            #print(delta,tau)
            tmp=np.exp(-2*tau/delta)
            """
            try: 
                if len(tmp)>=1:
                    tmp=tmp[0]
            except:
                tmp=tmp
            """
            R[i,j]=tmp  
            """ERROR: setting an array element with a sequence."""
    C=sig**2*R  
    res=(Data-mu*l)
    p0=(np.pi*2)**(-D_len/2)*(np.abs(np.linalg.det(C)))**(-0.5)*np.exp( -0.5 * (((res.T)@np.linalg.inv(C)@res)))
    return np.log10(p0)

with pm.Model() as model:
    mu = pm.Uniform('mu', lower=5, upper=17)
    sig = pm.Uniform('sig',lower=0, upper=4) 
    delta = pm.Uniform('delta',lower=0.5,upper=3)

    like = pm.Potential('like',likli(mu,sig,delta,Data))
    #like = pm.DensityDist('like', likli_dst(mu,sig,delta), observed=[Data])
    trace = pm.sample(500)

```

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/0/0006b4a67739aa9216184f51f2564dc06eb4ab03.png)  
Anyone who could be helpful will be appreciated, thanks a lot!

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [February 28, 2023, 4:30pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/2 "2023-02-28T16:30:48Z")

</div>

Hi Daniel!

I ran your code and didn’t hit the same error, so I’d double-check your inputs. That said, it seems like you are just using a multivariate normal likelihood, so you might be able to dispense with `pm.Potential` and instead just build the covariance matrix `C` inside the model block.

Consider this model:

```auto
import pymc as pm
import pytensor.tensor as pt

with pm.Model() as mod:
    data = pm.ConstantData('data', np.random.normal(size=(36,)))
    datac = pm.ConstantData('datac', np.random.normal(size=(36,)))
    
    mu = pm.Uniform('mu', lower=5, upper=17)
    sig = pm.Uniform('sig',lower=0, upper=4)
    delta = pm.Uniform('delta',lower=0.5,upper=3)
    
    # Broadcasting subtraction across column-datac and row-datac is equivalent to your double loop
    R = pt.exp(-2 * pt.abs(datac[:, None] - datac[None]) / delta)
    
    #pt.linalg.inv doesn't have a defined gradient, but pt.linalg.solve does
    cov = pm.Deterministic('cov', pt.linalg.solve(sig **2 * R, pt.eye(36), 
                                                  assume_a='pos', lower=False, check_finite=False))
    
    obs = pm.MvNormal('obs', mu=mu, cov=cov, observed=data)
    
    idata = pm.sample()

```

---

<div class="post-metadata">

**Author:** ![Daniel](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/daniel/32/6184_2.png) [@Daniel](https://discourse.pymc.io/u/Daniel)\
**Post date:** [March 1, 2023, 3:45pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/3 "2023-03-01T15:45:23Z")

</div>

Your answer is much more refined. Thanks a lot, Jesse.  
Regards, Daniel.

---

<div class="post-metadata">

**Author:** ![Daniel](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/daniel/32/6184_2.png) [@Daniel](https://discourse.pymc.io/u/Daniel)\
**Post date:** [March 1, 2023, 5:26pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/4 "2023-03-01T17:26:26Z")

</div>

Hi, Jesse. When executing the line of

R = pt.exp(-2 \* pt.abs(datac[:, None] - datac[None,:]) / delta)

I met the problem:  
NotImplementedError: Cannot convert Elemwise{sub,no\_inplace}.0 to a tensor variable.

 ![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/1/12c4965db7f43853ed028887b40d2ca5dcc10854.png)  
I run the code on Google Colab, and the versions of related packages are:  
#pymc 4.1.4  
#pytensor 2.10.1.  
#Numpy 1.20.3  
If possible, can you please show some hints for solving this?  
Thanks a lot.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [March 1, 2023, 5:33pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/5 "2023-03-01T17:33:50Z")

</div>

You can’t use pytensor with pymc v4. Prior to version 5, the computational backend was `aesara`. You just need to change the import from `import pytensor.tensor as pt` to `import aesara.tensor as at`, and find-replace `pt.` to `at.`

Also restart your environment and don’t install pytensor along side aesara (which comes pre-installed on the colab) because having both installed can cause unexpected things to happen.

---

<div class="post-metadata">

**Author:** ![Daniel](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/daniel/32/6184_2.png) [@Daniel](https://discourse.pymc.io/u/Daniel)\
**Post date:** [March 1, 2023, 8:26pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/6 "2023-03-01T20:26:57Z")

</div>

Thanks a lot, Jesse. You are amazing. Every thing works fine, except warning:  
aesara/tensor/slinalg.py:454: LinAlgWarning: Ill-conditioned matrix (rcond=7.00257e-17): result may not be accurate  
I think it’s the problem about mine cov matrix.  
Thanks again for your kind help.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [March 1, 2023, 8:36pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/7 "2023-03-01T20:36:16Z")

</div>

Inverting random matrices is non-trivial. You need to make sure that R is always full rank and, in your case, positive definite. I don’t know enough about your specific model to give recommendations here. You could try `at.linalg.pinv` if you think it’s completely impossible to impose more structure on R.

---

<div class="post-metadata">

**Author:** ![Daniel](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/daniel/32/6184_2.png) [@Daniel](https://discourse.pymc.io/u/Daniel)\
**Post date:** [March 2, 2023, 8:45am UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/8 "2023-03-02T08:45:55Z")

</div>

Thanks a lot for your kind explanation, Jesse. The problem is a Gaussian random field (people assumed it to be), which I think also could be treated as a 1D random signal. I don’t know why the actually measured data cannot guarantee an invertible R, maybe it’s because the measured data cannot be treated as a 1D gaussian field. Since autocorrelation density can be converted from power spectral density by the Fourie transform, I think maybe it’s better for me to generate the random field directly from power spectral density. As far as I know, some python packages exist for doing this.

---

<div class="post-metadata">

**Author:** ![jessegrabowski](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jessegrabowski/32/5010_2.png) [@jessegrabowski](https://discourse.pymc.io/u/jessegrabowski)\
**Post date:** [March 2, 2023, 12:24pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/9 "2023-03-02T12:24:16Z")

</div>

It might also be because the estimated `sig` is very small, so you multiply the whole matrix by almost-zero. You could try a prior that rules out zero, like Gamma(2,5)?

Another option for diagnosis is to record the rank of `R` in a deterministic variable by checking the number of singular values greater than zero:

```auto
S = at.linalg.svd(R, compute_uv=False)
pm.Deterministic('rank_R', (S > 1e-8).sum())

```

Then you can do some statistical analysis on the trace to see if low-rank samples are associated with values of other parameters.

---

<div class="post-metadata">

**Author:** ![Daniel](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/daniel/32/6184_2.png) [@Daniel](https://discourse.pymc.io/u/Daniel)\
**Post date:** [March 2, 2023, 8:01pm UTC](https://discourse.pymc.io/t/pm-potential-and-confliction-between-theano-and-numpy/11499/10 "2023-03-02T20:01:04Z")

</div>

Thanks for your kind help, Jesse. I will have a thorough check of my problem.