Computing log likelihood in generated quantities block

[jacquie...@gmail.com]

Hi all,

I'm trying to compute log likelihood for my model (reinforcement learning model) to compute the WAIC for model comparison. Everything is working, but I'm not sure whether there would be a more economical way to code it as at the moment I seem to have to repeat most of my model in the 'generated quantities' block - I'm not sure whether there is a way around this? This is a slight problem because it might make this quite prone to errors if I'm changing something in the model block, but then forget to change it in the generated quantities block without it producing any errors etc.

Many thanks

Jacquie

data{

  int nsub;

  int mtr;           // maximal number of trials per person (=240)

  int Dec[nsub,mtr]; // real choices

  real RMagO[nsub,mtr];

  real Prob[nsub,mtr];

  int Symbs[nsub,mtr];   // 1 or 2 are normal symbols, 3 and 4 should be ignored for learning

}

parameters{

  // Group level

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

  real<lower=0,upper=10> mu_beta;

  real<lower=0,upper=1>  mu_relre;    // relative weight of learnt reward to effort

  // Group level standard deviation

  real<lower=0> sd_alpha;

  real<lower=0> sd_beta;

  real<lower=0> sd_relre;

  // Individual-level parameters (will be drawn from group distributions

  real<lower=0,upper=1>  alpha[nsub]; 

  real<lower=0,upper=10> beta[nsub];

  real<lower=0,upper=1>  relre[nsub]; 

}

model{

  // Variables that will be computed in the model

  real RMagP[nsub,mtr,2];

  real EMagP[nsub,mtr,2];

  real Util[nsub,mtr,2];    // Utility     

  real DecP[nsub,mtr];      // Predicted choice (Utility transformed by softmax)

  real log_lik[nsub,mtr]; 

  // Priors - group

  mu_alpha ~ uniform(0,1);

  mu_beta  ~ uniform(0,10);

  mu_relre ~ uniform(0,1);

  sd_alpha ~ uniform(0,1.5); // value from pvl_d.stan model

  sd_beta  ~ uniform(0,1.5);

  sd_relre ~ uniform(0,1.5);

  // Single sub drawn from normal distr.

  for (s in 1:nsub){

    alpha[s] ~ normal(mu_alpha,sd_alpha);

    beta[s]  ~ normal(mu_beta,sd_beta);

    relre[s] ~ normal(mu_relre,sd_relre);

  }

  // Compute learning

  for (s in 1:nsub){

    for (t in 1:mtr){

      if (t==1 || t==121){  // detect first trial of set

RMagP[s,t+1,:]  <- RMagO[s,t,:];

      }

      else if (Dec[s,t]!=1000 && t!=120 && t!=240) { // actual learning

for (o in 1:2){

  if (Symbs[s,t,o]<3){

    RMagP[s,t+1,o] <- RMagP[s,t,o] + alpha[s]*(RMagO[s,t,o]-RMagP[s,t,o]);

          }

  else {

    RMagP[s,t+1,o] <- RMagP[s,t,o];

  }

}

      }

    }

  }

  // Compute utility and assign decisions (only on normal trials

  for (s in 1:nsub){

    for (t in 1:mtr){

      if (Dec[s,t]!=1000 && Symbs[s,t,1]<3 && Symbs[s,t,2]<3&& t!=1 &&t!=121){

for (o in 1:2){

  Util[s,t,o] <- relre[s]*RMagP[s,t,o] + Prob[s,t,o];

}

DecP[s,t] <- 1.0/(1.0+exp(beta[s]*(Util[s,t,2]-Util[s,t,1])));

        Dec[s,t] ~ bernoulli(DecP[s,t]);

      }

    }

  }

}

generated quantities{

  real RMagP[nsub,mtr,2];

  real Util[nsub,mtr,2];    // Utility     

  real DecP[nsub,mtr];      // Predicted choice (Utility transformed by softmax)

  real log_lik[nsub,mtr];

  // Compute learning

  for (s in 1:nsub){

    for (t in 1:mtr){

      if (t==1 || t==121){  // detect first trial of set

   RMagP[s,t+1,:]  <- RMagO[s,t,:];

      }

      else if (Dec[s,t]!=1000 && t!=120 && t!=240) { // actual learning

   for (o in 1:2){

     if (Symbs[s,t,o]<3){

       RMagP[s,t+1,o] <- RMagP[s,t,o] + alpha[s]*(RMagO[s,t,o]-RMagP[s,t,o]);

          }

     else {

       RMagP[s,t+1,o] <- RMagP[s,t,o];

     }

   }

      }

    }

  }

  // Compute utility and assign decisions (only on normal trials

  for (s in 1:nsub){

    for (t in 1:mtr){

      if (Dec[s,t]!=1000 && Symbs[s,t,1]<3 && Symbs[s,t,2]<3&& t!=1 &&t!=121){

   for (o in 1:2){

     Util[s,t,o] <- relre[s]*RMagP[s,t,o] + Prob[s,t,o];

   }

   DecP[s,t] <- 1.0/(1.0+exp(beta[s]*(Util[s,t,2]-Util[s,t,1])));

log_lik[s,t] <- bernoulli_log(Dec[s,t],DecP[s,t])

      }

    }

  }

}

[Jonah Gabry]

Hi Jacquie, 

I suggest creating a user-defined function in the functions block that computes the log-likelihood. Then you can use it in the model block and generated quantities. The "User Defined Functions" section of manual goes through how to do this, but if you run into any trouble let us know. 

