# pm.Potential and confliction between Theano and Numpy

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?

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

Anyone who could be helpful will be appreciated, thanks a lot!

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:

``````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()
``````
3 Likes

Regards, Daniel.

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.

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.

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.

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.

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.

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.

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:

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

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