problem with mixture GLM

[Chris]

Hi,

I'm trying to build a GLM as a mixture of two Gaussians in R and STAN.

here's the code:

# generate data

set.seed(1)

N <- 1000

df <- data.frame(x1=runif(N,-1,1), x2=runif(N,-1,1), x3=runif(N,-1,1))

X <- model.matrix(~.,df)          # the model matrix (contains intercept)

betas1 <- runif(4, -1,1)  

betas2 <- runif(4, 2,4)

r.norm.mix <- function(n,p,m1,m2,sigma) {    

  weight <- runif(n)<p

  rnorm(n,mean=ifelse(weight,m1,m2),sd=sigma)

}

y <- r.norm.mix(1000, p=0.8, m1=X%*%betas1, m2=X%*%betas2, sigma=1)

hist(y,30) # visualization

# STAN model

library(rstan)

model_string <-"

data {

  int<lower=1> K; // number of mixture components

  int<lower=1> N; // number of data points

  int<lower=1> M; //the number of columns in the model matrix

  matrix[N,M] X; //the model matrix

  real y[N]; // observations

}

parameters {

  simplex[K] theta; // mixing proportions

  vector[M] betas1;

  vector[M] betas2;

  real<lower=0> sigma[K]; // sigmas of mixture components

}

model {

  for(m in 1:M) {

    betas1[m] ~ normal(0,100);

    betas2[m] ~ normal(0,100);

  }

  sigma ~ cauchy(0,2.5);

  for (n in 1:N) {

    increment_log_prob(log_mix(theta,

      normal_log(y[n], X[n,]*betas1, sigma[1]),

      normal_log(y[n], X[n,]*betas2, sigma[2])));

  }

}

"

m <- stan(model_code = model_string, data = list(X=X, M=ncol(X), K=2, y = y, N = nrow(X)))

this model immediately crashes R and I'm out of ideas. Thanks

[David Westergaard]

Hi Chris,

If you read the manual, log_mix is defined for binary mixtures, and hence only takes a single real as input. I'm not sure if that's what is causing the crashes, but it might be related!

Cheers,

David

[David Muller]

Hi Chris,

I'm surprised that this is crashing R, but the first argument of log_mix needs to be a real, and you are passing it a vector (simplex).

Cheers,

David

[Chris]

thanks to both of you,

changing the simplex to

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

fixed it.

damn label switching though :)

[Daniel Lee]

Chris, can you tell us a little more about crashing R?

What version of R are you using?

What operating system?

What R GUI? (RStudio, command line, something else?)

When R crashes, what does it do? Just quits? Mentions a seg fault? Something else?

Thanks.

Daniel

--
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.

[Chris]

I'm using Microsoft R Open 3.2.3 on Ubuntu 15.04 with Rstudio (newest version)

I get the "bomb" symbol in Rstudio telling me to restart R immediately after executing the wrong Stan code.

Nothing else.

[Bob Carpenter]

Do you have a script that recreates the problem?

I'm not sure what's going on because Stan correctly rejects
your code with the correct error message:

--- Translating Stan model to C++ code ---
bin/stanc /Users/carp/temp2/chris.stan --o=/Users/carp/temp2/chris.hpp
Model name=chris_model
Input file=/Users/carp/temp2/chris.stan
Output file=/Users/carp/temp2/chris.hpp

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for:

log_mix(vector, real, real)

Available argument signatures for log_mix:

log_mix(real, real, real)


ERROR at line 24

22: increment_log_prob(log_mix(theta,
23: normal_log(y[n], X[n,]*betas1, sigma[1]),
24: normal_log(y[n], X[n,]*betas2, sigma[2])));
^
25: }

- Bob

P.S. What is Microsoft R Open and why would you use that other
than the standard R?

[Chris]

the script is basically what I posted above.

I did not compile it on the command line and I always use the model_code variant to pass the model string code.

But thanks for trying out my code, maybe that is a future approach for debugging crashes.

Microsoft R Open is a 100% compatible R version which uses the intel math kernel libraries for additional speed and parallelization.

Never had any problems with it.

[Bob Carpenter]

