Re: [stan-users] 'Matt trick' for logistic-binomial model

[Bob Carpenter]

A serious problem with your model is that you declare
sigma2 as having support (0,infinity)

real<lower=0> sigma2;

and then you use it with bounded support:

sigma2 ~ uniform(0,100);

What you need to do is declare sigma2 with support
that matches the model.

real<lower=0,upper=100> sigma2;

The re-encoding in model2 forgot the additive term for
the mean. Instead of

transformed parameters {
real u[G];
for (g in 1:G)
u[g] <- e_u[g] * sqrt(sigma2);
}

you need the assignment to include the mean u0:

for (g in 1:G)
u[g] <- u0 + e_u[g] * sqrt(sigma2);

but you really want to declare e_u as a vector:

vector[G] e_u;

and then vectorize the assignment:

u <- u0 + e_u * sqrt(sigma2);

In fact, for both speed and simplicity, I'd replace
prior for variance

sigma2 ~ uniform(0,100);

with a uniform prior on standard deviation (changing the
constraint on its support to upper=10, of course!):

sigma ~ uniform(0,10);

and then replace the vectorized assignment with

u <- u0 + e_u * sigma;

- Bob




On 6/7/13 5:50 AM, Adam Slez wrote:
> Hi All,
>
> I'm just getting started with rstan and I'm trying to implement a hierarchical logistic-binomial model. In an effort to
> speed things up a bit, I've been trying to incorporate the 'Matt trick' described here
> <https://groups.google.com/forum/#!msg/stan-users/4gv3fNCqSNk/J6ZItL2ZJ-IJ> (among other places). I'm not sure whether
> I'm missing something, but when I use the Matt trick I find that the estimates of the intercept/location and
> variance/scale^2 are way out of line with what I would expect given (a) results from glmer and (b) results from rstan
> based on the conventional parameterization of the model (i.e. without the Matt trick).
>
> Here are the results from glmer:
>
> Generalized linear mixed model fit by the Laplace approximation
> Formula: cbind(y, trials - y) ~ x + (1 | group)
> Data: dat
> AIC BIC logLik deviance
> 136717 136730 -68356 136711
> Random effects:
> Groups Name Variance Std.Dev.
> group (Intercept) 0.80362 0.89645
> Number of obs: 574, groups: group, 144
>
> Fixed effects:
> Estimate Std. Error z value Pr(>|z|)
> (Intercept) -2.7816015 0.0782625 -35.54 <2e-16 ***
> x 0.0591794 0.0006591 89.79 <2e-16 ***
> ---
> Signif. codes: 0 �***� 0.001 �**� 0.01 �*� 0.05 �.� 0.1 � � 1
>
> When I run the same model in rstan without using the Matt trick, the corresponding estimates are as follows:
>
> (Intercept) = -2.7816
> x = 0.0592
> Variance = 0.8036
>
> These are clearly in line with the results from glmer. When I run the same model with the Matt trick, I get the following:
>
> (Intercept) = 1.2287
> x = 0.0594
> Variance = 9.8941
>
> I used 4 chains with 5000 iterations apiece in each model. I've attached files of the output, along with copies of the
> code. Model 1 uses the conventional parameterization while model 2 uses the Matt trick.
>
> Interestingly, I had no problem getting the Matt trick to work when using a Bernoulli distribution to model some other
> (ungrouped) data. What am I missing about this case? I noticed that the effective sample size for the intercept and
> variance terms was pretty low when using the Matt trick. I was kind of surprised given that this is something that the
> Matt trick should help with. Rerunning the Matt trick model with an additional 5000 iterations didn't change the
> results. Of course, it is entirely possible that I just goofed up the code somewhere along the line.
>
> On a related note, what are the expectations with respect to the transformed parameters (e.g., the e_u terms in the
> attached ouput) which do not appear in the model? In particular, do we worry about the effective sample size for these
> things?
>
> Any help would be much appreciated.
>
> Adam
>
> --
> 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]

I don't think uniform(0,100) is such a good model. I'd prefer something smooth, also I'd prefer modeling sigma, not sigma^2.
A

>> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

[Adam Slez]

Bob and Andrew,

Thanks for the helpful suggestions.  Right after I posted I realized that I had forgotten to incude the u0 term in the transformed parameters block.  I tried to pull the post once I realized my mistake.  After I put the mean term back in, things worked more or less as expected in terms of the estimates from the posterior.

