I am using pystan. Shockingly, I am finding a mismatch between log probability calculated by pystan and the same calculated by hand. Suppose I have 300 data points: x, y(x) and y error \sigma_y(x), all three vectors have length 300. There is no correlation whatsoever.
At any point log-probability should be:
ln_p=-0.5*N*\ln(2 \pi)-\sum_i^{N} \ln({\sigma_y}_i)-\chi^2/2, where N=300. I am finding that this value is not matching with the lp_
value given by the pystan fit! The mismatch is small but quite significant. I do not believe it is coming from truncation error. In pystan, I am using
target+=normal_lpdf(y[i]|F[i],yerror[i])
where y[i], F[i], yerror[i]
represent data, model and error at the ith data point.
Where am I doing the mistake? Has anyone by hand calculated the log-probability and matched with the pystan value ever?
Edit: Because of the mismatch, the maximum lp_
does not correspond to minimum \chi^2 where the latter is calculated by hand at each sampled point.
Edit: Working example
import matplotlib.pyplot as plt
import pystan
import numpy as np
def function(x,a,b):
return a*x+b
a,b=5.0,1.0
Nd=300
x=np.linspace(0.5,5.5,Nd)
fx=function(x,a,b)
sigfx=np.full(Nd,1.0); #print(sigfx)
y=np.random.normal(fx,sigfx)
plt.errorbar(x,y, yerr=sigfx,c="r")
plt.plot(x,fx,c="b")
model = """
data {
int<lower=0> N;
vector[N] x;
vector[N] y;
vector[N] yerror;
}
parameters {
real<lower=0.0, upper=20> A;
real <lower=0, upper=10> B;
//real sig;
}
/*model {
for(i in 1:N) {
real F;
F=A*x[i]+B;
target+=normal_lpdf(y[i]|F,yerror[i]);
}
}*/
model{
vector[N] F;
F=A*x+B;
target+=normal_lpdf(y|F,yerror);
//y~normal(F,yerror);
}
"""
# Put our data in a dictionary
data = {'N': len(x), 'x': x, 'y': y, 'yerror':sigfx}
# Compile the model
sm = pystan.StanModel(model_code=model)
# Train the model and generate samples
fit = sm.sampling(data=data, iter=2000, chains=4, warmup=400, thin=3, seed=101)
print(fit)
def get_fit_bf(fit):
ex=np.array(fit.extract(permuted=False));
ex_flat=np.concatenate(ex,axis=0)
idx=np.argmax(ex_flat[:,-1])
MLE_pars=ex_flat[idx,:]
return[MLE_pars]
def get_chisq(pars):
rec_fx=function(x,pars[0],pars[1])
chisq=np.sum((y-rec_fx)**2/sigfx**2)
return chisq
print("Finding best-fit from the lnp..")
MLE_pars= get_fit_bf(fit)
print(MLE_pars[0])
Am,Bm,lpm=MLE_pars[0]
chisq_pyst=get_chisq(MLE_pars[0])
print("Chi^2=",chisq_pyst)
lnp=-0.5*Nd*np.log(2.0*np.pi)-np.sum(np.log(sigfx))-chisq_pyst/2.0
print("lnp for these values should be: ",lnp)