Jonah

[Andrew Gelman]

Yes, but this doesn’t always make sense because it can give up the efficiency that would be gained by vectorization. In my examples I found it made more sense to just use near-duplicate code, ugly as that is, with vectorized computation in the model block and point-by-point computation in the generated quantities block. We discuss this in our paper.
A

[Jonah Sol Gabry]

One way to get around that is to use a function to generate everything you need before computing the log-likelihood (e.g. the linear predictor) and then use it with the vectorized densities in the model block and loop over it in generated quantities.

--
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/auSQF5xMjyc/unsubscribe.
To unsubscribe from this group and all its topics, 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.

[Jonah Sol Gabry]

Or probably even better, write a vectorized version of the whole thing and a non-vectorized version as user defined functions, and use the first in model and second in generated quantities. If you're using RStan you can use expose_stan_functions to access the functions from R and test them. 

[Aki Vehtari]

On Thursday, February 4, 2016 at 12:32:17 AM UTC+2, Jonah Gabry wrote:

One way to get around that is to use a function to generate everything you need before computing the log-likelihood (e.g. the linear predictor) and then use it with the vectorized densities in the model block and loop over it in generated quantities.

How cool it would be, if in the future we could use a vector argument both in the model block and in the generated quantities block and have an option to get output as sum of log densities or vector of log densities..?

Aki

 

On Wednesday, February 3, 2016, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Yes, but this doesn’t always make sense because it can give up the efficiency that would be gained by vectorization.  In my examples I found it made more sense to just use near-duplicate code, ugly as that is, with vectorized computation in the model block and point-by-point computation in the generated quantities block.  We discuss this in our paper.
A

> On Feb 3, 2016, at 11:57 AM, Jonah Gabry <jga...@gmail.com> wrote:
>
> Hi Jacquie,
>
> I suggest creating a user-defined function in the functions block that computes the log-likelihood. Then you can use it in the model block and generated quantities. The "User Defined Functions" section of manual goes through how to do this, but if you run into any trouble let us know.
>
> Jonah


--
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/auSQF5xMjyc/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[jacquie...@gmail.com]

This is a naive question, but does it look like I could make a vectorized version of these for loops (for example in the section with 'compute learning')?  I assumed I couldn't because each entry in row t+1 of the array (e.g. RMagP) uses what is written in row t (and there are other 'if' conditions).

[Jacqueline Scholl]

actually, it turns out there is one more problem with expressing everything in the generated quantities block: as there are so many variables (of which I only want the log_lik), the output.csv files become huge - with 30000 iterations and 15000 warmup, thinning =1, (i've put refresh to 3000, but not sure this matters?), the file got to 21 GB for 1 chain (the other two chains were 6GB and 13GB). Is there any way to avoid this?

--
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/auSQF5xMjyc/unsubscribe.

To unsubscribe from this group and all its topics, 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.

--

Jacqueline Scholl

D.Phil Candidate, Department of Psychiatry, University of Oxford, OX3 7JX

[Jacqueline Scholl]

actually, using the function block as you advised solves this - only the 'returned' variable from the function is stored. thank you again for the tip!

[Jonah Sol Gabry]

On Thursday, February 4, 2016, Aki Vehtari <aki.v...@aalto.fi> wrote:

How cool it would be, if in the future we could use a vector argument both in the model block and in the generated quantities block and have an option to get output as sum of log densities or vector of log densities..?

Very cool indeed.

[Jonah Sol Gabry]

On Thursday, February 4, 2016, Jacqueline Scholl <jacquie...@gmail.com> wrote:

actually, using the function block as you advised solves this - only the 'returned' variable from the function is stored. thank you again for the tip!

Glad that's working for you. In general, you can also create local blocks using { } inside of generated quantities if you only want to store a subset of quantities you define. For example 

generated quantities {

  real x;

  x <- 1.5;

 {

  real y;

  y <- 0.5;

  x <- x + y;

 }

} 

will only save x (with a value of 2) and y will be dropped since it's defined in a local block. 

Hope that helps,

Jonah

[Bob Carpenter]

> On Feb 4, 2016, at 2:11 AM, Aki Vehtari <aki.v...@aalto.fi> wrote:
>
> On Thursday, February 4, 2016 at 12:32:17 AM UTC+2, Jonah Gabry wrote:
> One way to get around that is to use a function to generate everything you need before computing the log-likelihood (e.g. the linear predictor) and then use it with the vectorized densities in the model block and loop over it in generated quantities.
>
> How cool it would be, if in the future we could use a vector argument both in the model block and in the generated quantities block and have an option to get output as sum of log densities or vector of log densities..?

Until we get more C++ programmers willing to work on core
Stan infrastructure, there's no way this is going to happen.

- Bob

[Bob Carpenter]

:-) Should've read Jonah's response first!

- Bob

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

[Bob Carpenter]

You can also add local variables in blocks, as in:

generated quantities {
vector[N] theta; // printed
{
matrix[N, N] Sigma; // not printed
...

}
}

- Bob

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

[Jonah Sol Gabry]

On Thursday, February 4, 2016, Bob Carpenter <ca...@alias-i.com> wrote:

:-)  Should've read Jonah's response first!

- Bob

Just trying to chip in and lighten your stan-users question load. But you know, if you're just going to answer everything regardless ...