If you don't mind, I have a couple of follow up questions...

1.  Andrew, is Bob's suggestion of using uniform(0,10) for sigma more in line with what you had in mind or were you thinking of something else?

2.  The effective sample size on all of the terms involved with transformed parameter block are surprisingly low.  I've attached the code (updated based on your suggestions) and results from a run using 4 chains with 10K iterations apiece.  Any thoughts on what is going on here?  While I'm definitely concerned about the effective sample size associated with sigma and u0, is there any need to worry about the e_u terms, or can we safely set these aside as a byproduct of the Matt trick. In short, do we expect viable sample sizes for things like the e_u terms?

3.  I ran into the Matt trick while looking into ways of speeding up the model.  One possibility is to run this in JAGS which has a dedicated glm module.  In discussing this point in another thread, Matt Hoffman alluded to the possibility of using the Matt trick in JAGS.  I decided to check this out, but started getting totally off-the-wall estimates for sigma.  As suggested by my original post, it is totally possible that I just goofed up the code (see attached).  If you have any thoughts, I'd love to hear them, though you should feel free to tell me to take a hike given that this is a Stan list and I can understand that you might not want to accidentally encourage possible defection, even if only in the very short term.  Even if you don't feel comfortable commenting on the JAGS code, is there any reason to believe that the Matt trick would be viable in one program but another?  I could imagine situations in which it didn't do much with respect to speed, but I have a tougher time seeing how it could actually produce nonsensical results.

Thanks again for all of your help.  All of these questions aside, Stan is a fine piece of software and I look forward to seeing where it goes.

Adam 

[Andrew Gelman]

Adam

You can do uniform(0,10) on sigma but gneerally I'd prefer something smoother that drops off gradually, something like a half-cauchy.

Also, I can't see it making sense to do things in Jags. If Stan is running much slower (in time per effective sample size) that Jags, this is a problem we should be able to fix.  

A

<stanout_new.txt><model2.jags><model2.stan>

[Adam Slez]

Andrew,

I've attached a table which includes diagnostics from a series of models, two of which were run in Stan, one of which was run in JAGS, all called via R.  The first Stan model uses a conventional random intercept parameterization, while the second uses the Matt trick.  The JAGS model uses the conventional parameterization because, as mentioned before, the Matt trick code I wrote produces crazy estimates for sigma.  These estimates were pulled straight off the screen which means there may be a bit of an apples-and-oranges problem here when it comes to trying to compare effective sample size and R-hat between Stan and JAGS.  If you think this is a serious problem for the purposes of this discussion, I'm happy to go back and recalculate things.  Looking at these rough estimates, however, not only does JAGS seem to outpace Stan, but Stan without the Matt trick seems to better than Stan with the Matt trick.

Just to be clear about a couple of things...

1.  Times for Stan do not include the time taken to compile the models.

2.  I made sure to use set_cppo("fast"), though this type of thing wouldn't explain the divergence between alternative parameterizations in Stan.

What else could I be missing here?

Adam

[Bob Carpenter]

On 6/7/13 6:35 PM, Adam Slez wrote:
...

> 2. The effective sample size on all of the terms involved with transformed parameter block are surprisingly low. I've
> attached the code (updated based on your suggestions) and results from a run using 4 chains with 10K iterations apiece.

> Any thoughts on what is going on here? While I'm definitely concerned about the effective sample size associated with
> sigma and u0, is there any need to worry about the e_u terms, or can we safely set these aside as a byproduct of the
> Matt trick. In short, do we expect viable sample sizes for things like the e_u terms?

Are you doing this through bin/print or using RStan?

I'm curious about what's going on here conceptually.

With the Matt trick, you sample

e_u ~ normal(0,1)

and then set:

u <- u0 + sigma * e_u

u0, sigma, and e_u all have low n_eff (< 70)
whereas u has a high n_eff (2,000).

Is it the extra variability of the three variables combined
that gives u relatively low autocorrelation?

It also looks like some of the u0 have too high R-hat (1.16),
which I wouldn't expect after so many iterations.
So it's possible one of the chains just got stuck with your
initialization. You can try init=0, which is more like what
you get with JAGS if you don't initialize yourself. Could
you also send the traceplot for u0, bx, and sigma?

Do the traceplots show high correlation of the u0, sigma,
and e_u variables? Are the posterior variances the same
in Stan and JAGS?

