Problems with longer iterations in Mixture Model

[smockin morry]

Hi Stan Community.

The attached model and data seems to work fine on 2 chains only with 500 iterations but poor convergence. It throws an error whenever i attempt to extend iterations say to 1000 or even using more chains say 3.  

Where could be the problem, model specification?  Kindly help.

Thanks in advance.

 # Data
stanData= readRDS("mixtureModelData.rds")

summary(stanData$y)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.01953 0.01953 0.02237 0.05813 0.08937 0.22110

 #This works

mix1= stan(file = "getMixtures.stan", data = stanData, chains = 2, iter = 500)

 #This throws an error
mix2= stan(file = "getMixtures.stan", data = stanData, chains = 3, iter = 500)

 #This throws an error
mix3= stan(file = "getMixtures.stan", data = stanData, chains = 2, iter = 1000)


#Error

Error in checkForRemoteErrors(val) : 
  one node produced an error: invalid class “stanfit” object: The following variables have undefined values: 
 p[1,1],The following variables have undefined values: 
 p[2,1],The following variables have undefined values:
  p[3,1],The following variables have undefined values: 
 p[4,1],The following variables have undefined values: 
 p[5,1],The following variables have undefined values: 
 p[6,1],The following variables have undefined values: 
 p[7,1],The following variables have undefined values: 
 p[8,1],The following variables have undefined values: 
 p[9,1],The following variables have undefined values: 
 p[10,1],The following variables have undefined values: 
 p[11,1],The following variables have undefined values:
  p[12,1],The following variables have undefined values: 
 p[13,1],The following variables have undefined values: 
 p[14,1],The following variables have undefined values: 
 p[15,1],The following variables have undefined values: 
 p[16,1],The following variables have undefined 

Session info

R version 3.3.2 (2016-10-31)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS release 6.8 (Final)

locale:
[1] C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rstan_2.14.1         StanHeaders_2.14.0-1 ggplot2_2.2.1        magrittr_1.5        

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.9      codetools_0.2-15 digest_0.6.12    assertthat_0.1   grid_3.3.2       plyr_1.8.4       gtable_0.2.0     stats4_3.3.2    
 [9] scales_0.4.1     lazyeval_0.2.0   labeling_0.3     tools_3.3.2      munsell_0.4.3    parallel_3.3.2   inline_0.3.14    colorspace_1.3-2
[17] gridExtra_2.2.1  tibble_1.2 

[Bob Carpenter]

If you put RStan in the subject, you'll get the attention of
the RStan devs. We're getting buried in messages and they're
hard to triage.

Can you run one chain for 1000 iterations? If so, then you're
either having problem with memory or some interaction with RStan's
internals that you'll need an RStan dev to help disentangle.

- Bob

> R version 3.3.2 (2016-10-31)
> Platform: x86_64-redhat-linux-gnu (64-bit)
> Running under: CentOS release 6.8 (Final)
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] rstan_2.14.1 StanHeaders_2.14.0-1 ggplot2_2.2.1 magrittr_1.5
>
> loaded via a namespace (and not attached):
> [1] Rcpp_0.12.9 codetools_0.2-15 digest_0.6.12 assertthat_0.1 grid_3.3.2 plyr_1.8.4 gtable_0.2.0 stats4_3.3.2
> [9] scales_0.4.1 lazyeval_0.2.0 labeling_0.3 tools_3.3.2 munsell_0.4.3 parallel_3.3.2 inline_0.3.14 colorspace_1.3-2
> [17] gridExtra_2.2.1 tibble_1.2
>
>

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
> <getMixtures.stan><mixtureModelData.rds>

[smockin morry]

Hi Bob,

One chain with 1000 iterations works fine. Am using 64 cores machine with 214G of RAM. Probably Memory cannot be an issue here.

[Ben Goodrich]

On Monday, March 20, 2017 at 9:57:16 AM UTC-4, smockin morry wrote:

Where could be the problem, model specification?  Kindly help.

My guess is that this in generated quantities:

generated quantities {
  matrix[N, k_groups] p;
  for (n in 1:N){
    vector[k_groups] ps;
    for (K in 1:k_groups){
      ps[K] = theta[K]*exp(normal_lpdf(y[n] | mu[K], sigma));
    }
    for (k in 1:k_groups){
      p[n, k]=ps[k]/sum(ps);
    }
  }
}

