Re: [stan-users] hierarchical mixture modelling

[Bob Carpenter]

I moved this discussion to dev.

Don't we need a Jacobian on the
(z,phi) to (x,y) transform?

Neither kappa nor mu are declared.

parameters {
real x;
real y;
}
transformed parameters {
real phi;
real z;
phi <- atan2(x, y);
z <- x*x + y*y;
}
model {
phi ~ von_mises(mu, kappa);
z ~ gamma(0.5, 0.5);
}

Now there's just the stumbling block of the
the Bessel functions involved if kappa is
a parameter.

- Bob

[Ben Goodrich]

On Wednesday, September 19, 2012 2:02:19 PM UTC-4, Bob Carpenter wrote:

Don't we need a Jacobian on the
(z,phi) to (x,y) transform?

Right. From the user's perspective it would be awesome if in

parameters {
   real x;
   real y;

      real mu<lower=-pi,upper=pi>;
      real<lower=0> kappa;

}
model {
   atan2(x, y) ~ von_mises(mu, kappa);
   x*x + y*y ~ gamma(0.5, 0.5);
}

that the parser would generate the code for the Jacobian automagically when it sees a function of parameters on the LHS of the ~. But I don't know how doable that is.

Ben

[Matt Hoffman]

Oh, right. I was thinking that was automatic in Stan, but of course it's not.

Happily, the Jacobian of the transformation turns out to be...2. (It's
easiest to work out in reverse, via the transformation
x=sqrt(z)cos(phi); y=-sqrt(z)sin(phi) and using the classic sin^2 +
cos^2=1 identity.) You could make it 1 easily by saying:

z <- 0.5 * (x*x + y*y);

but it's just a constant term regardless.

Matt

> --
>
>

[Bob Carpenter]

We can provide any number of built-in auto-transforms.

parameters {
real y;
}

transformed parameters {
real exp_y;

exp_y <- exp_transform(y);
}

which would compile to:

exp_y = exp_transform(y,lp__);

where we define exp_transform(y,lp__) by:

lp__ += y; // log Jacobian of exp(y) transform
return exp(y);

This would be a cinch to build. We already
have the transforms Stan itself uses implemented
this way. It only takes a bit of work on the parser
side to add the lp__ argument automatically.

What kind of transforms would be useful?

- Bob

[Ben Goodrich]

On Wednesday, September 19, 2012 10:22:23 PM UTC-4, Bob Carpenter wrote:

This would be a cinch to build.  We already
have the transforms Stan itself uses implemented
this way.  It only takes a bit of work on the parser
side to add the lp__ argument automatically.

What kind of transforms would be useful?

 
It seems like it would be hard to enumerate all useful transforms, but definitely having exp_transform() and circle_transform() would be better than making users do their own log absolute Jacobian determinants. Is it technically hard for the parser to recognize a general function of parameters on the LHS of a ~ and adjust lp__?

Ben

[Bob Carpenter]

It's easy to recognize that the LHS of ~ involves
a parameter. But then what?

We don't have second-order auto-dif, so while we
can compute the Jacobian of a general function,
we can't then roll that into the log probability
function because it needs to be differentiated again.

For the special cases, we've just built in
handwritten Jacobians, which we then auto-dif along
with the rest of the log probability function.

- Bob

[Ben Goodrich]

On Thursday, September 20, 2012 2:10:20 AM UTC-4, Bob Carpenter wrote:

On 9/20/12 12:19 AM, Ben Goodrich wrote:
> On Wednesday, September 19, 2012 10:22:23 PM UTC-4, Bob Carpenter wrote:
We don't have second-order auto-dif, so while we
can compute the Jacobian of a general function,
we can't then roll that into the log probability
function because it needs to be differentiated again.

That is the part I was missing. I was hoping that the Jacobian determinant could be evaluated and then the gradient of the log posterior evaluated afterward, but I guess that doesn't work.

Ben

[Bob Carpenter]

With second-order auto-dif we could do it. But then we
probably wouldn't be able to figure out if the Jacobian
was triangular and save ourselves a bunch of auto-dif work
in building the below-diagonal entries.

The other alternative is to supply a bunch of basic transforms,
including multivariate ones, with built-in Jacobians. Then figure
out how they could be legally composed (i.e., when we can legally
form the composition's Jacobian by multiplying the Jacobians
of the components).

- Bob