ice example

[Ben Goodrich]

The top of src/models/bugs_examples/vol2/ice/ice.stan states:

# http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/Vol2.pdf
# Page 29: Ice: non-parametric smoothing in an age-cohort model

// The model is the same specified here as in JAGS example,
// but this example's JAGS version is a different from that of
// WinBUGS because JAGS does not support some feature.
// TODO:
// Maybe we should check if we could specify
// the model the same as WinBUGS?

// status: the results are farely cose to those
// from JAGS

As it turns out, the models that I have been trying to estimate for
Wawro are based on the ice example in WinBUGS and we ended up trying
to do it in Stan specifically because JAGS can't estimate that model.

So, I tried to replicate the ice example from WinBUGS, which I can
essentially do. But WinBUGS is getting a (mean) posterior standard
deviation of 0.05 whereas I am getting about 0.32 even with tons of
iterations. If anyone has time, can you check to see if there is
something wrong with vanilla_ice.stan or has an idea as to why our
results are less precise. It is revision 594.

Thanks,
Ben

[Bob Carpenter]

I'll check after I get done chuckling at the name,
"vanilla_ice.stan". Now I'm dealing with the cognitive
dissonance from vanilla_ice being more complex
than ice.

The doc in ice.stan says the results are fairly close
to those from JAGS. This may be, but they'll probably
be even closer when sigma is declared double(0,1)
rather than double(0,) -- it's sampled sigma ~ uniform(0,1),
so you need the constraint.

Also, I'm not sure if normal(0,1000)
was intended -- we're using deviation, not
variance, as the parameter, so this is pretty
flat! But then BUGS uses 1E-6, so this is
probably what you meant.

And why that crazy sigma * 1E3 scale param for
beta[1] and beta[2]. Does it really relate to
sigma? Why not just draw beta[1] and beta[2]
from something relatively flat? Is it that we
expect sigma to be close to 1 so this actually
does something to the inferences?

vanilla_ice.stan
----------------

This model's really confusing (like Wawro's!).

I don't see a sigma in your model at all. Is
it still going to be 1/sqrt(tau)? You can define
a "generated quantities" block to generate things
like this (more efficient than adding a derived parameter,
which gets calculated every leapfrog step within HMC)_.

And while I see tau getting sampled from a gamma(shape,d),
I don't see it touching any other data. It's used to
define betasd[k], but then betasd[k] isn't used
anywhere else. So I'd expect tau to just be
distributed as gamma(shape,d) in the posterior.

I doubt that's what you intended, because in
the BUGS model, tau is also used in defining
betaprec[k], which is then used as the precision
parameter for sampling beta[k].

- Bob

[Bob Carpenter]

Nice!

- Bob

On 1/26/12 5:00 PM, Daniel Lee wrote:
> Bob, I just wanted to let you know that I made a make target that will
> run the model through the gm, compile, and run it. When we can build
> all our models, that was the target I was going to use for Jenkins.
>
> make models/bugs_examples/vol2/ice/vanilla_ice.stan

[Ben Goodrich]

On Jan 26, 4:56 pm, Bob Carpenter <c...@alias-i.com> wrote:
> I'll check after I get done chuckling at the name,
> "vanilla_ice.stan".  Now I'm dealing with the cognitive
> dissonance from vanilla_ice being more complex
> than ice.

vanilla_ice is a very complex phenomenon. I haven't actually looked at
the JAGS example in ice.stan very carefully, so maybe Maverick can say
more about that.

I did fix the prior on the betas like you suggested and now am getting
the right sigma (which is just 1 / sqrt(tau)) although these log Risk
Ratios are still not super close to what WinBUGS produces. Wawro is
going to try to run the ice example in WinBUGS but monitor all the
parameters so that we can get a closer comparison.

Thanks,
Ben

> > #http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/Vol2.pdf

[Jiqiang Guo]

As I coded that example, before looking at the details, I could talk about some 

that might help. 

The feature that JAGS dose not support but Openbugs does is cycles. 

This happened again in vol3/jama example, in which actually I only writes 

the following now:

# http://www.openbugs.info/Examples/Jama.html
# Jama River Valley Ecuador - Radiocarbon calibration with phase information

# It turns out that the model specified in OpenBUGS
# has *cycles*, which is not support by Stan (and JAGS).
# In section `Directed cycles` of JAGS manual, it writes
#
# Directed cycles are forbidden in JAGS. There are two important instances
# where directed
# cycles are used in BUGS.
# * Defining autoregressive priors
# * Defining ordered priors

So my understanding it that in BUGS, we can specify model, for example, 

## a bi-variate normal for (x1, x2)'

## mu <- c(2, 3); 

## Sigma <- matrix(c(40, 10, 10, 60), ncol = 2) 

## (x1, x2)' ~ mnorm(mu, Sigma) 

model {

  x1 ~ dnorm(mu1, tau1);

  mu1 <- 2 + (x2 - 3) / 6;

  tau1 <- 3 / 115;

  x2 ~ dnorm(mu2, tau2);

  mu2 <- 3 + (x1 - 2) / 4;

  tau2 <- 2 / 115;

}

In Stan, this model would be like 

parameters {

  double x1;

  double x2;

}

model {

   x1 ~ normal(2 + (x2 - 3) / 6, sqrt(115 / 3));

   x2 ~ normal(3 + (x1 - 2) / 4, sqrt(115 / 2)); 

}

However, this would be problematic in Stan as Stan does not check if there is any cycle. 

And this get *directly* translated into a logposterior as 

normal(2 + (x2 - 3) / 6, sqrt(115 / 3)) * normal(3 + (x1 - 2) / 4, sqrt(115 / 2)),

which is not what we want.  

Actually, when we did the Wawro's example, I have noticed that and pointed out that

there are cycles in the model specification.  Back then, I did not how 

stan parse the model.  I am not sure that is a problem. 

--

Jiqiang 

[Ben Goodrich]

Ruh-roh ...

On Jan 26, 5:44 pm, Jiqiang Guo <guoj...@gmail.com> wrote:
> # It turns out that the model specified in OpenBUGS
> # has *cycles*, which is not support by Stan (and JAGS).

Just to make sure I understand you correctly, this version of the ice
example

http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html#SECTION000121000000000000000

has a Directed Cycle in the top right of the graph at the bottom,
which Stan cannot parse correctly yet?

If so, does that imply that vanilla_ice should be using this

http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#SECTION000122000000000000000

reparameterization, which as far as I can tell does not have a
Directed Cycle in its graph at the bottom?

Thanks,
Ben

[Jiqiang Guo]

On Thu, Jan 26, 2012 at 6:05 PM, Ben Goodrich <goodri...@gmail.com> wrote:

Ruh-roh ...

On Jan 26, 5:44 pm, Jiqiang Guo <guoj...@gmail.com> wrote:
> # It turns out that the model specified in OpenBUGS
> # has *cycles*, which is not support by Stan (and JAGS).

Just to make sure I understand you correctly, this version of the ice
example

http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html#SECTION000121000000000000000

has a Directed Cycle in the top right of the graph at the bottom,
which Stan cannot parse correctly yet?

This is the cycle I am talking, it is like Gibbs (or conditionally

specified models).  Stan could parse it, but the result is 

not correct.   

If so, does that imply that vanilla_ice should be using this

http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#SECTION000122000000000000000

reparameterization, which as far as I can tell does not have a
Directed Cycle in its graph at the bottom?

I have not figured out why the two forms of ice example are equivalent. 

If the model is conditionally specified correctly, we should have a 

joint distribution (an equivalent way), but the closed form might not 

be easy to obtain.  In other words, I am not sure we can always 

avoid the cycle.

Thanks,
Ben

[Ben Goodrich]

I'd be surprised if we could always reparameterize a directed cycle to
yield an undirected cycle, but are we in agreement that the alternate
parameterization is what I should be utilizing in Stan for Wawro's
models this week?

On Jan 26, 6:20 pm, Jiqiang Guo <guoj...@gmail.com> wrote:

> On Thu, Jan 26, 2012 at 6:05 PM, Ben Goodrich <goodrich....@gmail.com>wrote:
>
> > Ruh-roh ...
>
> > On Jan 26, 5:44 pm, Jiqiang Guo <guoj...@gmail.com> wrote:
> > > # It turns out that the model specified in OpenBUGS
> > > # has *cycles*, which is not support by Stan (and JAGS).
>
> > Just to make sure I understand you correctly, this version of the ice
> > example
>

> >http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html...

>
> > has a Directed Cycle in the top right of the graph at the bottom,
> > which Stan cannot parse correctly yet?
>
> This is the cycle I am talking, it is like Gibbs (or conditionally
> specified models).  Stan could parse it, but the result is
> not correct.
>
>
>
> > If so, does that imply that vanilla_ice should be using this
>

> >http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html...

[Jiqiang Guo]

Yes. 

[Ben Goodrich]

Good to know. Thanks. -- Ben

[Ben Goodrich]

In revision 598 I have committed pineapple_ice.stan, which attempts to
replicate the (thankfully, somewhat simpler) parameterization of the
ice model at

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#SECTION000122000000000000000

This parameterization does not have a directed cycle, but it does have
a redundancy where a positive scale parameter (sigma) multiplies a
vector of random intercepts.

If I specify

double(0,) sigma;

stan does not run very well. If it is lucky, it will escape the burnin
period but seems to not converge and yields odd results.

If I cheat a bit and specify

double(0.05,) sigma;

stan runs a lot better and the posterior mean for sigma is a little
more than 0.05, which is in line with the posterior mean of sigma on
page 31 of

http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/Vol2.pdf

However, our posterior log Risk Ratios run from -1.43 to 1.06, whereas
theirs run from -1.08 to 1.06, even when I do half a million
iterations.

Although there is a possibility that WinBUGS didn't converge, we've
done it with 350K iterations.

So, it seems that there are a few more things to iron out with this
example.

Ben

[Jiqiang Guo]

I guess I messed up when I talked about stuff related to 

graphs.  Sorry for that.  I had only agreed that we can only

parameterize the model in terms of directed *acyclic* graph (DAG). 

For undirected graph (in other words, bi-directed?), so there 

are cycles in terms of DAG.  Hence, I do not think Stan could handle 

undirected graph correctly. Then the pineapple_ice is problematic.  

It seems that I said there is a cycle in the graph,

http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html

I think I am wrong.  There is not.  