OK, I ran everything you gave me in the latest release
version of RStan on my Mac, and it worked the same as on the
command line:

> m <- stan(model_code = model_string, data = list(X=X, M=ncol(X), K=2, y = y, N = nrow(X)))

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

No matches for:

log_mix(vector, real, real)

Available argument signatures for log_mix:

log_mix(real, real, real)

ERROR at line 25

23: increment_log_prob(log_mix(theta,
24: normal_log(y[n], X[n,]*betas1, sigma[1]),
25: normal_log(y[n], X[n,]*betas2, sigma[2])));
^
26: }

Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model '832efd089eaa99364234f245b82fce72' due to the above error.
>


So it's not a language issue or a basic RStan issue.

Did you specialize the R and RStan config to use the Intel compiler?

- Bob

[Chris]

here is a screenshot:

http://i.imgur.com/tR12lS2.png

rstan is compiled with gcc, only MS R uses intel math kernel libraries.

Bob, can you help me extend my mixture model to a case where the number of mixture

components K is not fixed, but estimated as well.

this cannot be done with log_mix as it only accepts two components, correct?

How would this be done then?

[Daniel Lee]

Hi Chris,

On Feb 12, 2016, at 3:35 AM, Chris <wel...@gmail.com> wrote:

here is a screenshot:

http://i.imgur.com/tR12lS2.png

This looks like a weird interaction between Microsoft R Open and RStudio. If you wanted to help hunt it down, try running the code in Microsoft R Open without RStudio and see if the error message is properly displayed. 

If it doesn't, send a bug report there. If it does, send a bug report over to RStudio (try the release version before filing bugs to daily builds, though). 

rstan is compiled with gcc, only MS R uses intel math kernel libraries.

Bob, can you help me extend my mixture model to a case where the number of mixture

components K is not fixed, but estimated as well.

this cannot be done with log_mix as it only accepts two components, correct?

How would this be done then?

There's a section in the manual called Latent Discrete Parameters. It has explicit math and instructions on how to marginalize over a fixed K. 

[Chris]

just wanted to add that the crashing must be ubuntu+rstudio specific.

In the R command line I get no crashes, only the syntax errors.

On a windows machine I get no crashes in rstudio and the syntax errors.

[David Muller]

When I run your original code in GNU R from Rstudio on Ubuntu I get the expected error message, and no crash, so there's definitely something weird going on with MS R Open.

[Bob Carpenter]

Thanks for trying that out. I don't use Rstudio or Ubuntu,
so won't be much help there.

To answer the question about extending beyond two-component
mixtures, see the manual chapter on mixture modeling. Basic
idea is that lambda is a K-simplex and if lp is a K-vector of
densities for the components, then you want

increment_log_prob(log_sum_exp(log(lambda) + lp));

- Bob

[Chris]

Hi Bob and others,

I got a three component mixture working with

  for (n in 1:N) {

    increment_log_prob(log_sum_exp(log(theta) +

      normal_log(y[n], X[n,]*betas1, sigma[1]) +

      normal_log(y[n], X[n,]*betas2, sigma[2]) +

      normal_log(y[n], X[n,]*betas3, sigma[3])));

  }

Can you explain to me please how I would proceed to define a mixture where the number of mixture components K

is also estimated?

And how do I calculate the log likelihood for this model in the generated quantities part?

the log likelihood of a three component mixture with distributions f and mixing proportion p for n observations is log(prod_n(p1*f1 + p2*f2 + p3*f3)).

thanks

Chris

[Bob Carpenter]

That's

> On Mar 1, 2016, at 7:40 AM, Chris <wel...@gmail.com> wrote:
>
> Hi Bob and others,
>
> I got a three component mixture working with
>
> for (n in 1:N) {
> increment_log_prob(log_sum_exp(log(theta) +
> normal_log(y[n], X[n,]*betas1, sigma[1]) +
> normal_log(y[n], X[n,]*betas2, sigma[2]) +
> normal_log(y[n], X[n,]*betas3, sigma[3])));
> }
>

I'm afraid that's not a 3-component mixture. Assuming theta is
a 3-simplex, the mixture would look as follows (assuming betas is made
into a vector):