Andrew?

> 3. I ran into the Matt trick while looking into ways of speeding up the model. One possibility is to run this in JAGS
> which has a dedicated glm module. In discussing this point in another thread

> <https://groups.google.com/forum/?fromgroups#!topic/stan-users/QDumPPbRa4Q>, Matt Hoffman alluded to the possibility of

> using the Matt trick in JAGS. I decided to check this out, but started getting totally off-the-wall estimates for
> sigma.

I'd check in with Martyn Plummer on the JAGS list to see
if JAGS will sample properly through the assignments.
I don't see why not, but I don't know for sure it can
calculate the right Markov blanket.

I'd also check if the Matt trick will interfere with the GLM.

For Stan, there's nothing special done for GLMs and the Matt
trick works because of the way Stan works by just defining
the log probability function.

> As suggested by my original post, it is totally possible that I just goofed up the code (see attached).

It looks OK to me. And I don't see an obvious way to
make it more efficient.

None of the scales look that crazy, so I'm also wondering if
many of the predictors are correlated.

> I could imagine situations in which it [the Matt trick] didn't do much with

> respect to speed, but I have a tougher time seeing how it could actually produce nonsensical results.

There are so many ways things can go wrong and produce
nonsensical results.

- Bob

[Ben Goodrich]

On Friday, June 7, 2013 10:39:53 PM UTC-4, Bob Carpenter wrote:

I'm curious about what's going on here conceptually.

With the Matt trick, you sample

   e_u ~ normal(0,1)

and then set:

   u <- u0 + sigma * e_u

u0, sigma, and e_u all have low n_eff (< 70)
whereas u has a high n_eff (2,000).

Is it the extra variability of the three variables combined
that gives u relatively low autocorrelation?

This is odd. I would be interested to see the joint distribution of u0 and sigma, which can be done via rstan with something like

pairs(model2, pars = c("u0", "sigma"), chain_id = c(1,2))

in this case the results for chain 1 will be in the lower triangle and chain 2 will be in the upper triangle. I would bet at least one of the chains is messed up.

Ben

[Adam Slez]

Bob, Ben, and Andrew,

I redid the analysis used to produce the previous speed tests which, I forgot to mention, were only based on 2K iterations (hence why some of the R-hats were so high).  Here's what is going on in the setup of the new tests:

1.  I used 4 chains with 10K iterations.

2.  I estimated three additional models.  As per Bob's suggestion, the first uses the Matt trick as before, but initial values are set to zero.  The second uses the Matt trick, but the x variable has been mean centered via the use of a transformed data block.  The third model also uses a mean-centered x, but does not use the Matt trick.

3.  Output was produced using the monitor command so as to ensure comparability between Stan and JAGS, as per the discussion here.  This speaks to Bob's question regarding the calculation of effective sample size.  One thing I noticed is that the monitor results did not match the print results.  I might expect this when looking at the JAGS output, but this also happened in Stan as well.  My understanding is that the print output in Stan is based on the summary function.  How do the summary calculations differ from the monitor calculations?

Here's my reading of the findings:

1.  JAGS outpaces Stan by quite a bit, regardless of how the model is specified in Stan.

2.  Stan models which use the Matt trick seem to do worse than Stan models without the Matt trick.

3.  Models with and without the Matt trick differ in terms of which parameters they are able to efficiently estimate.  Models with the Matt trick do better with respect to bx as compared to u0 or sigma.  The opposite holds for models without the Matt trick.

4.  Mean centering doesn't seem to help much.

5.  To speak to Bob's earlier question, in terms of posterior variance Stan actually seems to be doing slightly better when we look at the standard deviation.  Things are obviously a bit different when we look at, say, the standard error of the mean which corrects for effective sample size. Overall, JAGS seems to be doing better here, in large part because of the fact that Stan does well with either u0 and sigma or bx, but not both.  More specifically, we see Stan without the Matt trick doing better than JAGS with to the standard error of the mean for u0 and sigma.  In contrast, if we either center x or set initial values to 0, we see Stan with the Matt trick doing better than JAGS with respect to the standard error of the mean for bx.

