# Understanding log likelihood of sums

**URL:** <https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060>\
**Category:** Modeling\
**Tags:** fitting-issues\
**Created:** [September 29, 2017, 3:40pm UTC](https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060 "2017-09-29T15:40:08Z")\
**Posts on this page:** 5\
**Page:** 2

<div class="post-metadata">

**Author:** ![finnholger](https://avatars.discourse-cdn.com/v4/letter/f/13edae/32.png) [@finnholger](https://discourse.mc-stan.org/u/finnholger)\
**Post date:** [October 2, 2017, 3:11pm UTC](https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060/21 "2017-10-02T15:11:18Z")

</div>

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:

![fig](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/6e7f9861fd1594cd704b1625da21fdf44137f032.png)

---

<div class="post-metadata">

**Author:** ![dlakelan](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dlakelan/32/2219_2.png) [@dlakelan](https://discourse.mc-stan.org/u/dlakelan)\
**Post date:** [October 2, 2017, 6:38pm UTC](https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060/22 "2017-10-02T18:38:48Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![finnholger](https://avatars.discourse-cdn.com/v4/letter/f/13edae/32.png) [@finnholger](https://discourse.mc-stan.org/u/finnholger)\
**Post date:** [October 3, 2017, 12:48pm UTC](https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060/23 "2017-10-03T12:48:20Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![finnholger](https://avatars.discourse-cdn.com/v4/letter/f/13edae/32.png) [@finnholger](https://discourse.mc-stan.org/u/finnholger)\
**Post date:** [October 3, 2017, 1:04pm UTC](https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060/24 "2017-10-03T13:04:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dlakelan](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dlakelan/32/2219_2.png) [@dlakelan](https://discourse.mc-stan.org/u/dlakelan)\
**Post date:** [October 3, 2017, 2:47pm UTC](https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060/25 "2017-10-03T14:47:52Z")

</div>

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

[Previous page](https://discourse.mc-stan.org/t/understanding-log-likelihood-of-sums/2060.md?page=1)