Gamma reparametrization into mode / sd - Question on jacobian transformation

[Luke]

Hello, I'm loving Stan!

I have a question on jacobians within Stan.

A pattern is emerging on the sd sqrt(alpha/beta^2) and mode (alpha-1)/beta of a gamma distributed parameter.

Now I can reparametrize the gamma distribution for mode and sd but then I have to link mode and sd in, for example:

sd[n] ~ normal( ang_coeff * mode, sigma[n] );  equal to 

-----------------------------------------------------------------------------------------------

sqrt(alpha[n]/beta[n]^2) ~ normal(ang_coeff *  (alpha[n]-1)/beta[n], sigma);

-----------------------------------------------------------------------------------------------

if alpha and beta are vectors[N] how can adjust with the jacobian?

P.S. The plot above in terms of pure alpha and beta is 

Thanks

[Bob Carpenter]

You need the Jacobian of the f:R^2 -> R^2 function

f(alpha, beta) = sqrt(alpha / beta^2), (alpha - 1) / beta

if you're going to declare alpha and beta as parameters and
put a prior on sqrt(alpha / beta^2) and (alpha - 1) / beta.
I'd also recommend sqrt(alpha) / beta as it does less arithmetic
and should also be more stable for small beta.

Alternatively, you can declare a parameter

vector[N] sqrt_alpha_over_beta;

and define

sqrt_alpha_over_beta ~ ...

and then derive the inverse of f to get alpha and beta back.
Then you won't need a Jacobian as you're putting the distribution
right on the parameter.

- Bob

> On Mar 31, 2017, at 12:28 AM, Luke <mangiol...@gmail.com> wrote:
>
> Hello, I'm loving Stan!
>
> I have a question on jacobians within Stan.
>
> A pattern is emerging on the sd sqrt(alpha/beta^2) and mode (alpha-1)/beta of a gamma distributed parameter.
>
>
>
>
>
>
>

> Now I can reparametrize the gamma distribution for mode and sd but then I have to link mode and sd in, for example:
>
> sd[n] ~ normal( ang_coeff * mode, sigma[n] ); equal to
>
> -----------------------------------------------------------------------------------------------
> sqrt(alpha[n]/beta[n]^2) ~ normal(ang_coeff * (alpha[n]-1)/beta[n], sigma);
> -----------------------------------------------------------------------------------------------
>
> if alpha and beta are vectors[N] how can adjust with the jacobian?
>
>
> P.S. The plot above in terms of pure alpha and beta is
>
>
>
>
>
>

> Thanks
>
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[stefano]

Thanks Bob,

1)

just to make sure that is correct, I share a possible solution.

Considering that my priors will be

alpha ~ gamma( .., ..)

sqrt(alpha)/beta ~ normal( ang_coeff *  (alpha[n]-1)/beta[n], sigma )

a possible solution will be adjusting for Jacobian with two variables

for(n in 1:N) 

target += log(alpha[n] + 1) - log2 - 3*log(beta[n]) - 0.5 * log(alpha[n])

2)

it does not make sense to add 

beta ~ gamma( .., .. )

right? Or it could be done? In this case alpha and beta would have overlapping priors.

Thanks a lot! 

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ejEfPYdLryA/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

The formulas didn't come through.

You want to compute the Jacobian for

f : (alpha, beta) -> (alpha, sqrt(alpha) / beta)

Then you can put priors on alpha and on sqrt(alpha) / beta.
You probably wouldn't also want a prior on beta, but if
you include one, it just gets multiplied into the density (added
to the log density).

- Bob

> On Apr 4, 2017, at 10:40 PM, stefano <mangiol...@gmail.com> wrote:
>
> Thanks Bob,
>
> 1)
>
> just to make sure that is correct, I share a possible solution.
>
> Considering that my priors will be
>
> alpha ~ gamma( .., ..)
> sqrt(alpha)/beta ~ normal( ang_coeff * (alpha[n]-1)/beta[n], sigma )
>
> a possible solution will be adjusting for Jacobian with two variables
>
>
>
>
>

> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ejEfPYdLryA/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --

> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[stefano]

Thanks Bob,

sorry for the images.

the new Jacobian that I came up with is from (see also attachment, same image)

for which the absolute determinant is sqrt(x)/y^2 + 1/(2 sqrt(x) y)

For which log increment is

for(p in 1:P) for(g in 1:G) 

    target += log(2*alpha[p,g] + beta[p,g]) - (0.5*log(alpha[p,g])) - (2*log(beta[p,g])) - log(2);

I hope it makes sanse.

Thanks for your support!

> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ejEfPYdLryA/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ejEfPYdLryA/unsubscribe.

To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

You should be able to test that it recovers simulated
data. You'll have to invert the function you computed
the Jacobian of to do that, I think.

- Bob

> On Apr 9, 2017, at 2:13 AM, stefano <mangiol...@gmail.com> wrote:
>
> Thanks Bob,
>

> sorry for the images.
>
> the new Jacobian that I came up with is from (see also attachment, same image)
>

> <Capture.PNG>

>
> for which the absolute determinant is sqrt(x)/y^2 + 1/(2 sqrt(x) y)
>

> <Capture2.PNG>

> > > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> > > To post to this group, send email to stan-...@googlegroups.com.
> > > For more options, visit https://groups.google.com/d/optout.
> >
> > --
> > You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> > To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ejEfPYdLryA/unsubscribe.

> > To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
> >
> >
> > --
> > You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ejEfPYdLryA/unsubscribe.

> To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

> <Capture.PNG>

[Ste]

Thanks Bob