And the original version coded in Stan is from classical bugs examples (http://goo.gl/dJlk2)

Be careful when I am talking about graph, I am not confident about that now. 

But if you cannot do the model in JAGS and JAGS tells you that is because of cycles.

If that is the real reason, I do not think Stan could handle it correctly. 

--

Jiqiang 

[Bob Carpenter]

Stan's parser allows cycles. But it's up to you
to make sure the probabilities are correct. For
instance, you can write:

y ~ normal(x,1);
x ~ normal(y,1);

and Stan won't complain. It'll get compiled to

lp__ += normal_log(y,x,1);
lp__ += normal_log(x,y,1);

This is a complete Stan model specifiying joint
probability:

p(x,y) = normal(x|y,1) * normal(y|x,1)

Probably not what we want unless we're hacking an
MI model directly into Stan.

The problem BUGS/JAGS have with cycles conceptually is
that they treat multivariate variable nodes as nodes
of the graph rather than the elements of the vectors.
So if they had:

x[1] ~ normal(0,1);
for (n in 2:N)
x[n] ~ normal(0.5 * x[n-1],1);

They'd think of that as a cycle from x to x because
the variable x shows up on the left-hand and right-hand
sides of the looped sampling.

But the actual random variables are x[1], x[2], ...
and there's no cycle there.

Writing the above kind of thing is totally kosher in
Stan and leads to the right log probability.

- Bob

> <http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#SECTION000122000000000000000>

[Bob Carpenter]

This all sounds like fodder for the convergence
monitoring mill. And some tweaks to adaptation.

Matt -- could you take a look at the model and
see if you see anything that looks fishy?

On 1/27/12 12:36 AM, Ben Goodrich wrote:
> In revision 598 I have committed pineapple_ice.stan, which attempts to
> replicate the (thankfully, somewhat simpler) parameterization of the
> ice model at
>
> www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#SECTION000122000000000000000
>
> This parameterization does not have a directed cycle, but it does have
> a redundancy where a positive scale parameter (sigma) multiplies a
> vector of random intercepts.
>
> If I specify
>
> double(0,) sigma;
>
> stan does not run very well. If it is lucky, it will escape the burnin
> period

Otherwise, it'll just hang during burnin? We want
to avoid this. At the very least, we should be dumping
out things like the step size that adaptation finds.
And we should make sure we don't hang during adaptation.

> but seems to not converge and yields odd results.

Did you apply any kind of convergence monitoring? We don't
have nice Coda-friendly output yet.

> If I cheat a bit and specify
>
> double(0.05,) sigma;
>
> stan runs a lot better and the posterior mean for sigma is a little
> more than 0.05, which is in line with the posterior mean of sigma on
> page 31 of
>
> http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/Vol2.pdf

It'd be more encouraging if you could specify

double(0.01,) sigma;

and still get a value of 0.05.

> However, our posterior log Risk Ratios run from -1.43 to 1.06, whereas
> theirs run from -1.08 to 1.06, even when I do half a million
> iterations.

Did you initialize both in the same place?

> Although there is a possibility that WinBUGS didn't converge, we've
> done it with 350K iterations.

I'd like to see traceplots for the variable in question.
I don't have a good way to run BUGS any more because I
don't have a Windows machine.

> So, it seems that there are a few more things to iron out with this
> example.

That is why we're building so many functional test cases.

- Bob

[Ben Goodrich]

OK, I've confused myself here and quite possibly other people too. So
let me see if I have this correct now.

According to the BUGSExamples R package, this model

https://r-forge.r-project.org/scm/viewvc.php/pkg/inst/examples/Icemodel.txt?view=markup&revision=2&root=bugs-examples

can be estimated in WinBUGS and OpenBUGS but not JAGS. That is
apparently due to the fact that it has what JAGS calls a "cycle".
However, this "cycle" is not actually a cycle in the sense of graph
theory because the graph at the bottom of

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html

is acyclic in the sense that you cannot get back to where you started
from anywhere in the mess in the top right corner. It is only
classified as a "cycle" by JAGS because the parameters are all
contained in the same vector named beta and JAGS does not distinguish
among elements of beta.

So, Stan should be able to do that model, correct?
But, Stan cannot do models with an undirected graph? Like

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#iceSigfig

If so, the question is "Why is vanilla_ice yielding estimates that are
slightly different from those in WinBUGS?" rather than "Why is
pineapple_ice yielding estimates that are slightly different from
those WinBUGS gets when estimating the vanilla_ice model?" Although
the question of why pineapple_ice hangs during burnin if sigma is not
constrained to be sufficiently positive is worth investigating (unless
the answer is simply that Stan doesn't do models that correspond to
undirected graphs).

Ben

> > On Fri, Jan 27, 2012 at 12:36 AM, Ben Goodrich <goodrich....@gmail.com
> > <mailto:goodrich....@gmail.com>> wrote:
>
> >     In revision 598 I have committed pineapple_ice.stan, which attempts to
> >     replicate the (thankfully, somewhat simpler) parameterization of the
> >     ice model at
>

> >    www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#SECTIO...
> >     <http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html...>

[Ben Goodrich]

Iff pineapple_ice (with the undirected graph) is a legit model in
Stan, then it really limps along when I specify

double(0.01,) sigma

but it does finish. And the posterior mean is about 0.015. So, this
may just be a bad example in the sense that they don't have much data
and putting a gamma(1,1) prior on sigma puts more mass at zero than
anywhere else.

I haven't been able to find out what BUGS used as the inits for the
pineapple_ice model. (we do know what they are for vanilla_ice).

Ben

On Jan 27, 12:03 pm, Ben Goodrich <goodrich....@gmail.com> wrote:
> OK, I've confused myself here and quite possibly other people too. So
> let me see if I have this correct now.
>
> According to the BUGSExamples R package, this model
>

> https://r-forge.r-project.org/scm/viewvc.php/pkg/inst/examples/Icemod...

[Jiqiang Guo]

For better or worse, I commented below. 

On Fri, Jan 27, 2012 at 12:03 PM, Ben Goodrich <goodri...@gmail.com> wrote:

OK, I've confused myself here and quite possibly other people too. So
let me see if I have this correct now.

According to the BUGSExamples R package, this model

https://r-forge.r-project.org/scm/viewvc.php/pkg/inst/examples/Icemodel.txt?view=markup&revision=2&root=bugs-examples

The code that the link points to is not the model represented by your following link.  

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html 

can be estimated in WinBUGS and OpenBUGS but not JAGS. That is
apparently due to the fact that it has what JAGS calls a "cycle".
However, this "cycle" is not actually a cycle in the sense of graph
theory because the graph at the bottom of

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html 

There is no cycles in terms of elements of beta vector.  So Stan

and JAGS can handle it correctly. 

is acyclic in the sense that you cannot get back to where you started
from anywhere in the mess in the top right corner. It is only
classified as a "cycle" by JAGS because the parameters are all
contained in the same vector named beta and JAGS does not distinguish
among elements of beta.

JAGS can handle this example,  

http://mcmc-jags.hg.sourceforge.net/hgweb/mcmc-jags/examples/file/1771d5e24e1b/classic-bugs/vol2/ice/icear.bug

So, Stan should be able to do that model, correct?
But, Stan cannot do models with an undirected graph? Like

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#iceSigfig

If so, the question is "Why is vanilla_ice yielding estimates that are
slightly different from those in WinBUGS?" rather than "Why is
pineapple_ice yielding estimates that are slightly different from
those WinBUGS gets when estimating the vanilla_ice model?" Although
the question of why pineapple_ice hangs during burnin if sigma is not
constrained to be sufficiently positive is worth investigating (unless
the answer is simply that Stan doesn't do models that correspond to
undirected graphs).

I think both vanilla and pineapple have cycles, so they are not correctly

handled in Stan as what we might want in BUGS. 

[Ben Goodrich]

On Jan 27, 12:31 pm, Jiqiang Guo <guoj...@gmail.com> wrote:
> > According to the BUGSExamples R package, this model
>

> >https://r-forge.r-project.org/scm/viewvc.php/pkg/inst/examples/Icemod...

>
> The code that the link points to is not the model represented by your
> following link.
>
> www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html

Right, the current WinBUGS code differs from the classic BUGS code,
but I think the same graph corresponds to both. At least, the text is
essentially the same describing the motivation for the model in both
the WinBUGS and BUGS manuals.

> JAGS can handle this example
>
> http://mcmc-jags.hg.sourceforge.net/hgweb/mcmc-jags/examples/file/1771d5e24e1b/classic-bugs/vol2/ice/icear.bug

Yes, but it isn't clear to me how the model that JAGS can handle
relates to the WinBUGS model it is trying to replicate. In any event,
the differences in the posterior estimates relative to the WinBUGS
estimates seem to be about the same magnitude as the differences Stan
gets whether for ice, vanilla_ice, or pineapple_ice.

> I think both vanilla and pineapple have cycles, so they are not correctly
> handled in Stan as what we might want in BUGS.

Hmm. I would have said that neither has a cycle in the sense of

http://en.wikipedia.org/wiki/Directed_acyclic_graph

but I can see (now that Bob explained it) how JAGS would call it a
"cycle" in the sense that one element affects another element within
the same beta vector.

Ben

[Bob Carpenter]

On 1/27/12 12:03 PM, Ben Goodrich wrote:
...

> But, Stan cannot do models with an undirected graph? Like
>
> www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node14.html#iceSigfig

Stan should be able to handle anything for which you
can:

a. code up the log probability function, and

b. transform all variables to be unconstrained.

Although we don't support undirected graphical models
directly, you can just write out their log probabilities
directly in Stan.

That's what BUGS did in its model. You just need
to be careful you get a proper density (up to a constant).

There are conceptual problems with other BUGS models, like
ordinal logit. To get ordered points c[1] and
c[2], they write something tha also looks circular:

c[1] ~ normal(0,1) I(,c[2]);
c[2] ~ normal(0,1) I(c[1],);

This implies an ill-formed density

p(c) = normal(c[1] | 0, 1) / (1 - pnorm(c[2]|0,1))
* normal(c[2] | 0, 1) / pnorm(c[1]|0,1)

but is OK for sampling using Gibbs as the
conditionals are correct.

What you really need in the model spec
is something like the way you'd write the "same"
thing in Stan:

c[1] ~ normal(0,1)
c[2] ~ normal(0,1) I(c[1],)

This provides a proper density.

BUGS can't handle this version, because there's a chance
you'll sample c[1] from a point above c[2] without
the constraint.

- Bob

[Bob Carpenter]

On 1/27/12 12:51 PM, Ben Goodrich wrote:

> Yes, but it isn't clear to me how the model that JAGS can handle
> relates to the WinBUGS model it is trying to replicate. In any event,
> the differences in the posterior estimates relative to the WinBUGS
> estimates seem to be about the same magnitude as the differences Stan
> gets whether for ice, vanilla_ice, or pineapple_ice.

That's encouraging, at least.

>> I think both vanilla and pineapple have cycles, so they are not correctly
>> handled in Stan as what we might want in BUGS.
>
> Hmm. I would have said that neither has a cycle in the sense of
>
> http://en.wikipedia.org/wiki/Directed_acyclic_graph
>
> but I can see (now that Bob explained it) how JAGS would call it a
> "cycle" in the sense that one element affects another element within
> the same beta vector.

Right -- these are no problem for Stan (again, assuming
the joint density defined makes sense).

- Bob

[Matt Hoffman]

I think this is not a correct way of specifying this model. I guess it somehow works in BUGS because they're only ever considering one variable at a time. (Although it seems like you have to exploit some really weird "features" of BUGS to get it to work.) But in general it's by no means true that the joint probability of a bunch of variables is proportional to the product of the conditional probabilities of those variables. For example, the conditional probability of β_1 should be

p(β_1|β_-1) = N(β_1; 2*β_2 - β_3, 1)

but in the model as specified in Stan I think it'll be something like

p(β_1|β_-1) = N(β_1; 2*β_2 - β_3 + 2/5*β_2 - 1/6*β_3, 1/(1+1/5+1/6)).

The issue is that there are a bunch of terms that are being double-counted.

So maybe it's not so surprising that it's giving weird results?

Matt

[Jiqiang Guo]

I am trying to find the corresponding model in BUGS for pineapple and vanilla later. 

But ice.stan and JAGS'  

http://mcmc-jags.hg.sourceforge.net/hgweb/mcmc-jags/examples/file/1771d5e24e1b/classic-bugs/vol2/ice/icear.bug

is trying to replicate this example in bugs:

http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node13.html

> I think both vanilla and pineapple have cycles, so they are not correctly
> handled in Stan as what we might want in BUGS.
Hmm. I would have said that neither has a cycle in the sense of
http://en.wikipedia.org/wiki/Directed_acyclic_graph
but I can see (now that Bob explained it) how JAGS would call it a
"cycle" in the sense that one element affects another element within
the same beta vector.

 

I think it has cycles because I take (not any definition from somewhere) that

undirected are bi-directed.  In other words, if I have say

x | y ~ normal(,), and y | x ~ normal(,), I would graph it as 

x --------- y or 

x---------> y    (x <------> y) 

 \<--------/

From the second graph, I think it is a cycle. 

By the way, I take Bob's explanation with a gain of salt.  Because 

I can use JAGS (which forbid cycles) to sample from the following model:

model {

  a[1] ~ dnorm(0, 1);

  a[2] ~ dnorm(a[1], 1)

}

Without talking about graphs, I dot not think the log-posterior in pineapple and 

vanilla specified by Stan would lead to the log-posterior we are really interested. 

The fact they have close results is just a coincidence, I guess. 

I would quote from p3057 of Lunn et al (2009)

http://www.stat.columbia.edu/~gelman/bayescomputation/lunnbugswithcomments.pdf

on the undirected graph since I could not find other reference: 

The majority of such distributions actually
correspond to undirected sub-graphs. However, they are easily accommodated within the directed
graphical framework of BUGS by simply thinking of them as multivariate quantities represented
by a single node in the graph with complex internal structure—if all internal nodes can be updated
together, then the internal structure is unimportant for the purposes of performing the encompassing
Gibbs scheme

-- 

Jiqiang 

[Bob Carpenter]

On 1/27/12 1:48 PM, Jiqiang Guo wrote:
>
> > I think both vanilla and pineapple have cycles, so they are not
> correctly
>

> I think it has cycles because I take (not any definition from somewhere)
> that
> undirected are bi-directed.

That's not how cycles are usually defined in undirected
graphs. If it were, every undirected graph with at least
one edge would have a cycle. Usually it means that
there is a non-empty, non-repeating path of edges from a
node back to itself.

> In other words, if I have say
> x | y ~ normal(,), and y | x ~ normal(,), I would graph it as
>
> x --------- y or
>
> x---------> y (x <------> y)
> \<--------/
>
> From the second graph, I think it is a cycle.

I think I'm seeing the confusion here.

There's a whole literature on "undirected graphical
models", which are used for Markov random fields.
(See, e.g., Koller's book "Probabilisitc Graphical
Models"), or the Wikipedia entry:

http://en.wikipedia.org/wiki/Markov_random_field

I now realize that BUGS is using the terminology
from both graph theory and graphical models in
exactly the non-standard way Jiqiang thought.
That is, they're using

{ a --- b }

as shorthand for

{ a --> b, a <-- b }

> By the way, I take Bob's explanation with a gain of salt. Because
> I can use JAGS (which forbid cycles) to sample from the following model:
>
> model {
> a[1] ~ dnorm(0, 1);
> a[2] ~ dnorm(a[1], 1)
> }

I would, too, if JAGS allows that model. What kind of cycle
does JAGS NOT allow?

> Without talking about graphs, I dot not think the log-posterior in
> pineapple and
> vanilla specified by Stan would lead to the log-posterior we are really
> interested.
> The fact they have close results is just a coincidence, I guess.

Is there a way to write down the log posterior you
are interested in up to an additive constant?
If so, it should be possible to just code it directly in
Stan. The Stan statement

a ~ foo(theta);

is just shorthand for:

lp__ += foo_log(a,theta);

You don't have to worry about how it gets computed, as
they seem to be doing in BUGS/JAGS with their conditionals.

Of course, as Matt pointed out, you aren't guaranteed a proper
joint model by writing out a product of conditionals (I'm
sure you guys are used to this discussion from MI!).

> I would quote from p3057 of Lunn et al (2009)
> http://www.stat.columbia.edu/~gelman/bayescomputation/lunnbugswithcomments.pdf
> on the undirected graph since I could not find other reference:
>
> The majority of such distributions actually
> correspond to undirected sub-graphs. However, they are easily
> accommodated within the directed
> graphical framework of BUGS by simply thinking of them as
> multivariate quantities represented

> by a single node in the graph with complex internal structure�if all

> internal nodes can be updated
> together, then the internal structure is unimportant for the
> purposes of performing the encompassing
> Gibbs scheme

I see what they're getting at for sampling. You do the
same thing for computing potential functions in MRFs,
looking at maximal cliques (largest fully connected subgraphs).

Grouping isn't relevant for Stan -- it does all params at
the same time like Metropolis. I'm not sure I understand what
Lunn et al. mean in terms of modeling, but I suspect it's
what Jiqiang said, using an undirected edge as shorthand for
both directed edges.

- Bob

[Ben Goodrich]

On Jan 27, 1:48 pm, Jiqiang Guo <guoj...@gmail.com> wrote:
> I am trying to find the corresponding model in BUGS for pineapple and
> vanilla later.

I said it wrong earlier. This is what I think is true at approximately
2:15 PM (subject to change at 2:16 PM or later):

Classic BUGS lists 3 ways to do the "same" model here:

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node12.html

The corresponding text files are available here:

http://www.stat.uiowa.edu/ftp/woodwort/BUGSXmpl/Classic%20Examples/EXAMPLES/VOL2/ICE/

The first ("Autoregressive smoothing of relative risks") is
implemented in ICEAR.BUG are corresponds to the example in JAGS.

The second ("An undirected model using an intrinsic prior for the
random effects") is implemented in ICECARSG.BUG and corresponds to
pineapple_ice.

The third ("An intrinsic prior with a hyperparameter") is implemented
in ICECAR.BUG and (in theory) corresponds to both the example in
WinBUGS and vanilla_ice.

When I said that the WinBUGS example and BUGS example refer to the
same graph, that was wrong. The WinBUGS example corresponds to this
graph

http://www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node15.html#iceCARfig

which they refer to as "undirected", presumably because there are no
arrows linking the betas (only arrows linking tau to the betas).

They also have a table of posterior estimates from the three versions
of the model. However, the results in the third column differ from
those given in the WinBUGS manual, from the other two columns, and
from the original article. For vanilla_ice, Stan gets close on sigma
but not on the logRRs.

Hopefully that addresses a part of the confusion.

Ben

[Matt Hoffman]

Note that my point wasn't that you aren't guaranteed a joint model by writing the product of conditionals, it was that taking a product of conditionals does not give you the joint to which those conditionals correspond. 

I'm not saying that the model in pineapple_ice.stan isn't a valid model, I'm saying that it's not equivalent to the model in ice.stan, or to the model as it's interpreted by BUGS, and so we shouldn't expect it to produce the same results. Because it's a different model.

Matt

   by a single node in the graph with complex internal structure—if all

   internal nodes can be updated
   together, then the internal structure is unimportant for the
   purposes of performing the encompassing
   Gibbs scheme

[Ben Goodrich]

Right. The pineapple_ice model is a bit of a distraction (that I
caused) at this point because that is neither the model we are trying
to estimate, nor the model that WinBUGS does estimate, nor the model
that BUGS recommends estimating. For the model vanilla_ice is trying
to replicate:

www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node15.html

speaks to Matt's point somewhat. They work out the full-conditional
distribution of tau and stick into the BUGS model. Apparently, BUGS
has some precedence heuristics where if you do that it adjusts its
behavior accordingly. I am guessing that part of the problem we are
getting with vanilla_ice in Stan is that we shouldn't just be
literally translating that BUGS code into Stan code because the bit
about the FCD of tau is out of place in Stan. If I am correct, then we
should just put a vague prior on tau and roll?

Ben

> >  http://en.wikipedia.org/wiki/**Markov_random_field<http://en.wikipedia.org/wiki/Markov_random_field>

> >>http://www.stat.columbia.edu/~**gelman/bayescomputation/**
> >> lunnbugswithcomments.pdf<http://www.stat.columbia.edu/~gelman/bayescomputation/lunnbugswithcom...>

[Bob Carpenter]

On 1/27/12 2:29 PM, Matt Hoffman wrote:
> Note that my point wasn't that you aren't guaranteed a joint model by
> writing the product of conditionals, it was that taking a product of
> conditionals does not give you the joint to which those conditionals
> correspond.

Thanks for clarifying -- there are two issues here.

One is that we want a properly specified model,
and multiplying a bunch of conditionals might
not do that.

Two is what Matt just said, and what I was
trying to get at with the ordinal logit cutpoint
example. The only cases I know where the product
of conditionals give you the joint is when there
is independence. So directly using the BUGS formulation
can lead you astray as they can get away with conditionals
specified up to a proportion.

Bottom line (literally): for Stan, the only thing
that matters is getting the correct log posterior (or log joint)
up to a proportion.

- Bob

[Matt Hoffman]

To reemphasize what Bob just said, the issue isn't just a question of vanilla versus pineapple. It's that there is no way in Stan to specify a distribution in terms of its conditionals. None. You have to specify the full joint via some kind of factorization that respects the chain rule of probability. So it's not a subtle question of the way that pineapple_ice or vanilla_ice handles tau, it's that there simply is no correct way to specify the ice model this way in Stan.

[Bob Carpenter]

On 1/27/12 2:37 PM, Ben Goodrich wrote:
> Right. The pineapple_ice model is a bit of a distraction (that I
> caused) at this point because that is neither the model we are trying
> to estimate, nor the model that WinBUGS does estimate, nor the model
> that BUGS recommends estimating. For the model vanilla_ice is trying
> to replicate:
>
> www.mrc-bsu.cam.ac.uk/bugs/documentation/exampVol2/node15.html
>
> speaks to Matt's point somewhat. They work out the full-conditional
> distribution of tau and stick into the BUGS model. Apparently, BUGS
> has some precedence heuristics where if you do that it adjusts its
> behavior accordingly. I am guessing that part of the problem we are
> getting with vanilla_ice in Stan is that we shouldn't just be
> literally translating that BUGS code into Stan code because the bit
> about the FCD of tau is out of place in Stan. If I am correct, then we
> should just put a vague prior on tau and roll?

You could certainly see if that gave you a similar
result. Or just constrain it to the value they get
and see if everything else works before trying to
estimate tau.

Or you can just write down their joint probability directly.
They give part of the derivation, but I still don't really
see what their model's trying to do.

- Bob

[Ben Goodrich]

Ergh. The part of my thinking that I neglected to actually write down
was that we should write in the .stan file a joint prior on the betas
that corresponds to that implied by WinBUGS's Conditional
AutoRegressive prior. Then put a vague prior on tau, then those priors
and the likelihood should yield a posterior that Stan can sample from.
I think (if I am reading the original article correctly) that the
joint prior on the betas is just a multivariate normal with a very
particular covariance matrix (even more particular than the usual AR1
covariance matrix). So, this seems to be more of an issue of (Win)BUGS
making things more complicated than they need to be because (Win)BUGS
is not that good with joint priors over multivariate parameters?

On Jan 27, 2:50 pm, Matt Hoffman <mdhof...@cs.princeton.edu> wrote:
> To reemphasize what Bob just said, the issue isn't just a question of
> vanilla versus pineapple. It's that there is no way in Stan to specify a
> distribution in terms of its conditionals. None. You have to specify the
> full joint via some kind of factorization that respects the chain rule of
> probability. So it's not a subtle question of the way that pineapple_ice or
> vanilla_ice handles tau, it's that there simply is no correct way to
> specify the ice model this way in Stan.
>

[Jiqiang Guo]

On Fri, Jan 27, 2012 at 2:15 PM, Bob Carpenter <ca...@alias-i.com> wrote:

By the way, I[Jiqiang] take Bob's explanation with a grain of salt.  Because

I can use JAGS (which forbid cycles) to sample from the following model:

model {
  a[1] ~ dnorm(0, 1);
  a[2] ~ dnorm(a[1], 1)
}

I would, too, if JAGS allows that model.  What kind of cycle
does JAGS NOT allow?

This is an example supported by OpenBUGS but NOT JAGS.  So 

what BUGS really supports is beyond what Lunn et al says in the 

paper, but I could not find the document. 

## (x1, x2)' ~ mutil_normal with 

## mean (0, 0) 

## var  ( 100, 24 )

##       ( 24,    9  )

model {

    x1 ~ dnorm(mu1, tau1);

    x2 ~ dnorm(mu2, tau2);

    mu1 <- 8 * x2 / 3;

    mu2 <- .24 * x1;

    tau1 <- 1 / 36;

    tau2 <- 1 / 3.24;

}

 

[Andrew Gelman]

Just to be clear: the example below _wouldn't_ actually hack an MI model. Rather, it's a "pseudo-likelihood" (which I think might be the same as "naive Bayes") model in which, instead of working with p(x,y), you work with p(x|y)*p(y|x), for computational convenience. We might be able to come up with a setting where it would make sense, but MI isn't it.
A

[Andrew Gelman]

We should remember to discuss some of these struggles in the article and manual.

[Bob Carpenter]

Thanks for the clarification. I'll just
cut out the analogies to MI, which I clearly
still don't understand.

I was erroneously assuming that if you had
predictors x1, x2 and x3, and x1 and x2
were missing, that you'd sample x2|x1,x3
and x1|x2,x3. Do you stage them instead
and do something like x2|x3 and x1|x2,x3
or just separate them and do x2|x3 and x1|x3?

In natural language processing and machine
learning classifiers, the term "naive" is used
for making (unwarranted) independence assumptions,
such as in this case assuming

p(x|y) = p(x), and

p(y|x) = p(y).

- Bob

[Andrew Gelman]

Bob:
Yes, in MI you would sample x|y and sample y|x from specified (if incoherent) joint densities. But that's not the same as what Stan does which is to multiply these two densities together.
A

[Andrew Gelman]

I'm also wondering if the model will converge under various improvements. Let me label some things:

Nuts 1: The original Nuts algorithm
Nuts 2: The improvement that Matt is planning to put into Stan 1.0, which has an improved warm-up procedure and diagonal adaptation of the mass matrix
Nuts 3: The new Nuts that Matt was talking about that uses 3rd derivative info, takes about 3x as long as Nuts 1 or 2, and works on the triangle, log-gamma, and Cauchy examples (all of which defeat Nuts 1 and 2)
Nuts 4: The newer Nuts that Matt was discussing with Mike B., which does some sort of Riemannian something or another.

Now consider this slow model we've all been emailing about. Stan now uses Nuts 1. One thing I'm wondering about is, would the model still run slow in Nuts 2? What about Nuts 3?

[Bob Carpenter]

Andrew -- you're several e-mails behind.
The conclusion was that it was a specification
problem, not a convergence problem.

All good questions, though.

Yes, all of this is going in the manual, in a
simplified, "here's what you should do" (and
maybe "here's what'll happen if you don't" form).

- Bob

[Ben Goodrich]

On Jan 27, 6:09 pm, Bob Carpenter <c...@alias-i.com> wrote:
> Andrew -- you're several e-mails behind.
> The conclusion was that it was a specification
> problem, not a convergence problem.

I am now pretty sure that the joint prior for the betas is just
multivariate normal with mean vector zero and a known, very weird,
precision matrix. So, I would need to implement that, which we had
talked about before a little bit.

Do you want a multi_normal_log() whose signature includes a matrix of
doubles rather than a templated matrix or do you want a function with
a different name like multi_normal_log_known_precision()? Also, in
this case, the precision matrix happens to be singular (thus it's an
improper prior) so it *can't* allow propto be true; the
log(det(Sigma)) part would be -Inf. What do you want to do about that?
I would lean toward making a multi_normal_log_improper() function that
ignores the constant always.

Ben

[Andrew Gelman]

Yes, in Stan it should definitely be possible to specify a joint multivariate distribution (such as a conditional autoregression) by just writing down the log-conditional multivariate density.
A

[Ben Goodrich]

Noooooooooooooooooooooooooooooo!

It took us 38 emails to get to the point where we agreed that
(Win)BUGS was confusing everything by mangling a multivariate normal
prior with known mean and precision matrix into a set of conditionals
so that it can do Gibbs, whereas Stan (or M-H for that matter) can
keep everything simple.

[Andrew Gelman]

B
I don't see why you say No (let alone Nooooooooo). I think we're saying the same thing: In San, we can specifiy a log-joint density for any multivariate distribution with known normalizing constant; we _don't_ need to specify it in terms of conditionals. For example, we can specify a condtional autoregression by mathematically working out its log-joint density and just programming it directly.
A

[Ben Goodrich]

In that case, Noooooooooooooooooooooooooooooo! means yes. I though
that by

"Yes, in Stan it should definitely be possible to specify a joint
multivariate distribution (such as a conditional autoregression) by
just writing down the log-conditional multivariate density."

you were talking about Stan's parser inferring the joint from
conditionals in the .stan file, which would be very, very hard.

In general, there are a lot of tricks to get (Win)BUGS to do something
interesting (LISREL is the same way) that should be a lot more
straightforward in Stan, in theory. But if people are coming from the
(Win)BUGS world and trying to translate their tricky (Win)BUGS code
(or Wawro's) into Stan, then they will be wrong and / or confused.

[Andrew Gelman]

Yeeeeeeeessssssssssss, I agree. It's a key benefit of Stan that we can program the model directly and forget the damn tricks.

[Bob Carpenter]

It's now our duty to spare others two score
or more e-mails on how to translate a BUGS model
into Stan. So to repeat what Andrew keeps saying,
we should doc all this. I think just the BUGS
model translations themselves will go a long way
for users who are example driven.

- Bob

[Bob Carpenter]

On 1/27/12 8:20 PM, Ben Goodrich wrote:

> I am now pretty sure that the joint prior for the betas is just
> multivariate normal with mean vector zero and a known, very weird,
> precision matrix. So, I would need to implement that, which we had
> talked about before a little bit.

> Do you want a multi_normal_log() whose signature includes a matrix of
> doubles rather than a templated matrix or do you want a function with
> a different name like multi_normal_log_known_precision()?

This is going to lead to an unwieldy profusion of
densities and functions if we're not careful.

We're eventually going to need a multi_normal_log
with a precision matrix instead of covariance matrix
parameter and then there's the full-matrix vs.
Cholesky factor split.

> Also, in
> this case, the precision matrix happens to be singular (thus it's an
> improper prior) so it *can't* allow propto be true;the
> log(det(Sigma)) part would be -Inf. What do you want to do about that?
> I would lean toward making a multi_normal_log_improper() function that
> ignores the constant always.

I think you mean we can't allow propto to be false.
propto=false means include every term in the log
probability function, whereas propto=true includes
only terms containing at least one stan::agrad::var.

The precision matrix will be involved in two terms,
the log(det(Tau)), which we wouldn't evaluate if Tau
was made up of doubles, and the product term with location (mu)
and outcome (y) vectors, which gets evaluated unless the outcome
and both parameters are constant.

What seems to be needed is an "improper" flag
that turns off the check for singularity. But then
there's still the issue of how to pass flags from
Stan to configure the error-checking on a
call-by-call basis.

One way to go would be as Ben suggests, splitting out
the functions, and using names like:

multi_normal_precision_improper_log()

The "_log" needs to go last so that you can refer to
the prefix in stan, as in:

y ~ multi_normal_precision_improper(mu,Tau);

On the other hand, maybe we could add constants to Stan
representing error-checking policies, so that you could
write

y ~ multi_normal_precision(mu,Tau,IMPROPER);

and the IMPROPER would get translated to an error policy
in C++ and passed to the multi_normal_precision function.

But this is problematic in that we only want to turn off
checks for the precision matrix -- we'd still check that
none of the mu or y were NaN and the mu were finite.

Speaking of this latter decision, Matt argued that
outcome y and location mu should be treated symmetrically,
as they're symmetric terms in the probability function.

Boost throws exceptions if mu is infinite when it
constructs a normal_distro(mu,sigma) (Boost doesn't do
multivariate). I thought this made sense, because
a normal_distro(mu,sigma) will always returns 0
probability.

On the other hand, if you have a finite mu and positive
finite sigma, then normal_distro(y|mu,sigma) returns
0. So I also see where Boost's coming from.

The question is, what should we do?

- Bob

[Ben Goodrich]

Yes, I meant propto has to (at least in effect) be true when the
precision matrix is singular.

Another possibility for the case where the precision matrix is known
is to precalculate the log determinant, pass it to the function, and
have the function check to see whether it is Inf. In other words, even
if the precision matrix is invertible and the user wants to normalize,
we don't want the function to be recalculating the determinant every
time the function is called anyway.

Relatedly, if the precision matrix is known, I think we would want to
precalculate its 3-factor Cholesky decomposition, Tau = LDL', set U :=
D^(1/2) * L', and pass U rather than Tau. This would save flops inside
the density function, which would basically just be dotself(U.template
triangularView<Eigen::Upper>() * (y - mu)).

So, I don't think we need a version that takes the whole Tau, but the
rest of Bob's questions remain. Also, maybe now is a good time to
reconsider whether we need a CovMatrix class, which could have
attributes like CheckMe, IsPSD, IsPD, IsInverted, LogDeterminant,
CholeskyFactorization, etc. That way the multi_normal_log() functions
could rely more on polymorphing and have the inside of the function
figure out what calculations it needs to do and what calculations it
can avoid by looking at the attributes.

Ben

[Bob Carpenter]

There are two benefits to allowing:

y ~ multi_normal_precision(mu,Tau,UNCHECKED);

1. Allowing improper precision matrix priors
(and presumably other improper calculations).

2. Letting us have a good way to have one
distribution call another without double
error checking. For instance, with infinite
degrees of freedom, we want the Student-t
to call the normal. Daniel also pointed out
that the errors need to be checked by the Student-t
so that the error messages make sense to the user.

- Bob

[Bob Carpenter]

On 1/28/12 4:17 PM, Ben Goodrich wrote:
> Yes, I meant propto has to (at least in effect) be true when the
> precision matrix is singular.
>
> Another possibility for the case where the precision matrix is known
> is to precalculate the log determinant, pass it to the function, and
> have the function check to see whether it is Inf. In other words, even
> if the precision matrix is invertible and the user wants to normalize,
> we don't want the function to be recalculating the determinant every
> time the function is called anyway.

I've thought about that. We can vectorize the
outcomes and save recalculating inverses and
determinants. But that still does them every
leapfrog step, and it'd be nice to have them
in data. Something like:

data {
cov_matrix(K) Sigma;
}
transformed data {
cov_matrix(K) Tau;
double log_det_Tau;

Tau <- inverse(Sigma);
log_det_Tau <- log(determinant(Tau));
}
model {
for (n in 1:N)
y[n] ~ multi_normal_prec_det(mu,Tau,log_det_Tau);
}

This'll open the possible error that the log_det_Tau
argument isn't actually log(det(Tau)).

We could also use this exact same strategy
if Sigma were a parameter -- just replace
"data" with "parameters" and "transformed data"
with "transformed parameters". Better yet, use
a local, as in:

parameters {
cov_matrix(K) Sigma;
}
model {
cov_matrix(K) Tau;
double log_det_Tau;

Tau <- inverse(Sigma);
log_det_Tau <- log(determinant(Tau));
for (n in 1:N)
y[n] ~ multi_normal_prec_det(mu,Tau,log_det_Tau);
}

because you probably don't want Tau and log_det_Tau
output as parameters as they would be if declared
as "transformed parameters".

Another minor point is that you'll have to code up
the improper prior with a matrix(K,K) type instead
of a cov_matrix(K) type, or you'll get an error reading
in the data.

> Relatedly, if the precision matrix is known, I think we would want to
> precalculate its 3-factor Cholesky decomposition, Tau = LDL', set U :=
> D^(1/2) * L', and pass U rather than Tau. This would save flops inside
> the density function, which would basically just be dotself(U.template
> triangularView<Eigen::Upper>() * (y - mu)).

Exactly. And it'd get coded the exact same way.

> So, I don't think we need a version that takes the whole Tau, but the
> rest of Bob's questions remain. Also, maybe now is a good time to
> reconsider whether we need a CovMatrix class, which could have
> attributes like CheckMe, IsPSD, IsPD, IsInverted, LogDeterminant,
> CholeskyFactorization, etc. That way the multi_normal_log() functions
> could rely more on polymorphing and have the inside of the function
> figure out what calculations it needs to do and what calculations it
> can avoid by looking at the attributes.

I'll have to think more about this.

I've been designing everything with a much more
functional than object orientation. It's much
simpler conceptually as there's no pre-conditions
or post-conditions in an object's state to maintain.

So I've been thinking much more of this kind of
optimization work would be done by rewriting the grammars,
which I have to admit, is far from ideal for readability.

I don't see how to maintain these states given
the flexibility of Eigen matrice's operator[],
which returns a REFERNENCE to an element:

Matrix<T,...>

T& operator[](size_t n) { ... }

If someone gets such a reference and changes it,
we can't guarantee any of our checks like isInverted
or isPosDef or whatever still hold.

So we'd have to change the way things are accessed
for matrices and disallow assignments like:

Sigma[m,n] <- S[m,n] * s[m] * s[n];

Or let the user explicitly apply tests, but that has
the same problem of uglifying the models.

- Bob

[Ben Goodrich]

I am sure there are use cases for a parameterized Tau. In fact, I
think Andrew sent me a paper about a distribution for precision
matrices. But I in the case below:

> We could also use this exact same strategy
> if Sigma were a parameter -- just replace
> "data" with "parameters" and "transformed data"
> with "transformed parameters".  Better yet, use
> a local, as in:
>
>    parameters {
>        cov_matrix(K) Sigma;
>    }
>    model {
>        cov_matrix(K) Tau;
>        double log_det_Tau;
>
>        Tau <- inverse(Sigma);
>        log_det_Tau <- log(determinant(Tau));
>        for (n in 1:N)
>            y[n] ~ multi_normal_prec_det(mu,Tau,log_det_Tau);
>    }

Why wouldn't the user just pass Sigma to multi_normal()? This way is
slower and risks them messing up log_det_Tau somehow.

> > So, I don't think we need a version that takes the whole Tau, but the
> > rest of Bob's questions remain. Also, maybe now is a good time to
> > reconsider whether we need a  CovMatrix class, which could have
> > attributes like CheckMe, IsPSD, IsPD, IsInverted, LogDeterminant,
> > CholeskyFactorization, etc. That way the multi_normal_log() functions
> > could rely more on polymorphing and have the inside of the function
> > figure out what calculations it needs to do and what calculations it
> > can avoid by looking at the attributes.
>
> I'll have to think more about this.
>
> I've been designing everything with a much more
> functional than object orientation.  It's much
> simpler conceptually as there's no pre-conditions
> or post-conditions in an object's state to maintain.

I don't mind functional. I even tolerate R, which is dysfunctional.
But I do kind of like the way the Matrix package in R does object-
orientation. And it is so difficult to get a covariance matrix right
that maybe it should be treated differently.

> I don't see how to maintain these states given
> the flexibility of Eigen matrice's operator[],
> which returns a REFERNENCE to an element:
>
>     Matrix<T,...>
>
>     T& operator[](size_t n) { ... }
>
> If someone gets such a reference and changes it,
> we can't guarantee any of our checks like isInverted
> or isPosDef or whatever still hold.

In the CovMatrix class in my head, there would be constructors that
would take parameters and output a CovMatrix with the various flags
set. But I can't think of a great reason why some rogue function
should subsequently be changing its elements before the posterior is
evaluated. I might even overload [] and () to throw errors. But if
there were some reason why a function should be changing its elements,
then I would say it is that function's responsibility to set CheckMe
to true so that downstream functions know to go through all the
checks.

Maybe I haven't emphasized the following enough:

-- There are a small number of ways to construct a covariance matrix
that guarantee it will be positive definite
-- And there are a great many ways of constructing a covariance matrix
that do not have a snowball's chance in hell of being positive
definite in all 1000 (or even 3) iterations. The subset of positive
definite matrices within the set of symmetric matrices with positive
diagonal elements is just too small.

So, although we probably can't absolutely prevent the following from
happening:

parameters {
double covariances[K * (K - 1) / 2 + K];
}
transformed parameters {
double Sigma[K,K];
int count;

count <- 1;
for (i in 1:K) {
for (j in i:K) {
Sigma[i,j] <- covariances[count];
Sigma[j,i] <- Sigma[i,j]
count <- count + 1;
}
}

We need to do everything we can to discourage the user from doing that
because Sigma will not be positive definite in the vast majority of
the iterations. Users should be using convenience functions in
their .stan files like makeLKJ(), makeAR1(), etc. that output valid
covariance matrices. And when those functions are used, we need some
way of indicating that the density should not check Sigma.

The downside is that maybe in the future somebody invents some cool
way of constructing a positive definite covariance matrix that a user
wants to use but isn't implemented in Stan. In that case, they would
have to do it themselves within the .stan file. But I would almost
rather them email a feature request.

Ben

[Bob Carpenter]

On 1/29/12 2:13 AM, Ben Goodrich wrote:
> I am sure there are use cases for a parameterized Tau. In fact, I
> think Andrew sent me a paper about a distribution for precision
> matrices. But I in the case below:
>
>> We could also use this exact same strategy
>> if Sigma were a parameter -- just replace
>> "data" with "parameters" and "transformed data"
>> with "transformed parameters". Better yet, use
>> a local, as in:
>>
>> parameters {
>> cov_matrix(K) Sigma;
>> }
>> model {
>> cov_matrix(K) Tau;
>> double log_det_Tau;
>>
>> Tau<- inverse(Sigma);
>> log_det_Tau<- log(determinant(Tau));
>> for (n in 1:N)
>> y[n] ~ multi_normal_prec_det(mu,Tau,log_det_Tau);
>> }
>
> Why wouldn't the user just pass Sigma to multi_normal()? This way is
> slower and risks them messing up log_det_Tau somehow.

With Sigma, you need to do all the solving business. Once
you have the precision matrix Tau and log(det(Tau)) it's
just multiply/add. You'd probably want to Cholesky factor
it, which you can also do in Stan.

In general, any sort of transformation runs the huge risk
that the user won't apply the right Jacobian. Here, you're
going to get the same answer back, so it should be OK.

Are you suggesting we extend Eigen with a CovMatrix class?

The nice thing about this is that we could have constructors
that only took the lower-triangular portion of the CovMatrix
efficiently. That'd save half the auto-dif variables compared
to how it's being done now (which could certainly be improved
to at least exploit the self-adjoint views, but that gets back
to bad interactions with heterogeneous operators in auto-dif
in Eigen).

It shouldn't be that hard to extend the Stan language
to treat cov-matrix as a different type than matrix and
declare the multi_normal densities to require a covar matrix
(or precision matrix). Right now, you can plug any matrix
in and the function checks the arguments are OK if error-handling
is turned on.

> But if
> there were some reason why a function should be changing its elements,
> then I would say it is that function's responsibility to set CheckMe
> to true so that downstream functions know to go through all the
> checks.

This is the kind of thing I don't see at all how to do inside
the Stan modeling language. But then we might not need to.

We could have the multivariate densities call a
a checkMe() method on a covariance matrix before it's used.
Then we keep a cached answer that's cleared whenever something
changes the matrix in an unsound way.

We can handle caching of Cholesky factors, log determinants,
etc., the same way.

> Maybe I haven't emphasized the following enough:
>
> -- There are a small number of ways to construct a covariance matrix
> that guarantee it will be positive definite

I get that. It's why we're testing, so we can let
the user know something's going wrong.

I also see the benefit of the kind of positive-definiteness
safety you're talking about.

In fact, there are lots of other places we can use it,
like transforming parameters. We don't want someone doing:

parameter { double y; }
model { log(y) ~ normal(0,1); }

thinking they're going to get a log normal. They can do this
and then add the (in this case very simple) log Jacobian determinant,
but if they just write the above, it's not going to do what
they might expect.

Another big potential source of error is not
declaring parameters with appropriate transforms.
The degenerate case of covariance matrixes is
the univariate case. There, you need a positive
scale parameter. There's nothing to prevent the
user from writing:

y ~ normal(0,-1);

Now we could write a positive-number class, that only
allowed positive number operations on it (square,
square root, log if the value's greater than 1, multiplying
by another positive number). Then require scale parameters
to belong to this class.

Maybe I've been inspired by C++ itself, but I've
given users a nice length rope with which to hang
themselves.

> -- And there are a great many ways of constructing a covariance matrix
> that do not have a snowball's chance in hell of being positive
> definite in all 1000 (or even 3) iterations. The subset of positive
> definite matrices within the set of symmetric matrices with positive
> diagonal elements is just too small.

Right. This is going to be a problem if someone
doesn't declare a matrix parameter as a covar_matrix.

If they initialize to a covariance matrix (we don't
like the snowball's chances on random init), then very
small movements should maintain positive definiteness,
but this is definitely not something we should encourage.

> So, although we probably can't absolutely prevent the following from
> happening:
>
> parameters {
> double covariances[K * (K - 1) / 2 + K];
> }
> transformed parameters {
> double Sigma[K,K];
> int count;
>
> count<- 1;
> for (i in 1:K) {
> for (j in i:K) {
> Sigma[i,j]<- covariances[count];
> Sigma[j,i]<- Sigma[i,j]
> count<- count + 1;
> }
> }
>
> We need to do everything we can to discourage the user from doing that
> because Sigma will not be positive definite in the vast majority of
> the iterations. Users should be using convenience functions in
> their .stan files like makeLKJ(), makeAR1(), etc. that output valid
> covariance matrices. And when those functions are used, we need some
> way of indicating that the density should not check Sigma.

The main method of discouragement at this point is that the
programs won't work.

Which reminds me that we need to instrument them to somehow
stop and fail under some conditions.

> The downside is that maybe in the future somebody invents some cool
> way of constructing a positive definite covariance matrix that a user
> wants to use but isn't implemented in Stan. In that case, they would
> have to do it themselves within the .stan file. But I would almost
> rather them email a feature request.

We can always do that. Either the way we
have with the LKJ business, or by adding new
density functions.

I hope to make Stan extensible soon so that
users can define their own functions and add
them to the symbol table so they play nicely with
all the (basic) type checking.

- Bob

[Ben Goodrich]

That makes sense iff Sigma were known, but in this case where Sigma is
parameterized, why not pass Sigma to multi_normal() so that it can do
solve(Sigma) and det(Sigma) in a fast and safe way?

> Are you suggesting we extend Eigen with a CovMatrix class?

Maybe. I think eventually Eigen will probably extend itself with all
kinds of structured matrices, but that is not on their agenda right
now. So, we could "extend" Eigen's Matrix class (and possibly deal
with the consequences if and when Eigen makes its own CovMatrix class)
or we could make a totally new class called CovMatrix that "embeds"
one or more Eigen Matrices as private attributes so that Stan can
pretty much keep the same API even if Eigen changes in major ways.

[Ben Goodrich]

Hit send too soon ...

I do kind of like the design of the Matrix package in R (as opposed to
the matrix that is in the base R package, which is basically just a
vector of numbers with a dimensions attribute). You can see its design
for a correlation matrix, for example, with

library(Matrix)
help("corMatrix-class")

A corMatrix has a slot (attribute) called "factors", which is a list
that caches any matrix factorizations like Cholesky and so forth. So
when some high-level operation (like det()) is done to a corMatrix,
the low-level internals are free to carry that out via the Cholesky
factor.

Ben

[Bob Carpenter]

Because you have to do the same solve N times.
An alternative is to vectorize, and write:

y ~ multi_normal(mu,Sigma);

for vector y. Then multi_normal_log does the same unfolding
inside its loop. If the mu are determined by regressions,
then we also need to allow passing in a sequence of mu
vectors.

for (n in 1:N)

y[n] ~ normal(alpha + beta * x[n], sigma);

with alpha, beta and x[n] simple scalars, turns into

y ~ normal(alpha + beta * x, sigma);

where y is a vector, x is a vector and alpha and
beta are still scalars, with the expression defined so that

(alpha + beta * x)[n]
= alpha + beta * x[n].

so that the whole is equivalent to the original statement
looping over the N.

So far, we've only done this for the regular normal,
and even then, we haven't allowed varying means. It'll
be a much bigger saving with covariance matrices.

>> Are you suggesting we extend Eigen with a CovMatrix class?
>
> Maybe. I think eventually Eigen will probably extend itself with all
> kinds of structured matrices, but that is not on their agenda right
> now.

Doug Bates and some other folks have been pushing
things like symmetric matrix and easier handling
of positive definiteness

> So, we could "extend" Eigen's Matrix class (and possibly deal
> with the consequences if and when Eigen makes its own CovMatrix class)
> or we could make a totally new class called CovMatrix that "embeds"
> one or more Eigen Matrices as private attributes so that Stan can
> pretty much keep the same API even if Eigen changes in major ways.

I was thinking of whether we'd just do something
in Stan's compiled language, whether we'd try to
do something with a nice C++ API to build into our
matrix class, or whether we'd just do something with
Stan's own modeling language.

Either way, it makes sense to me to extend Matrix.
After all, a CovMatrix is a Matrix. That lets us apply
all of the usual operators to it. But given C++'s semantics,
that'd really want us to allow the assignment and
access referent operators. If it doesn't extend Matrix,
we'll have a ton of work to do defining all the operators
for it. But then maybe it won't have many.

These kinds of extensions seem to be on Eigen's agenda,
just not highly prioritized. They're willing to help
with high-level design advice, but need someone to code
and test.

We already plan to go down this road in the near future.
Eigen needs someone to push through all the heterogeneous
operator traits (so that you could multiply a Matrix<var,...>
and a Matrix<double,...>). That's what we need for much more
efficient matrix auto-diff operations. Right now, we're
coding around their expressions and have to over-promote
doubles to auto-dif variables, which is costly in terms of
both time and memory.

Currently, Eigen has an example for Matrix<complex,...> and
Matrix<double,...>, but it's hard coded for those types.
That opens the way to an intermediate possibility where we
just follow what they did for complex rather than trying to
generalize all their matrix and operator traits.

- Bob

[Bob Carpenter]

You made the ops clear in the previous messages, and
they make a lot of sense. What's nice is that
they'll automatically handle the kind of lifting of
operations out of loops I was talking about in the
last message (at the cost of a very inexpensive cache
check and without requiring the programmer to do
it themselves in Stan, which is doable now, but both
error prone and ugly).

Implementing a Matrix itself may be hard to do
consistently given Eigen's use of the
curiously recurring template pattern and the way
the types get lost after expression assignment.

Eigen's LLT and LDLT classes have the kind of design
Ben's talking about. They take a matrix in their
constructors, and allow access to the factors through a
method. They even have a kind of update function, but I
don't think they maintain a full cache-like behavior.
So we could follow that design and I think implement
everything Ben's been suggesting.

All in all, this should be doable in a way that's
natural enough that we could contribute it to Eigen.

- Bob

[Ben Goodrich]

On Jan 29, 3:28 pm, Bob Carpenter <c...@alias-i.com> wrote:
> On 1/29/12 3:01 PM, Ben Goodrich wrote:
> > That makes sense iff Sigma were known, but in this case where Sigma is
> > parameterized, why not pass Sigma to multi_normal() so that it can do
> > solve(Sigma) and det(Sigma) in a fast and safe way?
>
> Because you have to do the same solve N times.
> An alternative is to vectorize, and write:

Aha! Now I see. Yes, we should definitely avoid recalculating an
inverse or a determinant. And, we should try to avoid making users
write optimized .stan files. Ideally, they should just write readable
(and correct) .stan files and Stan should figure out how to do it
fast.

With covariance matrices (more so than other parameters), numerical
algorithms need "the" covariance matrix in various forms. Sometimes
--- but rarely --- you need the literal, textbookish, KxK form,
sometimes you need the SelfAdjointView, sometimes you need a
correlation matrix and standard deviations, sometimes you need a
Cholesky factor, sometimes you need the inverse, etc. And on top of
that, you don't want to store more than you need or being doing the
same autodiff twice, so you have various storeage schemes, dense vs.
sparse, etc. All of that would seem to argue for making a CovMatrix
class (extending Matrix is fine with me) that is capable of
implementing all this but mostly hides the implementation.

Since Doug already went down this road for the Matrix library in R,
I'm sure we could coordinate with him on how best to do it with Eigen.

Ben