plan for variables
Bob Carpenter
language variable declarations and statement blocks:
--------------------------------
data {
<declarations>
}
transformed data {
<declarations>
<statements_no_sampling>
}
parameters {
<declarations>
}
transformed parameters {
<declarations>
<statements_no_sampling>
}
model {
<statements>
}
derived {
<declarations>
<statements_no_sampling>
}
-----------------------------------
<statements> is a sequence of statements, including
sampling statements (y ~ foo(theta)) and assignment
statments (y <- foo(theta)) and loops.
<statements_no_sampling> eliminates sampling statements.
In each block of statements, assignments are only
allowed to be made to local variables and variables
declared in the block.
Statement blocks allow local variables to be declared as
in a programming language. For example, we'll allow
{
double mu;
for (n in 1:N) {
mu <- x[n] * alpha + beta;
y[n] ~ normal(mu,1);
}
}
The log probability variable lp__ will only be
available in transformed parameters (to allow
Jacobians) and in the model.
(1) Input will read all declared data.
(2) sampling is defined over the declared parameters;
they're read by the model each leapfrog step.
(3) output will produce all of the parameters,
transformed parameters and derived quantities (but
not any local variables).
I know Andrew liked this idea, but I think it's a
really bad idea to have "y ~ foo(theta)" be ambiguous
and act as "y <- r_foo(theta)" in one context and
"lp__ += foo_log(y,theta)" in others.
- Bob
Andrew Gelman
A
Matt Hoffman
Matt
Andrew Gelman
"data" = input
"transformed data" = computed before compilation
"parameters" = updated every leapfrog step
"transformed parameters" = computations performed every leapfrog step
"model" = computations of log-posterior performed every leapfrog step
"derived" = computations and simulations performed at the end of every iteration
By the way, why can't we have sampling in the "derived" step? Unless I'm misunderstanding, that would seem to make sense to me. For example, that would be the place in the algorithm to do forward sampling of leaves and to simulate fake data, replicated data, etc.
A
On Dec 31, 2011, at 1:17 PM, Bob Carpenter wrote:
Bob Carpenter
> "transformed data" = computed before compilation
> "parameters" = updated every leapfrog step
> "transformed parameters" = computations performed every leapfrog step
> "model" = computations of log-posterior performed every leapfrog step
> "derived" = computations and simulations performed at the end of
every iteration
Right.
Also, every leapfrog step, constrained parameters have
their constraints and log absolute Jacobian determinants
applied. Once per sample, there's an additional constraining
of parameters for output.
Partial derivatives only get calculated for parameters.
> By the way, why can't we have sampling in the "derived" step? Unless
> I'm misunderstanding, that would seem to make sense to me. For
> example, that would be the place in the algorithm to do forward
> sampling of leaves and to simulate fake data, replicated data, etc.
As to no sampling in derived blocks, it's just a matter
of notation. I think it's too confusing to have
y ~ foo(theta)
mean
lp__ += log(d_foo(y,theta));
in the model and
y <- random_foo(theta);
in derived quantities. So users are now required to write
out the assignment (y <- random_foo(theta)) expicitly
for derived quantities. So we can do the forward sampling
to simulate data.
Eventually, I'd like to write a proper forward-sample
compiler. It would require all the upper parameters
to be fixed, so that every time we sampled, all the
parameters (in the usual sense) would be known. A simple
case is sampling leaves of a directed graphical model.
With the compiler, the user could use the same model code
to do the simulations.
I know you (Andrew) like the overloading of sampling notation
(y ~ foo(theta)). That's probably a CS/stats cultural
difference, as I find the overloading in standard stats
notation and in Bayesian practice very confusing.
I agree with you (Andrew) that writing "random_foo"
is tacky and that R/BUGS users will start writing "d_foo" and
"r_foo" by mistake. I'll at least be able to catch that with
the compiler and flag the location of the problem.
Not much I can do about the tackiness.
> These labels are as good as any, I suppose,
If you have better suggestions for the labels, bring
them on! I'm not happy with "parameters" -- it
matches BUGS's usage, but it doesn't match statistics
conventions where for normal(y|mu,sigma), the mu
and sigma are "parameters" and y is "data". Hierarchical
Bayes screws all this up, because you can
have normal(mu|mu_mu,sigma_mu) and normal(y|mu,sigma)
where mu is both a "parameter" and "data". I sort
of like "knowns" and "unknowns", but it just doesn't
seem to scan in the notation. Or "constants" and
"variables", but that's also at odds with statistical
convention because some of the "modeled data" are
treated like variables even though their value is
known.
- Bob
Andrew Gelman
Let's just avoid dnorm() and rnorm() entirely.
A
Bob Carpenter
sampling in the model is:
y ~ normal(m,s);
For generating random variates, it'll be:
y <- random_normal(m,s);
Would sim_normal(m,s) be clearer or sample_normal(m,s)?
The sampling notation compiles to:
lp__ += normal_log(y,m,s);
The user can refer to the function normal_log(y,m,s)
the same way as they refer to functions like sqrt(t).
They also have access to the lp__ accumulator (but
only in the model part), so they could write:
lp__ <- lp__ + normal_log(y,m,s);
which compiles to:
lp__ = lp__ + normal_log(y,m,s);
- Bob