the last line overflows to infinity / infinity. You would be better off doing

ps[K] = log(theta[K]) + normal_lpdf(y[n] | mu[K], sigma);

and then

p[n, ] = transpose(softmax(ps));

Ben

[smockin morry]

Thanks Ben. That fixed the problem.

[smockin morry]

Hi everyone

Attached traceplot seem weird. 2 chains mix nicely but they seem parallel to the 3rd chain? why such a behavior even after running 20000 iterations?

Thanks in advance.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/KEowxLcVK1U/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[stijn masschelein]

That looks like an identification problem which is to be expected in mixture models. Have a look at Michael's excellent case study for potential solutions. http://mc-stan.org/documentation/case-studies/identifying_mixture_models.html

[smockin morry]

Thanks Masschelein for the response.

 I have already implemented one of the resolutions, that is by enforcing the ordering of means of components in the model as  suggested 

by Michael in the link you provided. Just to understand where the problem could be, i have changed Michael's  binary mixture model to a general mixture model as 

provided in the Stan manual  and both approaches provides similar estimates of parameters as shown below. 

But the problem of label identification seems apparent to a general mixture model when

 the number of components K is greater than 2 as shown below. 

1. Does it mean that solutions to label switching only works in a binary mixture model?  

2. How can the general mixture model be tweaked to avoid label switching menace when is say K>3

Binary mixture model approach

binary_mix_ordered<-"

data {

 int<lower = 0> N;

 vector[N] y;

}

parameters {

  ordered[2] mu;

  real<lower=0> sigma[2];

  real<lower=0, upper=1> theta;

}

model {

 sigma ~ normal(0, 2);

 mu ~ normal(0, 2);

 theta ~ beta(5, 5);

 for (n in 1:N)

   target += log_mix(theta,

                     normal_lpdf(y[n] | mu[1], sigma[1]),

                     normal_lpdf(y[n] | mu[2], sigma[2]));

}"

binary_mix_fit <- stan(model_code =binary_mix_ordered

                    ,chains=4

                    , seed=48389292

                    , refresh=2000)

traceplot(binary_mix_fit)

print(binary_mix_fit)

             mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
mu[1]       -2.73    0.00 0.04    -2.82    -2.76    -2.73    -2.71    -2.65  3355    1
mu[2]        2.87    0.00 0.05     2.77     2.83     2.87     2.90     2.97  4000    1
sigma[1]     1.03    0.00 0.03     0.97     1.01     1.03     1.05     1.09  4000    1
sigma[2]     1.02    0.00 0.04     0.95     1.00     1.02     1.05     1.10  4000    1
theta        0.62    0.00 0.02     0.59     0.61     0.62     0.63     0.65  4000    1
lp__     -2106.86    0.04 1.58 -2110.69 -2107.67 -2106.56 -2105.73 -2104.78  1896    1

General mixture model approach on binary mixture (K=2)

#General mixture according to the manual

general_mix_ordered="

data {

  int<lower = 0> N;

  vector[N] y;

  int<lower=1> K;

}

parameters {

  ordered[K] mu;

  real<lower=0> sigma[K];

  //real<lower=0, upper=1> theta[K];

  simplex[K] theta;

}

model {

  real contr[K];

  sigma ~ normal(0, 3);

  mu ~ normal(0, 3);

  //theta ~ beta(5, 5);

  theta ~ dirichlet(rep_vector(5, K));

  for (n in 1:N){

    for(k in 1:K){

      contr[k]=log(theta[k])+normal_lpdf(y[n] | mu[k], sigma[k]);

    }

    target += log_sum_exp(contr);

  }

}"

K=2

general_2.mix_ordered_Fit <- stan(model_code =general_mix_ordered

                      ,chains=4

                      , seed=48389292

                      , refresh=2000)

traceplot(general_2.mix_ordered_Fit)

print(general_2.mix_ordered_Fit)

             mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff Rhat
