a Waic example!

[Andrew Gelman]

Hi all. I've programmed a simple example of WAIC for the example on p.165 and p.177 of BDA3. Attached are three files:
waic_example.R
lm_waic.stan
hibbs.dat

I have some thoughts on how we could alter Stan to do this automatically. The simplest approach would be to just put a function waic() in rstan that looks in the output for the variable "log_lik" and, if such a variable exists, calculates Waic. A more comprehensive approach would allow tagging of the Stan model so that certain lines of code would be identified with the loglikelihood so that "log_lik" would get calculated automatically.

By the way: in calculating Waic, we do want to use the full loglik (including the normalizing factor); this gets done automatically when we use increment_log_prob in Stan. If we did it using tagging, we'd just have to be sure to remember to include those factors.

---

[Andrew Gelman]

We could do Waic in rstan as is (just using the function I already wrote!), conditional on the user creating the log_lik vector correctly in the Stan program. But maybe it makes sense to get it working in regular Stan using tagging. That sort of user intervention (tagging lines in the model as components in the loglik vector) would be necessary in any case. There's no way that it could be done automatically given only the log posterior function with no labeling.

I do think, however, that we should be able to do it without setting up intermediate log probability variables. I think it should be possible (at least for many many common examples) to just have the user tag the appropriate lines in the model as being log-likleihood terms.

A

On Dec 28, 2013, at 10:59 PM, Bob Carpenter wrote:

>
>
> On 12/28/13, 4:49 PM, Shige Song wrote:
>> Will be great if it can be incorporated into the rstan package.
>
> It'll have to be changed in basic Stan first so that
> there's a way to specify the log probabilities in the algorithm.
>
> As Andrew brought up, it's a bit tricky now because of the interpretation
> of up-to-a-proportion calcs. One thing I could do is template out the
> propto template parameter for both the sampling statements and
> log probability functions and let that be specified at compile time.
> But there's still the issue of how to divide the data up and specify the
> model. It won't be something we can do without user intervention in
> setting up all the intermediate log probability variables needed.
>
> - 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.
> For more options, visit https://groups.google.com/groups/opt_out.

[Andrew Gelman]

Hi--I don't know why waic is inf, that sounds like a mistake or a problem, possibly something where the expected value is 0 but the data are nonzero?  Also I think there's a mistake in your model in that it looks like you're counting the data twice.  If you include the increment_log_prob thing, you don't need to also express your model using ~.

On Jan 6, 2014, at 8:40 PM, Linas Mockus wrote:

Thanks for posting the example. 

I am trying to compare two models - gamma and poisson for count data.  waic for poisson model is Inf while waic for gamma model is 2424.  I assume that it means that gamma model is much better with regards to predictive capability.  My question is why waic for the poisson model is Inf.  Data is under-dispersed and therefore I know that the poisson model is not suited well.

Code is attached.

Thank you,
Linas

<gamma.stan><waic.r><poisson.stan>

[Andrew Gelman]

We discuss this in the (not quite released) Waic paper.  The short answer is I like to do it like this:

data {

  int N;

  int J;

  vector[N] y;

  matrix[N,J] X;

}

parameters {

  vector[J] b;

  real<lower=0> sigma;

}

model {

  y ~ normal(X*b, sigma);

}

generated quantities {

  vector[N] log_lik;

  for (n in 1:N){

    log_lik[n] <- normal_log(y[n], X[n]*b, sigma);

  }

}

I have two reasons for preferring the above form.  First, it has the advantage that the Waic-conputing code is just at the end, it’s just added to the end of an existing model and I don’t need to mess with the model itself.  Second, by not messing with the model, I can do the inference in vectorized form which is faster.  In any case I need to save log_lik as a vector not a scalar because Waic needs all the components separately.

A

On Apr 18, 2014, at 11:49 PM, Brian Hayden <bha...@lbl.gov> wrote:

What is the difference between doing log_lik in the transformed parameters block and using increment_log_prob in the model block, versus doing it in a generated quantities block like this:

generated quantities {
    real log_lik;
    log_lik <- 0;
    for (n in 1:N){
        log_lik <- log_lik + normal_log(y[n], X[n]*b, sigma);
    }
}


I am asking mainly because I am trying to implement waic and I want to make sure the generated quantities version doesn't do something different.

[Andrew Gelman]

Hi—how many data points do you have in your example?

On May 1, 2014, at 6:23 PM, Brian Hayden <bha...@lbl.gov> wrote:

I've been working on something else for the last couple of weeks, but in thinking about my stan-based project today, I had some other thoughts about WAIC. I am doing essentially a model comparison project, and I had been using the AIC as my metric of comparison. The stan "optimizing" function was not giving me reliable results from testing on simulated data; so I was getting the maximum likelihood estimate by simply taking way more samples than needed for estimating posterior statistics, and picking out the maximum likelihood from there. 

For WAIC, tracking all of the pointwise likelihoods becomes a memory problem very quickly. However, it occurs to me that with waic, one shouldn't need many more samples than what is required to obtain good statistics on the parameters of the model. Since waic is based on the mean and variance of each pointwise likelihood, when I am satisfied that the "error" in waic (simply due to the error of the mean and variance that goes into waic) is sufficiently low, I don't need more samples. Is this true?

-Brian

[Andrew Gelman]

400 points doesn’t seem like so much, I don’t see why this would overwhelm your memory?

On May 1, 2014, at 8:00 PM, Brian Hayden <bha...@lbl.gov> wrote:

There are ~400 data points with ~800 parameters defined in the parameter block. Each data point has certain measurements associated with it and we draw a "true" value for these measurements from a prior distribution (also fit by the model) each sample.

-Brian

[Andrew Gelman]

There is no way you should be doing 500,000 iterations in Stan (I think).  What is your effective sample size?

On May 1, 2014, at 8:14 PM, Brian Hayden <bha...@lbl.gov> wrote:

Oh, because I was sampling ~500,000 times in order to get a relatively reliable maximum likelihood estimate for the AIC. I am using pystan, and when I go to calculate WAIC, python is spitting out a memory error when calculating the mean and variance required.

There is a related issue, which is that in pystan, there is currently a bug where the user cannot tell stan to 'ignore' any parameters, so sampling from a large model runs out of memory pretty quickly (in this case the amount of memory available to me is ~3 GB on a compute node).

Regardless, my main question is: is it true that with WAIC I am only really concerned with how well I measure the mean and variance terms from BDA (lppd and pwaic2), which basically requires the same number of samples as is required to get good statistics on the posterior of a parameter?

-Brian