cyclic DAGs and coding BUGS models

[Bob Carpenter]

I'm starting to work through the BUGS models
that still aren't working.

The Jama River Valley model seems like a real hack
with mostly ordered variables with some play in
the middle. And it's encoded using a cyclic graphical
model. But I can reduce what's going on to a simple
example for the manual.

BUGS uses cyclic constraints to do things
like select two ordered values in the interval
(low,high), as in:

a ~ dunif(low,b);
b ~ dunif(a,high);

It's cyclic because a depends on b and b depends on a.

BUGS will get the wrong joint probability here,
but the conditionals are right, so it presumably
samples correctly.

The right joint probabilty function is:

p(a,b) = 1/(high - low) * 1/(high - a)

We can code this model in Stan as follows:

parameters {
real(low,high) a;
real(a,high) b;
}
model {
a ~ uniform(low,high);
b ~ uniform(a,high);
}

Daniel and I just filled up our board working
out that the transforms produce the right
multivariate Jacobian.

The transform for

real(L,H) z;

is to take an unbounded z_raw to

z = L + (H-L) * inv_logit(z_raw).

Thus the unbounded probabilities work out to:

p(a_raw,b_raw)

= 1/(high - low)

* 1/(high - a)

* (high - low) * inv_logit(a_raw) * (1-inv_logit(a_raw))

* (high - a) * inv_logit(b_raw) * (1-inv_logit(b_raw))

The first two terms are for the uniforms and the next two
terms for the transforms. The first one's redundant when
high and low are data, but not the second one. The last
two terms give you the Jacobian determinant (the full Jacobian
is lower triangular).

- Bob

[Andrew Gelman]

Yes, I recall that Bugs has an awkward way of expressing constraints. It's good that we can do this in a cleaner fashion. This will be a useful contribution of stan.
A

[Jiqiang Guo]

Do we get the joint of (a, b) correctly?

I suppose you meant that the joint of (a, b) is 

 p(a,b) = 1/(high - low) * 1/(high - a) I(low < a < b < high)

If it is correct, we then have conditional dsn of a given b as 

p(a|b) \propto 1/(high -a) I(low < a < b),

which is not as specified by a ~ dunif(low, b) [if this means a given b]. 

The joint of (a, b) could be derived by

p(a,b) = p(b|a) / (\int p(b|a) / p(a|b) db).

--

Jiqiang 

[Bob Carpenter]

On 4/26/12 9:23 PM, Jiqiang Guo wrote:
> Do we get the joint of (a, b) correctly?
>
> I suppose you meant that the joint of (a, b) is
> p(a,b) = 1/(high - low) * 1/(high - a) I(low < a < b < high)

Indeed -- I just assumed the support meeting the constraint.

> If it is correct, we then have conditional dsnof a given b as

> p(a|b) \propto 1/(high -a) I(low < a < b),
>
> which is not as specified by a ~ dunif(low, b) [if this means a given b].

That's why BUGS will get the wrong probability.

> The joint of (a, b) could be derived by
>
> p(a,b) = p(b|a) / (\int p(b|a) / p(a|b) db).

Intuitively this makes no sense, but that integral
sure does evaluate to 1/p(a).

But this can't possibly be what BUGS is doing.

Having said that, we're only going to get the probability
up to a proportion.

- Bob

[Bob Carpenter]

Jiqiang was right -- I was double counting
the constraint and shouldn't have included
the b ~ uniform(a,1) statement.

The right probability function for a uniform
distribution over a, b such that 0 < a < b < 1 is:

p(a,b) = 2 * I(0 < a && a < b && b < 1)

I generalized Stan's transforms to allow
params (not pushed yet), so that the following
model implements the constant density with
constraint 0 < a < b < 1.

parameters {
real(0,1) a;
real(a,1) b;
}
model {
}

That's it. Nothing in the model because the
density p(a,b) is constant. The parameters
take care of the constraint (and Jacobian).

We can use this same idiom for constraint encoding
for the Jama model from BUGS Vol 3.

I also coded the JAGS-style model (at least I
think I recall seeing this in JAGS) with sorting
rather than circular DAGs:

parameters {
real(0,1) t1;
real(0,1) t2;
}
transformed parameters {
real(0,1) a;
real(a,1) b;
a <- fmin(t1,t2);
b <- fmax(t1,t2);
}
model {
}


- Bob

[Andrew Gelman]

Perhaps "Jiqiang is always right" has to go in with "Matt3 is always right" as one of our guiding principles.

[Bob Carpenter]

On 4/27/12 4:19 PM, Andrew Gelman wrote:
> Perhaps "Jiqiang is always right" has to go in with
> "Matt3 is always right" as one of our guiding principles.

I think we can just simplify to "Bob is always wrong (about
stats)".

The samplers did just what I said, not what I meant.

At least I have enough sense to do independent tests.

Anyway, the problem's resolved and this means we can
do truncation and censoring with unknown censoring
points as well as enforce the kind of crazy
ordering constraints in the BUGS Vol 3 Jama ordinal
regression type model.

- Bob