for (n in 1:N) {

vector[3] lp;
lp <- log(theta);
for (k in 1:3)
lp[k] <- normal_log(y[n], X[n] * betas[k], sigma[k]);
increment_log_prob(log_sum_exp(lp));
}

There's a discussion of why the algebra is this way
in the manual chapter on mixtures.

> Can you explain to me please how I would proceed to define a mixture where the number of mixture components K
> is also estimated?

If you're thinking of Dirichlet processes, we don't really do
that in Stan. You could technically make a mixture of mixtures
and evaluate this, but it's usually not worth the effort.

> And how do I calculate the log likelihood for this model in the generated quantities part?

Just follow the likelihood definition in the model block,
only use foo_log(y, ...) rather than y ~ foo(...).

> the log likelihood of a three component mixture with distributions f and mixing proportion p for n observations is log(prod_n(p1*f1 + p2*f2 + p3*f3)).

I don't know what prod_n is.

- Bob

[Chris]

Am Mittwoch, 2. März 2016 00:14:58 UTC+1 schrieb Bob Carpenter:

That's
> On Mar 1, 2016, at 7:40 AM, Chris <wel...@gmail.com> wrote:
>
> Hi Bob and others,
>
> I got a three component mixture working with
>
>   for (n in 1:N) {
>     increment_log_prob(log_sum_exp(log(theta) +
>       normal_log(y[n], X[n,]*betas1, sigma[1]) +
>       normal_log(y[n], X[n,]*betas2, sigma[2]) +
>       normal_log(y[n], X[n,]*betas3, sigma[3])));
>   }
>

I'm afraid that's not a 3-component mixture.  Assuming theta is
a 3-simplex, the mixture would look as follows (assuming betas is made
into a vector):

  for (n in 1:N) {
    vector[3] lp;
    lp <- log(theta);
    for (k in 1:3)
      lp[k] <- normal_log(y[n], X[n] * betas[k], sigma[k]);
    increment_log_prob(log_sum_exp(lp));
  }

There's a discussion of why the algebra is this way
in the manual chapter on mixtures.

thanks, I understand that now.

However, in your model definition the same betas are used for every component.

Assuming I have K mixture components and M betas each, I did

  vector[M] betas[K];

It does not work with a matrix[K,M] betas; 

 

> Can you explain to me please how I would proceed to define a mixture where the number of mixture components K
> is also estimated?

If you're thinking of Dirichlet processes, we don't really do
that in Stan.  You could technically make a mixture of mixtures
and evaluate this, but it's usually not worth the effort.

> And how do I calculate the log likelihood for this model in the generated quantities part?

Just follow the likelihood definition in the model block,
only use foo_log(y, ...) rather than y ~ foo(...).

so in the above model definition where do I need another log, I already have normal_log. 

Is it log(log_sum_exp(lp)) ?

 

thanks

[Bob Carpenter]

> On Mar 2, 2016, at 4:10 AM, Chris <wel...@gmail.com> wrote:
>
>
>
> Am Mittwoch, 2. März 2016 00:14:58 UTC+1 schrieb Bob Carpenter:
> That's
> > On Mar 1, 2016, at 7:40 AM, Chris <wel...@gmail.com> wrote:
> >
> > Hi Bob and others,
> >
> > I got a three component mixture working with
> >
> > for (n in 1:N) {
> > increment_log_prob(log_sum_exp(log(theta) +
> > normal_log(y[n], X[n,]*betas1, sigma[1]) +
> > normal_log(y[n], X[n,]*betas2, sigma[2]) +
> > normal_log(y[n], X[n,]*betas3, sigma[3])));
> > }
> >
>
> I'm afraid that's not a 3-component mixture. Assuming theta is
> a 3-simplex, the mixture would look as follows (assuming betas is made
> into a vector):
>
> for (n in 1:N) {
> vector[3] lp;
> lp <- log(theta);
> for (k in 1:3)
> lp[k] <- normal_log(y[n], X[n] * betas[k], sigma[k]);
> increment_log_prob(log_sum_exp(lp));
> }
>
> There's a discussion of why the algebra is this way
> in the manual chapter on mixtures.
>
> thanks, I understand that now.
> However, in your model definition the same betas are used for every component.

