Mixture model in stan

Modeling
1

Suppose \beta \sim (1-\pi)N(0,1) + \pi \delta_{0}(\beta), where \delta_{0}(\beta) denotes a point mass at 0, equally, p(\beta = 0) = 1.

My problem is , how can I write this model into stan?

I can use target function to achieve this, my code is as following,

data{
  int<lower=0> N; // number of obs
  int<lower=0> p; // number of predictors
  int<lower=0> K;
  int<lower=0,upper=1> y[N,K];
  matrix[N, p] X; // predictor variables
  real<lower=0> a[p];
  real<lower=0> b[p];
  real<lower=3> d;
  cov_matrix[K] Q;
  real<lower=0> api;
  real<lower=0> bpi;
}
parameters{
  // vector[p] beta[K];
  matrix[K,p] beta;
  // vector<lower=0>[p] tau;
  real<lower=0> tau[p];
  real<lower=0> lambda[p];
  cov_matrix[K] Sigma;
  real<lower=0,upper=1> pi;
}
transformed parameters{
  vector[N] mu[K];
  for(k in 1:K){
    mu[k] = X*beta[k,]';
    mu[k] = Phi(mu[k]);
  }
}
model{
  for(k in 1:K){
    target += bernoulli_lpmf(y[,k] | mu[k]);
  }
  for(j in 1:p){
    target += gamma_lpdf(tau[j] | (K+1)/2 , lambda[j]/2);
    target += gamma_lpdf(lambda[j] | a[j],b[j]);
  }
  for(j in 1:p){
    target += log_mix(1-pi , multi_normal_lpdf(beta[,j] | rep_vector(0,K) , tau[j]*Sigma ) , 0);
  }
  target += beta_lpdf(pi | api,bpi);
  target += inv_wishart_lpdf(Sigma | d , Q);
  
  
}

particularly, I achieve mixture model part by

  for(j in 1:p){
    target += log_mix(1-pi , multi_normal_lpdf(beta[,j] | rep_vector(0,K) , tau[j]*Sigma ) , 0);
  }

where log\_mix(\theta,\lambda_{1},\lambda_{2}) = log(\theta*exp(\lambda_{1}) + (1-\theta)*exp(\lambda_{2})).

Now my problem is, since \beta is a parameter generate directly by stan, I can’t restrict it equally to zero.

I’m new in stan and I can’t get through it, can anybody help me ?

Thank you!!

2

Yeah, you won’t be able to sample a parameter with a distribution like that.

You’ll need to integrate it all the way out. So write your data conditional on \beta and add em’ up:

p(y | \theta) = (1 - \pi) p(y | \theta, \beta \neq 0) + \pi p(y | \theta, \beta = 0)

Evaluating p(y | \theta, \beta = 0) is easy enough – just set \beta = 0 in your eqs and have fun. For the other one, sample another variable around that has the right distribution (\alpha \sim N(0, 1) or whatever) and just use that for the value of \beta in evaluating the p(y | \theta, \beta \neq 0) term.

I hope that’s right at least haha.

3

Thanks for your answer,

But \beta can’t be integrate out. I need to sample \beta from that distribution.

I try to have a indicator \gamma \sim bernoulli(1 - \pi) and if \gamma = 1, I sample \beta from N(0,1), or I set \beta = 0 if \gamma = 0.

But it’s hard to achieve in stan, because sampling such a \gamma I have to write user define functions which ending in _rng. But this function can be only used in transformed data block, generated quantities block.

Is it possible that I can do that in stan?

4

Ack, that’s me being bad with terminology and equations simultaneously. Sorry.

You can’t sample a discrete indicator variable like that. You can integrate out such a variable :D.

That is what I was saying in the last post, I just screwed it up.

So you have an indicator variable, \gamma, you want to integrate out.

p(y, \beta) = p(y | \beta) p(\beta) = \sum_i p(y | \beta) p(\beta | \gamma_i) p(\gamma_i)

So that’s just:

p(y, \beta) = p(y | \beta) p(\beta | \gamma = 0) p(\gamma = 0) + p(y | \beta) p(\beta | \gamma = 0) p(\gamma = 1)

And you’ve specified that p(\beta = 0 | \gamma = 0) = 1 and p(\beta | \gamma = 1) \sim N(0, 1).

Sorry for the confusion! Is that clearer?

edit: Updated eq. to get fancy latex right

6

Sorry, maybe I didn’t declare my question clearly and make you confused.

You mention

That is what I have done in my original question(in my code \beta is a vector). That is ,

    target += log_mix(1-pi , multi_normal_lpdf(beta| rep_vector(0,K) , Sigma ) , 0);

But this can’t restrict \beta = 0 if \gamma = 0 in stan. The code I write can be used to any mixture model in this form \beta \sim (1-\pi)N(0,1) + \pi \delta_{c}(\beta), c can be any real number. But in my model, I have to restrict it as 0 instead of any real number.

In the begging, I did try to generate a \gamma such that I can do it in following way,

  for(j in 1:p){ 
    if(gamma == 1){
      beta[,j] ~ multi_normal(rep_vector(0,K) , tau[j]*Sigma);
    }else{
      beta[,j]  = rep_vector(0,K);
    }
  }

But as what you have say, I have to integrate out \gamma, which is my original code.

I just want to know how to restrict \beta = 0

7

Oooh sorry, try:

parameters {
  real beta1;
}

model {
  real beta;
  beta1 ~ normal(0, 1);
  { // evaluate gamma == 1 part of lpdf
    beta = beta1;
    ...
  }

  { // evaluate gamma == 0 part of lpdf
    beta = 0;
    ...
  }
}

That make sense? I just made beta a scalar for simplicity.