# Bayesian finite mixture model in stan
**URL:** https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667
**Category:** Modeling
**Created:** [August 24, 2017, 5:54am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667 "2017-08-24T05:54:49Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [August 24, 2017, 5:54am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/1 "2017-08-24T05:54:49Z")
Dear fellow stan users:
My model is a 2-component mixture of normal with a constraint such that the overall mean is modelled by a cyclic b-spline.
I am trying to model hour of day temporal effect where my y\_it (i=1:13, t=1:24, in other words, there are 13 IID observations per hour of day) is simulated from
```
##simulate a single day total PNC
TT = rep(1:24)
T = length(TT)
N = 13
##first mixture component hourly temporal effect
sinT = 0.5*sin(2*pi*TT/24)
##second mixture component hourly temporal effect
cosT = 0.5*cos(2*pi*TT/24)+10
mu <-matrix(c(sinT, cosT),T,2);
sigma <- c(0.1);
set.seed(100)
lambda = runif(1,0,1) ## when simulating the data, mixture weights are not time-dependent but when modelling the
##dynamic processes, we allow mixture weights to be time dependent to capture the correlation between observations.
##simulate component indicator
z <-matrix(0, 13, T)
for(i in 1:T){
z[,i] <- rbinom(N, 1, lambda) + 1;
}
y <- matrix(0, N, T)
set.seed(100)
for(i in 1:T){
y[,i] <- mu[i,z[,i]] + rnorm(N,0, sigma)
}
```
I would like to estimate time-dependent mixing proportions theta (a vector of 24 elements), time-dependent location mu (a vector of 24 elements ) and a global sigma.
Specifically, I want to use a cyclic B-spline to capture the hourly temporal effect and I rather than modelling mixture component as 2 cyclic B-spline, I model the overall mean
that is theta[t] \* mu[t,1] +(1-theta[t]) \* mu[t,2] = f(x\_t) where f(x\_t) = B\_T (cyclic bspline basis with 6 bases) \* theta\_cyclic (its corresponding spline coefficients to be estimated as well).
In addition, I impose a second-order penalty prior on theta\_cyclic to ensure that the spline fit is not too wiggly.
With all the above specification, my model cannot run and it gives the following error message:
Initialization between (-2, 2) failed after 100 attempts.
Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
I am assuming it has something to do with my highly constrained parameters or there are something I overlook?
Any suggestions would be greatly appreciate it.
Thanks, in advance[data.r](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/3cfcc6fb1bb0beef49a1fb754b9571f59583e669.r) (1.5 KB)
[classical\_gauss\_mix.stan](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/d599d0470fc8ec3234ac48eaa842aa354debb841.stan) (3.8 KB)
[bspline basis construction function.r](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/37b3503958c33b3823da69a8e1cc86b807ffabe0.r) (3.8 KB)
last file contains the functions for constructing cyclic bspline basis functions and 2nd order penalty prior for theta\_cyclic parameters.
[edit: escaped code blocks]
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [August 28, 2017, 8:17pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/2 "2017-08-28T20:17:16Z")
> [@Jaslene\_Lin](#):
>
> With all the above specification, my model cannot run and it gives the following error message:
>
> Initialization between (-2, 2) failed after 100 attempts.
>
> Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
>
> I am assuming it has something to do with my highly constrained parameters or there are something I overlook?
This almost always arises for one of two reasons:
1. missing constraint on a parameter—every legal value of the parameters matching the constraints should lead to a finite log density (i.e., non-zero density)
2. the problem is so tightly constrained that random inits take you so far into the tail that there are numerical issues causing rejections.
You should be seeing error messages about which things failed, which should give you some hints. If you’re using a pre-2.16 Stan, you should upgrade because we fixed the warning message output during initialization.
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [August 28, 2017, 8:25pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/3 "2017-08-28T20:25:20Z")
You can also save yourself a whole bunch of cutting and pasting by just declaring
```
ordered[2] mu_order1; //ordered mean
ordered[2] mu_order2;
ordered[2] mu_order3;
ordered[2] mu_order4;
ordered[2] mu_order5;
ordered[2] mu_order6;
ordered[2] mu_order7;
ordered[2] mu_order8;
ordered[2] mu_order9;
ordered[2] mu_order10;
ordered[2] mu_order11;
ordered[2] mu_order12;
ordered[2] mu_order13;
ordered[2] mu_order14;
ordered[2] mu_order15;
ordered[2] mu_order16;
ordered[2] mu_order17;
ordered[2] mu_order18;
ordered[2] mu_order19;
ordered[2] mu_order20;
ordered[2] mu_order21;
ordered[2] mu_order22;
ordered[2] mu_order23;
ordered[2] mu_order24;
```
as
```
ordered[2] mu_order[24];
```
and then in transformed parameters, if you really need a 24x2 matrix for something, doing
```
matrxi[24, 2] mu;
for (n in 1:24) mu[n] = mu_ordered[n]';
```
This
```
for(t in 1:T){
mu[t,1] ~ normal(0,2);
mu[t,2] ~normal(10,2);
}
```
can be just:
```
mu[, 1] ~ normal(0, 2);
mu[, 2] ~ normal(10, 2);
```
You should also check that `tau_y * lambda_temp * K_T` is always positive definite.
---
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [September 2, 2017, 12:24am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/4 "2017-09-02T00:24:38Z")
Hi Bob, thanks for your reply. I have fixed some of the errors in my model and updated the latest version and now I see the error message. The error message indicates that my mixing weights lambda\_t is not between (0,1). However, given the error message, I am still unable to get my model running.
Essentially, my model has two parts.
one is the 2-component mixture per time period. At a given point in time, the overall mean of y\_t\_mix = lambda\_t \* mu1\_t + (1-lambda\_t) \*mu2\_t;
the second part is the cyclic b-spline where at a given point in time, the overall mean of y\_t\_spline = cyclic B-spline basis (constructed outside and input as data) %\*% theta\_cyclic (which is the coefficients to be estimated).
To link this 2 parts, I impose y\_t\_mix = y\_t\_spline and derive an equation on lambda\_t.
lambda\_t = (y\_t\_spline - mu2\_t) / (mu1\_t - mu2\_t ) and the model likelihood is on the mixture model.
To figure out which part of my model goes wrong, I fit a model on bspline part only (where the likelihood function is on the spline part, hence response\_t ~ normal( y\_t\_spline, sigma) but estimate lambda\_t as above, importing simulated mu1\_t and mu2\_t and places no prior distribution or any constraints on lambda\_t ( it should have a constraint between (0,1)). In this case, lambda\_t can be estimated quite accurately as what I simulated. and the y\_mix\_spline also recovers my simulated data.
However, once I add on the (0,1) constraint and a customised prior on lambda\_t , the model fails.
The customised prior is a penalised prior which penalises a large change in the mixing weights over time.
I am wondering are there any things that I overlook or simply my model is too restrictive to be modelled by stan?
Thanks in advance for any suggestions.
Jaslene
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [September 3, 2017, 12:27pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/5 "2017-09-03T12:27:45Z")
What guarantees that
```
lambda_t = (y_t_spline - mu2_t) / (mu1_t - mu2_t )`
```
is between zero and one? For a start, you are going to constrain
```
mu1_t - mu2_t > 0
y_t_spline > mu2_t
y_t_spline - mu2_t < mu1_t - mu_2_t
```
So you need to declare all these parameters with constraints that make that true.
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [September 3, 2017, 12:28pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/6 "2017-09-03T12:28:18Z")
The alternative is to make some predictor unconstrained then apply a logit link (that is, apply `inv_logit(x)` to some unconstrained value `x`, meaning you model the log odds directly rather than the probability).
---
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [September 4, 2017, 12:10am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/7 "2017-09-04T00:10:28Z")
Hi Bob:
thanks for your suggestion. I used the second option and now it runs.
Kindest regards,
Jaslene
---
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [September 14, 2017, 5:51am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/8 "2017-09-14T05:51:01Z")
Hi Bob and others stan users:
My model is currently running (at least not experience any numerical issues) however, still the outputs are not entirely as expected.
to recap, the model has two structures.
firstly, a 2-component mixture Gaussian model lambda\_t \* N(mu1, sigma1) + (1-lambda\_t) \* N(mu2, sigma2) where the mixing weights are time-dependent while the component locations and scales are the same for all time periods. My data y\_(it) are simulated from this mixture model. Currently, y is a matrix of dimension 20\*24, in other words, there are 20 IID observations per hour of the day.
The second part of my model is that I enforce the mean of the weighted mixtures of Gaussian to be modelled by a cyclic bspline model f(x\_t) = B\_T ( that is a cyclic bspline basis functions) \* beta ( bspline coefficients). beta is a vector of 6 elements as I predetermine the number of basis in my cyclic bspline to be 6.
The enforced relationship is
B\_T \* beta = f(x\_t) = lambda\_t \*mu1 + (1-lambda\_t) \*mu2.
However, I cannot let f(x\_t) equal to B\_T \* beta and at the same time equal to lambda\_t \*mu1 + (1-lambda\_t) \*mu2 and thus I rewrite the above relationship in terms of a function of lambda\_t
for(t in 1:24){
lambda\_t = (f(x\_t)-mu2) / (mu1-mu2) }
Also I enforce the variance of the weighted mixtures of Gaussian in the model to be
for(t in 1:24){
y\_variance[t] = lambda[t] \*(pow(sigma[1],2)+pow(mu[1],2)) + (1-lambda[t]) \*(pow(sigma[2],2)+
pow(mu[2],2)) - pow((lambda[t]\*mu[1]+(1-lambda[t])\*mu[2]),2);
}
I assign uninformative priors on lambda\_t. mu1, mu2, sigma1 and sigma2 and beta.
for the model likelihood, I use the log\_mix
for(t in 1:24){
for(i in 1:N){
target += log\_mix(theta[t], normal\_lpdf(y[i,t] | mu[1], sigma[1]),
normal\_lpdf(y[i,t] | mu[2], sigma[2]));
```
}
```
}
I would expect to see the posterior samples of f(x\_t) to be exactly the same as samples constructed from lambda\_t\*mu1 + (1-lambda\_t) \*mu2.
However, this is not the case, while the posterior mean of the mixture model follows closely with the true simulated data mean trend, the bspline model tends to be more variable than the true simulated data mean trend.
What I really want is the equation of f(x\_t) = lambda\_t\*mu1 + (1-lambda\_t) \*mu2 to hold while at the same time f(x\_t) is constructed by B\_T (input data ) \* beta
I think my model makes logical sense but somehow I probably coded wrongly so my model is not producing what I want.
Any suggestions or feedback are greatly appreciated.
Thanks in advance,
Jaslene
---
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [September 15, 2017, 5:08am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/9 "2017-09-15T05:08:54Z")
```
functions {
//Define log probability density function
real my_penalised_beta_lpdf( vector theta, real phi, int period) {
vector[period] loglikelihood;
for(t in 1:period){
loglikelihood[t] = (beta_lpdf(theta[t] | 1,1))
- ((1/phi)*sum((theta[2:period]-theta[1:(period-1)]) .* (theta[2:period]-theta[1:(period-1)])))/period-
((1/phi)*sum(((1-theta[2:period])-(1-theta[1:(period-1)])) .* ((1-theta[2:period])-(1-theta[1:(period-1)]))))/period ;
}
return
sum(loglikelihood);
}
}
data {
int tt;//total number of time periods =96, repeat of 4 days
int N;//obs per day=20
int num_basis; //currently 6
matrix[tt,num_basis+1] B_overall; //cyclic bspline for hourly temporal effect
matrix[N,tt] y;
matrix[num_basis,num_basis] K_T;
//real phi;
}
parameters{
vector[num_basis] theta_cyclic;
real lambda_temp;//penalty parameter for temporal effect
real phi;//penalty parameter for temporal effect
ordered[2] mu; //ordered mean
vector[2] sigma;//component sigma
}
transformed parameters{
vector[tt] theta;//matrix containing p_t for all N observations over time
vector [tt] y_mean;
real deltatheta_cyclic;
vector[num_basis+1] theta_overall;
vector[tt] y_variance;//also model on the overall variance of the bspline
deltatheta_cyclic = sum(B_overall[,2:(num_basis+1)] * theta_cyclic)/
sum(B_overall[,2:(num_basis+1)]);
theta_overall[2:(num_basis+1)] = theta_cyclic - deltatheta_cyclic;
theta_overall[1] = deltatheta_cyclic;
y_mean = (B_overall * theta_overall);
for(t in 1:tt){
theta[t] = inv_logit((y_mean[t]-mu[2])/(mu[1]-mu[2]));
}
for(t in 1:tt){
y_variance[t] = theta[t] *(pow(sigma[1],2)+pow(mu[1],2)) + (1-theta[t]) *(pow(sigma[2],2)+
pow(mu[2],2)) - pow((theta[t]*mu[1]+(1-theta[t])*mu[2]),2);
}
}
model{
int period;
matrix[num_basis+1, num_basis+1] A0;
period = tt;
sigma~cauchy(0,2);
A0=diag_matrix(to_vector(rep_array(0, num_basis+1)));
A0[1,1] = 0.1;
A0[2:(num_basis+1),2:(num_basis+1)]= lambda_temp * K_T;
lambda_temp~gamma(1,10); //in STAN, the parameterisation is shape and rate and thus the prior rate = 10
phi~gamma(1,10); //in STAN, the parameterisation is shape and rate and thus the prior rate = 10
//penalty prior on theta_overall
target += multi_normal_prec_lpdf (theta_overall | to_vector(rep_array(0, (num_basis+1))),A0);
//penalised prior for time-dependent weight
target += my_penalised_beta_lpdf( theta|phi, period);
//prior on mixture component location
to_row_vector(mu) ~normal(0,10);
//loglikelihood on the mixture model
for(t in 1:tt){
for(i in 1:N){
target += log_mix(theta[t], normal_lpdf(y[i,t] | mu[1], sigma[1]),
normal_lpdf(y[i,t] | mu[2], sigma[2]));
}
} //checked total log likelihood is calculated as expected
// for(t in 1:tt){
// for(i in 1:N){
//y[i,t] ~ normal(y_mean[t], sqrt(y_variance[t]));
// }
// } //checked total log likelihood is calculated as expected
}
generated quantities{
vector[tt] y_mix;
for(t in 1:tt){
y_mix[t] = theta[t] *mu[1]+(1-theta[t])*mu[2];
}
}
```
**The above is my stan code**
**Here is my simulated data**
```
##simulate an example where we have time-dependent cyclic mixing weights lambda
##time-independent component location and scale parameters
TT = rep(c(1:24),1)
tt = length(TT)
N = 20
##first mixture component hourly temporal effect
lambda = 0.45*sin(2*pi*TT/24)+.5
mu=c(5,10)
sigma=c(1,1.5)
set.seed(100)
res = matrix(0,N,tt)
z = matrix(0, N, tt)
set.seed(100)
for(i in 1:N)
{
for(t in 1:tt){
z[i,t] <- sample(1:2, size = 1, prob = c(lambda[t], 1-lambda[t]), replace = T)
res[i,t] = mu[z[i,t]] + sigma[z[i,t]] * rnorm(1)
}
}
num_basis=6
B_T<-bspline(TT,8,2,1,xl=min(TT),xr=max(TT)) ##basis function for f(hour_t)
K_T<-make.Crw2(6)+diag(num_basis)*0.0001 ##2nd order penalty matrix for f(hour_t)
```
[bspline basis construction function.r](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/37b3503958c33b3823da69a8e1cc86b807ffabe0.r) (3.8 KB)
**with the attached function to create cyclic bspline basis function**
[edit: escaped code]
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [September 15, 2017, 4:30pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/10 "2017-09-15T16:30:53Z")
I’m afraid these really complicated programs are difficult to debug. It doesn’t help that I’ve never understood spline models.
You may want to start with something simpler and try to build up to something this complex stepwise. It can be faster as when you introduce a wrong step, you’ll know where it came from.
I couldn’t quite tell from the description if you were testing posterior interval coverage. Just comparing the posterior means shouldn’t matter much.
Also, Michael Betancourt wrote a tutorial on mixture models and some of the problems with identifying them—it’s on our web site under users \>\> documentation \>\> case studies.
---
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [September 15, 2017, 9:31pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/11 "2017-09-15T21:31:28Z")
Hi Bob:
thanks for your reply. Without getting into the details of my model, I think what I want to achieve in this model is simply modelling y\_it coming from a mixture of 2 Normal distributions with mixing weights lambda\_t, component 1 parameters (mu1, sigma1) and component 2 parameters (mu2, sigma2). I have implemented this simple model in STAN and it works fine.
The only twist to the above model is that I would want the mean of the overall mixture distribution to equal to a function f(x\_t) modelled by bsplines, but we can just think of as a function of x\_t for now.
In other words, I want **f(x\_t) = lambda\_t\*mu1 + (1-lambda\_t) \*mu2 to hold** while **f(x\_t) = B\*beta**. B is just a data design matrix and beta is to be estimated.
I would want to write in the transformed parameters block something like
f(x\_t) = Bbeta= lambda\_t\*mu1 + (1-lambda\_t) \*mu2,
but it is not allowed and f(x\_t) will be overwritten.
Therefore, I write lambda\_t as a function of f(x\_t) and mu1, mu2 such that the above relationship is satisfied.
Finally, I assign priors on lambda\_t, mu1, mu2 and beta.
In my posterior samples, I would expect to see B beta= lambda\_t\*mu1 + (1-lambda\_t) \*mu2 should hold, however, it does not. that is what bugs me and I am trying to resolve at the moment.
Thanks again for your time,
Jaslene
---
### Author:  [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)
#### Post date: [September 17, 2017, 4:02pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/12 "2017-09-17T16:02:27Z")
> [@Jaslene\_Lin](#):
>
> The only twist to the above model is that I would want the mean of the overall mixture distribution to equal to a function f(x\_t) modelled by bsplines
It’s usually much simpler and better to stick to simple generative processes—without knowing how your data comes about, it’s hard to say more.
If you really want to go the constrained route, the easiest way to do this is to declare `mu_1` as a parameter and then solve for `mu_2` given the spline value. That gives you the one degree of freedom (relative to the spline) that you need. Then you can put a prior on `mu_1`, but if you want to put a prior on `mu_2` you’ll need the Jacobian for the constraint solution. Instead, you should be putting priors on the parameters for the spline, not on `mu_2`, which is just a constrained solution, not a free parameter.
---
### Author:  [@stijn](https://discourse.mc-stan.org/u/stijn)
#### Post date: [September 18, 2017, 6:57am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/13 "2017-09-18T06:57:10Z")
> [@Jaslene\_Lin](#):
>
> In other words, I want f(x\_t) = lambda\_t\*mu1 + (1-lambda\_t) \*mu2 to hold while f(x\_t) = B\*beta.
I think you can do what you want if you model y as a mixture of N(B \* beta + nu1, sigma1) and (B \* Beta + nu2, sigma2), where nu2 = - lambda\_t / (1 - lambda\_t).
mu1 = B\*beta + nu1, mu2 = B\*beta + nu2
---
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [September 18, 2017, 8:43pm UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/14 "2017-09-18T20:43:08Z")
Hi stijn:
Thanks for your suggestion, I will give it a try.
---
### Author:  [@Jaslene\_Lin](https://discourse.mc-stan.org/u/Jaslene_Lin)
#### Post date: [September 19, 2017, 7:11am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/15 "2017-09-19T07:11:02Z")
I have implemented what you have suggested however this does not seem to work. however, just to be assured that I did understand what you have suggested.
In this case I need f(x\_t) = lambda\_t\*mu\_1t+ (1-lambda\_t) _mu\_2t to hold while f(x\_t) = B_beta.
mixing weights lambda\_t. component locations mu\_1t, mu\_2t are time dependent, while beta is a vector of size (number of cyclic bspline basis function \* 1) and B is a cyclic bspline basis function of size ( t \* num\_basis).
Following your suggestion, I declare in the parameter block
vector[num\_basis] theta\_cyclic;
vector\ [tt] lambda;//matrix containing p\_t for all N observations over time
vector\[2] sigma;//component sigma
To ensure identifiability, I force mu\_2t \> mu\_1t is valid for all time periods, which means in the model block I have
for(t in 1:tt){
for(i in 1:N){
target += log\_mix(lambda[t], normal\_lpdf(y[i,t] | (y\_mean[t]-1), sigma[1]),
normal\_lpdf(y[i,t] | y\_mean[t]+nu2[t], sigma[2]));
```
}
```
} //checked total log likelihood is calculated as expected
where y\_mean = B \* beta and this is equivalent to your suggestion of modelling the mixture
lambda\_t \* N(B \* beta -1 , sigma1) + (1-lambda\_t) \* (B \* Beta + nu2, sigma2), where nu2 = lambda\_t / (1 - lambda\_t).
Also to model mu\_1t and mu\_2t, since one of them now will not be a free parameter. I think I can write it deterministically.
I have attempted to
declare mu\_1t in the transformed parameter block as mu\_1t = (B \* beta) [t] -1 .
since I have assigned priors on beta then no further priors are assigned.
I do not know the appropriate place where I should declare mu\_2t and its relationship to mu\_1t.
right now, it has been declared in the generated quantities block, I declare mu\_2t = mu\_1t + 1+ nu2[t].
However, my model does not recover the “true” values of mu\_1t and mu\_2t I simulated, the model however, roughly recover the mean trend modelled by y\_mean = B \* beta.
I have read Stan website case study regarding the fitting a mixture model. If I remove the bspline part fits the data only with the mixture model then with the ordered exchangeable prior on mu\_1t and mu\_2t then the model can correctly estimate and recover the simulated values. However, it does not seem to work in the case where I impose a constraint that the overall mean of the mixture model should equal to a smooth function modelled by the bspline. I am still a bit confused and any ideas are greatly appreciated.
---
### Author:  [@stijn](https://discourse.mc-stan.org/u/stijn)
#### Post date: [September 19, 2017, 7:54am UTC](https://discourse.mc-stan.org/t/bayesian-finite-mixture-model-in-stan/1667/16 "2017-09-19T07:54:51Z")
What I meant was the following but there was a typo in my post.
```
parameter{
vector[tt] nu1;
vector [tt] lambda;
...
}
transformed parameter{
vector[tt] nu2;
nu2 = - nu1 * lambda/(1 - lambda);
...
}
model{
...
for(t in 1:tt){
for(i in 1:N){
target += log_mix(lambda[t], normal_lpdf(y[i,t] | y_mean[t] + nu1[t]), sigma[1]),
normal_lpdf(y[i,t] | y_mean[t] + nu2[t], sigma[2]));
}
}
}
```
The idea is that with the restriction on nu2, the expected value of lambda \* nu1 + (1 - lambda) \* nu2 = 0. So the expected value for y[i,t] is y\_mean[t].