mu[1]       -2.73    0.00 0.04    -2.82    -2.76    -2.73    -2.70    -2.65  3501    1
mu[2]        2.87    0.00 0.05     2.77     2.83     2.87     2.91     2.98  4000    1
sigma[1]     1.03    0.00 0.03     0.97     1.01     1.03     1.05     1.10  4000    1
sigma[2]     1.02    0.00 0.04     0.95     1.00     1.02     1.05     1.11  4000    1
theta[1]     0.62    0.00 0.02     0.59     0.61     0.62     0.63     0.65  4000    1
theta[2]     0.38    0.00 0.02     0.35     0.37     0.38     0.39     0.41  4000    1
lp__     -2105.70    0.04 1.62 -2109.80 -2106.54 -2105.39 -2104.49 -2103.56  2113    1

General mixture model approach using 3 Components: Traceplot exhibit label switching problem 

K=3

general_3.mix_ordered_Fit <- stan(model_code =general_mix_ordered

                    ,chains=4

                    , seed=48389292

                    , refresh=2000)

traceplot(general_3.mix_ordered_Fit)

print(general_3.mix_ordered_Fit)

             mean se_mean   sd     2.5%      25%      50%      75%    97.5% n_eff  Rhat
mu[1]       -2.91    0.08 0.18    -3.29    -3.04    -2.86    -2.76    -2.68     5  1.30
mu[2]       -1.15    1.58 2.24    -2.74    -2.59    -2.39    -0.63     2.86     2 11.90
mu[3]        2.91    0.06 0.12     2.78     2.84     2.89     2.94     3.26     4  1.31
sigma[1]     1.02    0.01 0.09     0.85     0.97     1.02     1.06     1.21    96  1.04
sigma[2]     0.95    0.02 0.16     0.66     0.86     0.95     1.03     1.30   106  1.04
sigma[3]     1.02    0.00 0.09     0.84     0.99     1.02     1.06     1.22  1030  1.00
theta[1]     0.40    0.09 0.15     0.15     0.29     0.38     0.58     0.64     3  1.83
theta[2]     0.27    0.03 0.10     0.10     0.19     0.26     0.34     0.47     9  1.12
theta[3]     0.33    0.06 0.09     0.11     0.33     0.37     0.38     0.41     2  3.02
lp__     -2118.00    0.98 2.55 -2123.61 -2119.56 -2117.79 -2116.08 -2113.77     7  1.20

Can someone help here ,

Thanks.

[Bob Carpenter]

No, the ordering should work for any number of components, but
as you can see in Michael's case study, there are lots of identifiability
issues lurking if the means being ordered are close to one another.

- Bob

> rstan_options(auto_write=T)
> options(mc.cores=parallel::detectCores()*2)

> rstan_options(auto_write=T)
> options(mc.cores=parallel::detectCores()*2)

>
> general_3.mix_ordered_Fit <- stan(model_code =general_mix_ordered
> ,chains=4
> , seed=48389292
> , refresh=2000)
>
> traceplot(general_3.mix_ordered_Fit)
> print(general_3.mix_ordered_Fit)
> mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
> mu[1] -2.91 0.08 0.18 -3.29 -3.04 -2.86 -2.76 -2.68 5 1.30
> mu[2] -1.15 1.58 2.24 -2.74 -2.59 -2.39 -0.63 2.86 2 11.90
> mu[3] 2.91 0.06 0.12 2.78 2.84 2.89 2.94 3.26 4 1.31
> sigma[1] 1.02 0.01 0.09 0.85 0.97 1.02 1.06 1.21 96 1.04
> sigma[2] 0.95 0.02 0.16 0.66 0.86 0.95 1.03 1.30 106 1.04
> sigma[3] 1.02 0.00 0.09 0.84 0.99 1.02 1.06 1.22 1030 1.00
> theta[1] 0.40 0.09 0.15 0.15 0.29 0.38 0.58 0.64 3 1.83
> theta[2] 0.27 0.03 0.10 0.10 0.19 0.26 0.34 0.47 9 1.12
> theta[3] 0.33 0.06 0.09 0.11 0.33 0.37 0.38 0.41 2 3.02
> lp__ -2118.00 0.98 2.55 -2123.61 -2119.56 -2117.79 -2116.08 -2113.77 7 1.20
>
>

> Thanks,

>
> On Thursday, March 30, 2017 at 5:02:57 AM UTC+3, stijn masschelein wrote:
> That looks like an identification problem which is to be expected in mixture models. Have a look at Michael's excellent case study for potential solutions. http://mc-stan.org/documentation/case-studies/identifying_mixture_models.html
>

> --

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.