I have also attached the diagnostics requested by Bob and Ben.  These include both a traceplot for u0, sigma, and bx, as well as the pairs graph depicting chain-specific results (chains 1 and 2 are on the first page, while chains 3 and 4 are on the second).  To show the correlation between parameter estimates, I have also included a scatterplot matrix depicting draws for u0, sigma, bx, and five randomly chosen group effects (both transformed and untransformed). I'm not sure whether this is the best approach for answering Bob's question as to whether the traceplots show high correlation of the u0, sigma, and e_u variables.  In the first place, I'm not sure that what I did even counts as a traceplot in the usual sense. It would also be preferable to look at all of the u and e_u terms, but I couldn't think of a good way to do this given that I have 144 groups and, hence, 288 + 3 = 291 parameters to graph.  All graphs correspond to the row in the table labeled "Stan w/Matt" (i.e. I used the Matt trick but did not center x or force initial values to 0).

My read is that the chains for u0 and sigma aren't mixing especially well (at least not compared to bx) and that there is some correlation between sigma, u, and e_u, though maybe not a crazy amount.  The u terms seem okay, there is some evidence of correlation between sigma and the e_u terms, and there is definitely correlation among the e_u terms themselves.  Again, this approach is less than perfect because I'm only looking at a handful of u/e_u terms.  Things could be worse elsewhere in the model.

Adam

[Andrew Gelman]

Hmm, this is disturbing.  We should figure out what is going on here.

A

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

<new_comparison.xls><pairs.pdf><scatter.jpg><trace.pdf>

[Andrew Gelman]

Adam

Could you please send again your stan code, jags code, and r code?  I'm not sure how you called monitor after running the stan code; it's possible that by accident it threw away half the simulations.

Andrew

On Jun 9, 2013, at 8:59 AM, Adam Slez wrote:

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

<new_comparison.xls><pairs.pdf><scatter.jpg><trace.pdf>

[Adam Slez]

Andrew,

Everything is attached.  I apologize for not labeling the models in the order that they were run, but you should be able to see what's what from the R code.  

Good call on losing simulations.  That makes a lot of sense given that we see things maxing out at 10K effective samples.  With 4 chains at 10K iterations a piece and 5K warmup runs per chain, we'd expect to get 20K simulations under completely independent sampling. Obviously we don't usually expect perfectly independent draws, but I think you can see where I'm going here.  Outside of getting the maximum 20K, getting something like an even 10K is a little suspicious.  

I just looked at output produced via print and I'm getting effective sample sizes of 20K in some places which would explain the difference I saw between print and monitor.  Do you think the monitor issue alone is enough to explain the apparent disparities between Stan and JAGS?  I was seeing similar results when using the print command, though I don't think the sample size calculations were perfectly comparable in that context.

Adam

[Adam Slez]

Andrew,

Now that you have all the code in front of you, you'll see that I'm still using uniform priors for the variance components.  I want to be clear that I'm definitely not ignoring your advice on the half-Cauchy.  I'm still getting up to speed on how to implement this, so I decided to set it aside for the moment.  My thought was that the uniform prior would still be preferable to, say, an inverse-gamma.  I just ran across the half-Cauchy code on p.430 of Gelman and Hill (as well as the nice discussion of multiplicative redundant parameters on p.424-425), so I don't see any problem working this in very soon.  

My understanding is that the half-Cauchy is even easier in Stan, in that you can safely ignore the half part of things, recognizing, of course, that calling the Cauchy distribution directly might not be the most computationally efficient approach.

This is obviously a separate issue from the speed and monitor concerns.  I should probably stop hijacking my own thread...

Adam

[Bob Carpenter]

Cauchy's just like any other prior. It's in the manual
section on regression and also documented in the distributions.
In Stan:

parameters {
real<lower=0> sigma;
...
model {
sigma ~ cauchy(0,5);
...

In JAGS, you just need to truncate, so it'd be this
if they have the cauchy

sigma ~ dcauchy(0,5) T[0,];

and otherwise you'll need to replace dcauchy with the
T distribution with 1 degree of freedom.

You need to truncate for the uniform in JAGS, too.

- Bob

>> <javascript:>.
>> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.
>>
>>
>> <new_comparison.xls><pairs.pdf><scatter.jpg><trace.pdf>

[Bob Carpenter]

Thanks for all the careful evaluation work! More below.

On 6/9/13 8:59 AM, Adam Slez wrote:
> Bob, Ben, and Andrew,
>
> I redid the analysis used to produce the previous speed tests which, I forgot to mention, were only based on 2K
> iterations (hence why some of the R-hats were so high). Here's what is going on in the setup of the new tests:
>
> 1. I used 4 chains with 10K iterations.

It's really too bad that many iterations are necessary, but
I can see in the traceplots some parameters aren't mixing well.

> 2. I estimated three additional models. As per Bob's suggestion, the first uses the Matt trick as before, but initial
> values are set to zero. The second uses the Matt trick, but the x variable has been mean centered via the use of a
> transformed data block. The third model also uses a mean-centered x, but does not use the Matt trick.

The traceplots show pretty clearly that this isn't an initialization problem.

> 3. Output was produced using the monitor command so as to ensure comparability between Stan and JAGS, as per the
> discussion here. This speaks to Bob's question regarding the calculation of effective sample size. One thing I noticed
> is that the monitor results did not match the print results. I might expect this when looking at the JAGS output, but
> this also happened in Stan as well. My understanding is that the print output in Stan is based on the summary function.
> How do the summary calculations differ from the monitor calculations?

Very good question. Jiqiang and Daniel should know.

One issue is going to be whether the warmup iterations get included.
They should NOT be included in estimating means and variances and
MCMC error.

> Here's my reading of the findings:
>
> 1. JAGS outpaces Stan by quite a bit, regardless of how the model is specified in Stan.

Can you quantify "quite a bit"?

And you're sure that Stan itself was compiled using -O3? It
should be automatic if you're using the command-line version
and just call "make". RStan's a bit more probelmatic because your
current R config determines how binary packages get compiled.
We're trying to make this more transparent for the next release.

> 2. Stan models which use the Matt trick seem to do worse than Stan models without the Matt trick.

I have no idea what's going on here, as this isn't what we've
found in models like this that we've run.

> 3. Models with and without the Matt trick differ in terms of which parameters they are able to efficiently estimate.
> Models with the Matt trick do better with respect to bx as compared to u0 or sigma. The opposite holds for models
> without the Matt trick.

That's very interesting. The traceplot you sent shows bx mixing well, but u0 and
sigma not mixing well, so that must've been with the transform.

> 4. Mean centering doesn't seem to help much.

It shouldn't really help much at all given the refactorization.

> 5. To speak to Bob's earlier question, in terms of posterior variance Stan actually seems to be doing slightly better
> when we look at the standard deviation.

Better in that it's estimate is closer to the true values?

Was this simulated data you used or real data?

> Things are obviously a bit different when we look at, say, the standard error
> of the mean which corrects for effective sample size.

I'm not sure what you mean here. I was talking about estimating the variance,
so it'd be the standard error of the variance estimate we'd care about.

If the variance itself is overestimated or underestimated, it'll throw
off MCMC error by the same amount, of course.

The problem that can show up with either Gibbs or HMC is that the
sampler won't sample into the tails of the posterior, so variance
gets underestimated.

> Overall, JAGS seems to be doing better here, in large part because
> of the fact that Stan does well with either u0 and sigma or bx, but not both. More specifically, we see Stan without
> the Matt trick doing better than JAGS with to the standard error of the mean for u0 and sigma. In contrast, if we
> either center x or set initial values to 0, we see Stan with the Matt trick doing better than JAGS with respect to the
> standard error of the mean for bx.

> I have also attached the diagnostics requested by Bob and Ben. These include both a traceplot for u0, sigma, and bx, as
> well as the pairs graph depicting chain-specific results (chains 1 and 2 are on the first page, while chains 3 and 4 are
> on the second). To show the correlation between parameter estimates, I have also included a scatterplot matrix
> depicting draws for u0, sigma, bx, and five randomly chosen group effects (both transformed and untransformed).

The diagnostics look great. Surprisingly, there seems to be very little posterior
correlation among u0, sigma, and bx.

The blown up pairs.pdf indicates that we perhaps haven't even
fully covered the posterior density properly after 10K samples,
especially looking at u0 vs. bx. But I'm confused as to why
the u0 x sigma scatterplot doesn't look like the transpose of the
sigma x u0 scatterplot.

On the other hand, the coefficients themselves are highly correlated not
only with u0, but with each other.

This problem's small enough that you could deal with this kind of systematic
correlation using a full mass matrix. (It runs slower per iteration typically,
but can mix better.)

> I'm not
> sure whether this is the best approach for answering Bob's question as to whether the traceplots show high correlation
> of the u0, sigma, and e_u variables.

I was thinking more of just calculating the correlation in
the extracted values.

> In the first place, I'm not sure that what I did even counts as a traceplot in the
> usual sense. It would also be preferable to look at all of the u and e_u terms, but I couldn't think of a good way to do
> this given that I have 144 groups and, hence, 288 + 3 = 291 parameters to graph. All graphs correspond to the row in
> the table labeled "Stan w/Matt" (i.e. I used the Matt trick but did not center x or force initial values to 0).
>
> My read is that the chains for u0 and sigma aren't mixing especially well (at least not compared to bx) and that there
> is some correlation between sigma, u, and e_u, though maybe not a crazy amount.

What's surprising is that's not what it looks like in your scatterplot.

> The u terms seem okay, there is some
> evidence of correlation between sigma and the e_u terms, and there is definitely correlation among the e_u terms
> themselves. Again, this approach is less than perfect because I'm only looking at a handful of u/e_u terms. Things
> could be worse elsewhere in the model.

I'd guess that other u terms would look similar.

- Bob

[Ben Goodrich]

On Monday, June 10, 2013 10:44:31 AM UTC-4, Bob Carpenter wrote:

The blown up pairs.pdf indicates that we perhaps haven't even
fully covered the posterior density properly after 10K samples,
especially looking at u0 vs. bx.   But I'm confused as to why
the u0 x sigma scatterplot doesn't look like the transpose of the
sigma x u0 scatterplot.

TFM says for the chain_id argument to pairs():

By default, NULL, which will plot roughly half of the chains in the lower panel and the rest in the upper panel. Otherwise, an integer vector can be passed to select some subset of the chains, of which roughly half will be plotted in the lower panel and the rest in the upper panel. Or a list of two integer vectors, each specifying a subset of the chains to be plotted in the lower and upper panels respectively. The latter option can be useful if you suspect that one chain has failed to converge.

So, the i-j plot will never be identical to the j-i plot, unless you do something wasteful like chain_id = c(1,1). I'm not sure exactly what was passed for this plot, but if the i-j plot looks substantially different than the j-i plot, then it didn't converge.

Ben

[Adam Slez]

Bob,

Some follow-up comments....

3.  Output was produced using the monitor command so as to ensure comparability between Stan and JAGS, as per the
discussion here.  This speaks to Bob's question regarding the calculation of effective sample size.  One thing I noticed
is that the monitor results did not match the print results.  I might expect this when looking at the JAGS output, but
this also happened in Stan as well.  My understanding is that the print output in Stan is based on the summary function.
  How do the summary calculations differ from the monitor calculations?

Very good question.  Jiqiang and Daniel should know.

I think this may have already been answered by Andrew.  Monitor seems to be dropping half of the post-warmup iterations.  This does not seem to be the case with print.

Here's my reading of the findings:

1.  JAGS outpaces Stan by quite a bit, regardless of how the model is specified in Stan.

Can you quantify "quite a bit"?

I'm just looking at the effective samples per second calculations I included in the model comparison table ("new_comparison.xls").  In terms of u0 and sigma, JAGS appears to be producing at least twice as many effective samples per second.  In terms of bx, JAGS appears to be doing at least seven times as many effective samples per second.  This doesn't account for the monitor issue.  Having looked at the print output as well, I'm not sure correcting the monitor problem is going to radically change this outcome.

 

And you're sure that Stan itself was compiled using -O3?  It
should be automatic if you're using the command-line version
and just call "make".  RStan's a bit more probelmatic because your
current R config determines how binary packages get compiled.
We're trying to make this more transparent for the next release.

I'll be honest and say I'm not sure.  I've been running everything using RStan which I installed using the steps shown on the Getting Started page.  For the most part I've been using 64-bit Rstan on a MacBook Pro with set_cppo(mode = "fast").  I've also tried running things on a 32-bit Windows machine, though that seemed to do worse, not to mention the fact that I get crazy long warnings starting with MSDOS paths being detected and ending with 'void stan::agrad::free_memory()' defined but not used.  Though everything runs exactly the same despite these warnings, I've mostly stuck to running things on the Mac as a result.

 

5.  To speak to Bob's earlier question, in terms of posterior variance Stan actually seems to be doing slightly better
when we look at the standard deviation.

Better in that it's estimate is closer to the true values?

No.  I just meant that the estimated standard deviation of the posterior (as reported via monitor) was smaller in the Stan output than in the JAGS output.  I wasn't sure what you are asking for and this was my best guess.  Based on your original email and your response above, I see that you were looking for variance as opposed to the standard deviation.  I was being a bit lazy and working with the default output.  It wouldn't be too hard to go back and calculate the variance instead, assuming that I correctly understand what you are asking for.

Was this simulated data you used or real data?

This is real data.

Things are obviously a bit different when we look at, say, the standard error
of the mean which corrects for effective sample size.

I'm not sure what you mean here.  I was talking about estimating the variance,
so it'd be the standard error of the variance estimate we'd care about.

Again, I wasn't totally sure what you were looking for and this was a bad, nonsensical guess on my part. My thought was basically that the standard of the mean is going to be a function of the standard deviation of the posterior and, hence, might tell us about how spread out things are (assuming that a high standard deviation might indicate that a chain has run off in the opposite direction), though I suspect that this is not the right way to think about the problem given that similarities and differences in the standard error of the mean could be driven the effective sample size as opposed to the standard deviation.

Adam  

[Bob Carpenter]

I've attached:

-- Stan model (basic parameterization)
-- JAGS model (matching parameterization)
-- fake data simulator and fitter
The models are basically Adam Selz's with a bit of reformatting.
Here's the results:
JAGS Stan

time 30s 30s(*)
n_eff(**)
u0 180 4000
bx 90 4000
sigma 290 4000
R-hat(**)
max (hier) 1.04 1.000
low (low) 2.00 1.000
---------------------
(*) Stan time does not include model compilation.
(**) n_eff and R-hat calculated using R2jags for JAGS and
RStan for Stan; Stan's should be more conservative
Both systems estimated posterior means that matched the
simulated values well.
I ran a few times and the pattern was consistent.
Adam -- would you mind sending me your data so I
can try that? The above is with data simulated from
the model.
- Bob

>

[Bob Carpenter]

> I'll be honest and say I'm not sure. I've been running everything using RStan which I installed using the steps shown
> on the Getting Started page. For the most part I've been using 64-bit Rstan on a MacBook Pro with set_cppo(mode =
> "fast").

The problem is that you can only call set_cppo() after Stan's
installed. What I was asking about is what the situation was when
Stan was installed. Or you can reinstall after using set_cppo(), which
changes the R binary install configuration.

> I've also tried running things on a 32-bit Windows machine, though that seemed to do worse, not to mention the
> fact that I get crazy long warnings starting with MSDOS paths being detected and ending with 'void
> stan::agrad::free_memory()' defined but not used.

These are harmless warnings we've not yet fully suppressed. We need
to suppress them because they're clearly confusing our users.

> Though everything runs exactly the same despite these warnings, I've
> mostly stuck to running things on the Mac as a result.

Stan should work the same on the Mac and on Windows and on
Linux.

> 5. To speak to Bob's earlier question, in terms of posterior variance Stan actually seems to be doing slightly
> better
> when we look at the standard deviation.
>
>
> Better in that it's estimate is closer to the true values?
>
>
> No. I just meant that the estimated standard deviation of the posterior (as reported via monitor) was smaller in the
> Stan output than in the JAGS output.

We want to know what the right answer is for the
deviation (or variance). We estimate MCMC error by dividing
deviation by sqrt of the effective sample size.

> This is real data.

OK. I'll run some simulations and see how things look.

- Bob

[Adam Slez]

Bob,

More on the point below...

I'll be honest and say I'm not sure.  I've been running everything using RStan which I installed using the steps shown
on the Getting Started page.  For the most part I've been using 64-bit Rstan on a MacBook Pro with set_cppo(mode =
"fast").

The problem is that you can only call set_cppo() after Stan's
installed.  What I was asking about is what the situation was when
Stan was installed.  Or you can reinstall after using set_cppo(), which
changes the R binary install configuration.

I followed the instructions on the Getting Start page for reinstalling Stan after setting set_cppo() and found no appreciable difference in speed.

Adam

[Jiqiang Guo]

If following the getting started wiki, rstan would be installed under O3 mostly.  It is OK that set_cppo is not used the first time of installing rstan when set_cppo is not available; it uses O3 set in the package's Makefile.  For installing rstan the first time, other optimization levels other than O3 might be used, but it should be very rare.

Jiqiang