My mistake --- I meant "assuming betas is made into a matrix", but what
I should've said is array of vectors so we don't need an extra transpose
and it will work more efficiently by reference (see the manual chapter on
arrays vs. vectors):

vector[M] betas[3];

Then X[n] is a row vector and beta[k] is a column vector, so
the product is a scalar.

- Bob

[Chris]

My mistake --- I meant "assuming betas is made into a matrix", but what
I should've said is array of vectors so we don't need an extra transpose
and it will work more efficiently by reference (see the manual chapter on
arrays vs. vectors):

  vector[M] betas[3];

Then X[n] is a row vector and beta[k] is a column vector, so
the product is a scalar.

- Bob

thanks for making that more clear. Could you please advise how I turn my likelihood into a log-likelihood in the generated quantities section.

I'm a little confused because I have so many "logs" in my likelihood already.

This is the model likelihood:

[Jonah Gabry]

On Wednesday, March 2, 2016 at 2:52:31 PM UTC-5, Chris wrote:

thanks for making that more clear. Could you please advise how I turn my likelihood into a log-likelihood in the generated quantities section.

I'm a little confused because I have so many "logs" in my likelihood already.

This is the model likelihood:

>   for (n in 1:N) { 
>     vector[3] lp; 
>     lp <- log(theta); 
>     for (k in 1:3) 
>       lp[k] <- normal_log(y[n], X[n] * betas[k], sigma[k]); 
>     increment_log_prob(log_sum_exp(lp)); 
>   } 

 

I might be missing something but I think you want to add to lp[k] instead of redefining it. Right now you set lp equal to log(theta) but then set it equal to something else when you set lp[k] to the result of normal_log(...). So I think you want lp[k] <- lp[k] + normal_log(...). 

In generated quantities you should be able to do exactly the same thing you do in the model block but instead of incrementing the log probability you just want to save its value.

Jonah

[Bob Carpenter]

Thanks, Jonah --- that should definitely be

lp[k] <- lp[k] + normal_log(...);

I probably messed up the original suggestion.

You don't want to repeat the code, though. Break
the likelihood out into a function and just call it
in both the model block and generated quantities block.

Or you could define the likelihood as a transformed
parameter and then just use increment_log_prob(...) on
it in the model block.

- Bob

[Chris]

thanks, I corrected that.

for a start I just copy pasted the model part to the generated quantities part to get the loglikelihood;

generated quantities {

  vector[N] log_lik;

  for (n in 1:N) {

    vector[K] lp;

    for (k in 1:K) {

      lp[k] <- lp[k] + normal_log(y[n], X[n] * betas[k], sigma[k]);

    }

    log_lik[n] <- log_sum_exp(lp);

  }

}

the model part evaluated without problems but finally fails with these errors:

...

Chain 1, Iteration: 1800 / 2000 [ 90%]  (Sampling)

Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)# 

#  Elapsed Time: 370.208 seconds (Warm-up)

#                203.205 seconds (Sampling)

#                573.413 seconds (Total)

# 

Error in validObject(.Object) : 

  invalid class “stanfit” object: The following variables have undefined values:  log_lik[1],The following variables have undefined values:  log_lik[2],The following variables have undefined values:  log_lik[3],The following variables have undefined values:  log_lik[4],The following variables have undefined values:  log_lik[5],The following variables have undefined values:  log_lik[6],The following variables have undefined values:  log_lik[7],The following variables have undefined values:  log_lik[8],The following variables have undefined values:  log_lik[9],The following variables have undefined values:  log_lik[10],The

I tried "log_lik[n] <- sum(lp);" but that fails as well.

[Bob Carpenter]

Local variables need to be defined before they're used,
so your problem is this line

lp[k] <- lp[k] + ...

which is executed before lp[] gets a value. You're probably
missing the mixing proportions:

log_lik <- log(theta);

or if they're uniform, you want to just set, not increment.

- Bob

[Chris]

thanks, I got it now