# Two component mixture model, strange result from stan

**URL:** <https://discourse.mc-stan.org/t/two-component-mixture-model-strange-result-from-stan/8722>\
**Category:** Modeling\
**Created:** [May 5, 2019, 12:29am UTC](https://discourse.mc-stan.org/t/two-component-mixture-model-strange-result-from-stan/8722 "2019-05-05T00:29:17Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jg\_liao](https://avatars.discourse-cdn.com/v4/letter/j/13edae/32.png) [@jg\_liao](https://discourse.mc-stan.org/u/jg_liao)\
**Post date:** [May 5, 2019, 12:29am UTC](https://discourse.mc-stan.org/t/two-component-mixture-model-strange-result-from-stan/8722/1 "2019-05-05T00:29:17Z")

</div>

I have been struggling with the following code for a while and decided to seek help from experts here.

I am fitting a two component mixture model, one component normal and one component laplace. When I fit data that is generated from normal distribution with a large sample size of N=1000 or N=5000, I would expect the weight on the normal component to be close to 1. However, stan result always gives a weight close to 0.65. The same is true if the data is generated from a laplace distribution.

The code and result are

liao1 = "data {  
int N; // number of observations  
real theta[N];  
}

parameters {  
real\<lower=0.0, upper= 20/sqrt(2)\> tau1; //sqrt(2)\*tau1 is the standard deviation of laplace distribution  
real\<lower=0.0, upper=20.0\> tau2; //tau2 is the standard deviation of normal distribution  
real\<lower=0.0, upper=1.0\> lambda;  
}

model {

real term1;  
real term2;

tau1 ~ uniform(0.0, 20/sqrt(2.));  
tau2 ~ uniform(0.0, 20.0);

lambda ~ uniform(0.0, 1); //lambda\<0.0 and lambda\>=0 represents binary choice of laplace distribution and normal distribution

target += log\_mix(lambda, double\_exponential\_lpdf(theta |0.0, tau1), normal\_lpdf(theta | 0.0, tau2));  
}  
"

N = 1000  
theta = rnorm(N, 0, 2)

my.data = list(N=N, theta=theta)

stan(model\_code = liao1, data=my.data, iter=10000, chains=1)

Inference for Stan model: ae78845bd99245d1341f7cb022b70ce0.  
1 chains, each with iter=10000; warmup=5000; thin=1;  
post-warmup draws per chain=5000, total post-warmup draws=5000.

```
       mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat

```

tau1 7.00 0.06 4.04 0.36 3.52 6.95 10.48 13.75 5000 1  
tau2 2.01 0.00 0.05 1.93 1.98 2.01 2.04 2.10 4437 1  
lambda 0.33 0.00 0.24 0.01 0.12 0.28 0.50 0.85 4709 1  
lp\_\_ -2118.42 0.03 1.38 -2121.92 -2119.08 -2118.05 -2117.39 -2116.82 2027 1

Samples were drawn using NUTS(diag\_e) at Sat May 04 20:25:49 2019.  
For each parameter, n\_eff is a crude measure of effective sample size,  
and Rhat is the potential scale reduction factor on split chains (at  
convergence, Rhat=1).

Thanks in advance.

Jason Liao

---

<div class="post-metadata">

**Author:** ![Stephen\_Martin](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stephen_martin/32/7879_2.png) [@Stephen\_Martin](https://discourse.mc-stan.org/u/Stephen_Martin)\
**Post date:** [May 7, 2019, 12:17am UTC](https://discourse.mc-stan.org/t/two-component-mixture-model-strange-result-from-stan/8722/2 "2019-05-07T00:17:42Z")

</div>

I’m fairly certain log\_mix does not do what you think it is doing. Iirc, in the stan manual in the section on mixtures/clustering models, it talks about this.

You need to loop over the thetas.  
Currently, you are estimating:

p(\theta) = (1-\lambda)\prod\_{i=1}^Np(\theta|\tau\_1) + \lambda\prod\_{i=1}^Np(\theta|\tau\_2)

What you probably want is:

p(\theta) = \prod\_{i=1}^N \left[(1-\lambda)p(\theta\_i|\tau\_1) + \lambda p(\theta\_i|\tau\_2)\right]

Which means you want to loop over each theta and call log\_mix. E.g.,

```
for(n in 1:N){
  target += log_mix(lambda, double_exponential_lpdf(theta[n]|0, tau1), normal_lpdf(theta[n]|0, tau2));
}

```

Someone correct me if I’m wrong.

---

<div class="post-metadata">

**Author:** ![jg\_liao](https://avatars.discourse-cdn.com/v4/letter/j/13edae/32.png) [@jg\_liao](https://discourse.mc-stan.org/u/jg_liao)\
**Post date:** [May 7, 2019, 3:10am UTC](https://discourse.mc-stan.org/t/two-component-mixture-model-strange-result-from-stan/8722/3 "2019-05-07T03:10:35Z")

</div>

Stephen, this is amazing. I changed the code to what you suggested and got the answer I desperately wanted.

liao1 = ’

data {  
int N; // number of observations  
real theta[N];  
}

parameters {  
real\<lower=0.0, upper= 20/sqrt(2)\> tau1; //sqrt(2)\*tau1 is the standard deviation of laplace distribution  
real\<lower=0.0, upper=20.0\> tau2; //tau2 is the standard deviation of normal distribution  
real\<lower=0.0, upper=1.0\> lambda;  
}

model {

real term1;  
real term2;

tau1 ~ uniform(0.0, 20/sqrt(2.));  
tau2 ~ uniform(0.0, 20.0);

lambda ~ uniform(0.0, 1); //lambda\<0.0 and lambda\>=0 represents binary choice of laplace distribution and normal distribution

for(n in 1:N){  
target += log\_mix(lambda, double\_exponential\_lpdf(theta[n]|0, tau1), normal\_lpdf(theta[n]|0, tau2));  
}  
}  
’  
N = 1000  
theta = rnorm(N, 0, 2)

#theta = rexp(N)  
#theta = ifelse(rnorm(N)\>0, theta, -theta)

my.data = list(N=N, theta=theta)

stan(model\_code = liao1, data=my.data, iter=10000, chains=1)

```
     mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff

```

tau1 1.84 0.09 2.39 0.13 0.70 1.15 1.69 10.57 655  
tau2 1.95 0.00 0.05 1.86 1.92 1.95 1.99 2.06 1236  
lambda 0.05 0.00 0.05 0.00 0.01 0.04 0.08 0.18 767  
lp\_\_ -2082.41 0.06 1.58 -2086.27 -2083.28 -2082.10 -2081.16 -2080.35 678

One million thanks.

Jason