"official" Eight Schools model

[Bob Carpenter]

Andrew:

Here's what I checked into Stan:

data {
int[0,] J; // number of schools
real y[J]; // estimated treatment effect (school j)
real[0,] sigma_y[J]; // std dev of effect estimate (school j)
}
parameters {
real theta[J];
real mu_theta;
real[0,1000] sigma_theta;
}
model {
mu_theta ~ normal(0,1000);
theta ~ normal(mu_theta, sigma_theta);
y ~ normal(theta,sigma_y);
}


It uses vectorization for the mu_theta, theta, and y.
If you would prefer, I can write the model out with loops:

for (j in 1:J)
theta[j] ~ normal(mu_theta, sigma_theta);
for (j in 1:J)
y[j] ~ normal(theta[j], sigma_y[j]);

Or with a single loop:

for (j in 1:J) {
theta[j] ~ normal(mu_theta, sigma_theta);
y[j] ~ normal(theta[j], sigma_y[j]);
}

I could also add in the redundant uniform sampling model
for sigma_theta, if you'd prefer:

sigma_theta ~ uniform(0,1000);

- Bob

[Jiqiang Guo]

I will change the wiki page to this version.  Essentially, in my version, I do not have prior statement for mu_theta and sigma_theta, as presented in chapter 5 of BDA2.  Here we have those diffuse prior for mu_theta and sigma_theta just because BUGS and JAGS must need a prior statement for parameters.  So I think it might be a bit better that we use 

real sigma_theta;

in the parameters block so that we could either have 

mu_theta ~ normal(0,1000);

sigma_theta ~ uniform(0, 1000);

or leave them out in the model block as Stan allow it. 

--

Jiqiang 

--

[Bob Carpenter]

I want to highlight something important before
getting back to Eight Schools below:

If we define sigma_theta to sample from a uniform

sigma_theta ~ uniform(0, 1000);

then it should be declared with matching bounds, as in

real[0,1000] sigma_theta;

It might work w/o bounds if you happen to get lucky and start
with a value in [0,1000] -- otherwise, you're in NaN land
outside of the support for sigma_theta. You can also land
outside the support due to the random momentum of each iteration,
which is always Normal(0,1), no matter which dimension it is.

On 8/13/12 9:58 PM, Jiqiang Guo wrote:
> I will change the wiki page to this version. Essentially, in my version, I do not have prior statement for mu_theta and
> sigma_theta, as presented in chapter 5 of BDA2.

I forgot about that one. So that's three
versions just in BDA.

I don't know where the Cauchy prior came from
because the ones in src/models were from the appendix
of BDA. I had played around with other priors like
the Cauchy and Gamma before (following another paper of
Andrew et al.'s on hierarchical variance priors).

> Here we have those diffuse prior for mu_theta and sigma_theta just
> because BUGS and JAGS must need a prior statement for parameters. So I think it might be a bit better that we use
>
> real sigma_theta;
>
> in the parameters block so that we could either have
>
> mu_theta ~ normal(0,1000);
> sigma_theta ~ uniform(0, 1000);
>
> or leave them out in the model block as Stan allow it.

Given that we need the declaration to be:

real sigma_theta[0,1000];

then it's just na issue of whether we include the redundant:

sigma_theta ~ uniform(0,1000);

It's not really even any less efficient because the implementation's
clever

[Jiqiang Guo]

Yes they need to be highlighted.  I forgot that already. Thanks.

Jiqiang 

[Andrew Gelman]

J

I prefer the version that was on the webpage, that you already had before, using the notation mu, tau, and with no prior on tau.

A

--
 
 

[Bob Carpenter]

n 8/14/12 11:07 AM, Andrew Gelman wrote:
> J
> I prefer the version that was on the webpage, that you already had before, using

I'm not sure what web page you're talking about.

Given that you (Andrew) seem to care the most about
which version of Eight Schools we use, could you
perhaps write it out in Stan notation for us?
When you get back we can give you a tutorial on
Git and our web pages and get you checking things in
yourself so that we can cut out the middleman.

> the notation mu, tau, and with no prior
> on tau.

To get started, here's a repeat of what I last
sent with un-suffixed variable names and with tau
interpreted as a precision as in the model in the appendix
of BDA Andrew said we should follow.

data {
int[0,] J; // number of schools
real y[J]; // estimated treatment effect (school j)

real[0,] sigma[J]; // std dev of effect estimate (school j)
}
parameters {
real theta[J];
real mu;
real[0,1000] tau;
}
model {
mu ~ normal(0,1000);
theta ~ normal(mu, 1/sqrt(tau)); // ??? precision to scale
y ~ normal(theta,sigma);
}

This is the vectorized version. Unvectorized, you'd replace the
model body with:

mu ~ normal(0,1000);

for (j in 1:J) {

theta ~ normal(mu, 1/sqrt(tau));
y ~ normal(theta,sigma);
}

The reason I wrote it as I originally did:

data {
int[0,] J; // number of schools
real y[J]; // estimated treatment effect (school j)
real[0,] sigma_y[J]; // std dev of effect estimate (school j)
}
parameters {
real theta[J];

real mu;

real[0,1000] sigma_theta;
}
model {
mu_theta ~ normal(0,1000);

theta ~ normal(mu, sigma_theta);
y ~ normal(theta,sigma_y);
}

is that I was trying to be consistent in how we
named variables in our models. It helps a lot with
understanding programs and hence coding them and
maintaining them. For instance, the lack of such a
convention is why I can't tell if the tau you want
us to use is a scale or a precision. Varying notation
and terminology is especially hard on beginners.

BDA/ARM introduced a convention of subscripting location
and scale/precision parameters with the name of the variate.
BDA/ARM also used sigma for scale parameters and tau
for precisions, though I believe in some cases tau were also
used for scale parameters, which is why I don't know which
one to use just based on the name "tau".

Are there other coding conventions? I'll add an appendix
to the manual at some point.

- Bob

[Jiqiang Guo]

It is just  http://code.google.com/p/stan/wiki/RStanGettingStarted 

P.S. I have not changed it and I won't change it.

Jiqiang