Still, this doesn’t seem to give the correct result. Here’s the code:
stan_model.stan
data {
int<lower=1> N;
real<lower=0> D;
real<lower=0> noise;
real s[N];
}
parameters {
real x1[N];
real x2[N];
real x3[N];
}
model {
// Intialise
for (n in 1:N) {
target += normal_lpdf(x1[n] | 0, D);
target += normal_lpdf(x2[n] | 0, D);
target += normal_lpdf(x3[n] | 0, D);
}
// Measurements
for (n in 1:N) {
s[n] ~ normal(x1[n] + x2[n] + x3[n], noise);
}
}
mcmc.py
import numpy as np
import matplotlib.pyplot as plt
import pickle
from pystan import StanModel
if True:
sm = StanModel(file='stan_model.stan')
with open('model.pkl', 'wb') as f:
pickle.dump(sm, f)
else:
sm = pickle.load(open('model.pkl', 'rb'))
N = 50
D = 5.
noise = 0.01
D_list = np.linspace(0.2*D, 1.5*D, 20)
mean_sampling = np.zeros_like(D_list)
x1 = D*np.random.randn(N)
x2 = D*np.random.randn(N)
x3 = D*np.random.randn(N)
s_true = x1 + x2 + x3
s_mes = s_true + np.random.randn(N)*noise
# MCMC
mean_sampling = np.zeros_like(D_list)
for i, D_guess in enumerate(D_list):
print i
data = dict(s=s_mes, N=N, D=D_guess, noise=noise)
fit = sm.sampling(data=data)
if i == 0:
M = np.max(fit['lp__'])
mean_sampling[i] = np.mean(np.exp(fit['lp__'] - M))
plt.plot(D_list, np.log(mean_sampling), '-rx')
plt.xlabel('D')
plt.ylabel('log P')
plt.savefig('fig.png')
plt.show()
Although the true value of D is D = 5, the maximum is found for lower:
ideally what you’d like is:
arg max over D: integrate(p(D|x1,x2,x3,s,noise) prior(x1,x2,x3 |s,noise) dx1dx2dx3)
that’s what you get by sampling and then pooling D values and then KDE estimate for the argmax
the quantity you are calculating is:
arg max D integrate(p(x1,x2,x3 | s,noise,D) dx1dx2dx3)
which doesn’t seem to be the same thing, and this difference is I think what makes it not work.
EDIT: hmm I threw this comment out there in between doing a few other things, but now I think I should be more careful:
integrate(p(x1,x2,x3 | s,noise,D)p(D | prior) dx1dx2dx3) is the posterior marginal density on D and you’re looking to maximize that. The only difference is that what you’re working with doesn’t have p(D|prior) so what prior did you use?
Also in your numerical graph, I’d use the M value from the 5th point or something if the log of your mean value is 300 that’s a really big number you’re taking the log of… double floats only go up to 10^308 or so for example
I used D ~ normal(0, 100). But I doubt that this can change anything over the range of D in [0, 7]…
I’ll try using M from the fifth point.
Using the fifth point still gives the same result…
Code:
import numpy as np
import matplotlib.pyplot as plt
import pickle
from pystan import StanModel
if False:
sm = StanModel(file='stan_model.stan')
with open('model.pkl', 'wb') as f:
pickle.dump(sm, f)
else:
sm = pickle.load(open('model.pkl', 'rb'))
N = 50
D = 5.
noise = 0.01
D_list = np.linspace(0.2*D, 1.5*D, 20)
mean_sampling = np.zeros_like(D_list)
x1 = D*np.random.randn(N)
x2 = D*np.random.randn(N)
x3 = D*np.random.randn(N)
s_true = x1 + x2 + x3
s_mes = s_true + np.random.randn(N)*noise
# MCMC
mean_sampling = np.zeros_like(D_list)
iD = list(enumerate(D_list))
iD = [iD[5]] + iD[:5] + iD[6:] # use index 5 to calc M
for i, D_guess in iD:
print i
data = dict(s=s_mes, N=N, D=D_guess, noise=noise)
fit = sm.sampling(data=data)
if i == iD[0][0]:
M = np.max(fit['lp__'])
mean_sampling[i] = np.mean(np.exp(fit['lp__'] - M))
plt.plot(D_list, np.log(mean_sampling), '-rx')
plt.xlabel('D')
plt.ylabel('log P')
plt.savefig('fig.png')
plt.show()
Ok, then I’m a bit confused why you would get a different answer from doing my method vs doing your method. Can you post a plot of the log of the posterior KDE after doing my method for comparison? Do they look very similar? When an objective function is very flat near the peak, and there’s noise, you can wind up with the actual value of the arg-max being kind of random. Perhaps the KDE just got lucky? Looking at the curve would help detect if we’re not noticing something seriously different about the calculation, or if we’re just talking about noisy estimates that might have gotten lucky