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