bivariate normal-gamma distribution in stan

[shota...@gmail.com]

Hi all,

Edit: the previous code attached is my older version and has a typo. It is now fixed.

   I have a trouble with bivariate normal-gamma distribution in stan.

I wrote a simple code to sample random numbers based on bivariate normla-gamma distribution for my applications.

The following is my stan code and compiled in Ubuntu 13.10.

transformed data
{
   real mu;
   real<lower=0> kappa;
   real<lower=1> alpha;
   real<lower=0> beta;
   mu     <- -4.3;
   kappa  <-  3.1;
   alpha  <-  1.1;
   beta   <-  0.5;
}
parameters
{
   real<lower=0> lambda;
   real x;        
}
transformed parameters
{
   real<lower=0> sigma;
   sigma <- 1/sqrt(lambda*kappa);
}
model
{
   // Normal-Gamma distribution
   increment_log_prob(   gamma_log(lambda,alpha,beta)
                       + normal_log(x,mu,sigma)
                     );
}

I ran the executable with 

./normal-gamma.exe sample thin=10 num_warmup=1000000 num_samples=10000000

and I got a slightly off value which is

lambda = 2.251 with MCSE 2.323e-3

which is supposed to be 2.200 according to wikipedia.

Could anyone point me what is wrong with my code? Thank you for your help and kindness.

Note: My notation is based on this paper.

[Bob Carpenter]

I would check that

a) the model's converging properly,

b) you're using the same parameterization of gamma (or should
that be inverse gamma?), and

c) whatever you're checking against is giving you the posterior
mean, not the posterior mode.

Could you include the Wikipedia reference?

Also, it's more efficient to use the sampling statements where
possible, replacing

> increment_log_prob( gamma_log(lambda,alpha,beta)
> + normal_log(x,mu,sigma)
> );

with

lambda ~ gamma(alpha,beta);
x ~ normal(mu,sigma);

The reason is that the sampling statement versions drop all
the unnecessary normalizing constants. Of course, if you need
those, then you need to use the increment_log_prob form.

Also, it's probably easier to break the increments apart,
as

increment_log_prob(gamma_log(lambda,alpha,beta));
increment_log_prob(normal_log(x,mu,sigma));

It's also a bit more efficient because increment_log_prob
uses an accumulator whereas the addition operator adds a new
node to the expression tree for auto-diff.

- Bob


On Mar 13, 2014, at 3:29 AM, shota...@gmail.com wrote:

> Hi all,

>
> I have a trouble with bivariate normal-gamma distribution in stan.
>
> I wrote a simple code to sample random numbers based on bivariate normla-gamma distribution for my applications.
>
> The following is my stan code and compiled in Ubuntu 13.10.
>
> transformed data
> {
> real mu;
> real<lower=0> kappa;
> real<lower=1> alpha;
> real<lower=0> beta;
> mu <- -4.3;
> kappa <- 3.1;
> alpha <- 1.1;
> beta <- 0.5;
> }
> parameters
> {
> real<lower=0> lambda;
> real x;
> }
> transformed parameters
> {
> real<lower=0> sigma;

> sigma <- sqrt(lambda*kappa);

> }
> model
> {
> // Normal-Gamma distribution
> increment_log_prob( gamma_log(lambda,alpha,beta)
> + normal_log(x,mu,sigma)
> );
> }
>
> I ran the executable with
>
> ./normal-gamma.exe sample thin=10 num_warmup=1000000 num_samples=10000000
>
> and I got a slightly off value which is
>
> lambda = 2.251 with MCSE 2.323e-3
>
> which is supposed to be 2.200 according to wikipedia.
>
> Could anyone point me what is wrong with my code? Thank you for your help and kindness.
>
> Note: My notation is based on this paper.
>

> --
> 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.
> For more options, visit https://groups.google.com/d/optout.
> <normal-gamma.stan>

[shota...@gmail.com]

Hi Bob,

Edit: fixed some typos.

   Thank you for your quick response.

The following uses Wikipedia notation for convenience (sorry for this confusion).

Here are links to wikipedia.

http://en.wikipedia.org/wiki/Normal-gamma_distribution

http://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution

a) the model's converging properly, 

I have run the model several times and I got an almost same answer and R_hat = 1.000, so I thought it converged well.

b) you're using the same parameterization of gamma (or should that be inverse gamma?), and 

For normal distribution with unknown mean and unknown precision, normal-gamma is a conjugate prior. For unknown mean and unknown variance, normal-inverse-gamma is a conjugate prior. They are actually just parameter transformations.

I am interested in normal-gamma, but I can switch to normal-inverse-gamma for covenience. 

c) whatever you're checking against is giving you the posterior mean, not the posterior mode. 

Yes, I checked againt mean not mode. I should have checked the mode, too. I will check it when I get back to my office.

Also, it's more efficient to use the sampling statements where 
possible, replacing 

>    increment_log_prob(   gamma_log(tau,alpha,beta) 
>                        + normal_log(x,mu,sigma) 
>                      ); 
with 
  tau ~ gamma(alpha,beta); 
  x ~ normal(mu,sigma); 

I understand your point, but I fully do not understand 

  tau ~ gamma(alpha,beta); 
  x ~ normal(mu,sigma); 

is correct or not (but both of them returned almost same answers). I saw many BUGS people said that the same approach works with normal-inverse-gamma in BUGS, and it returned the expected answer in BUGS. However, I still not sure it is really correct way in the Gibbs sampler because given x, Normal-Gamma (or Normal-Inverse-Gamma) is not Gamma (or Inverse-Gamma). Given x, Normal-Gamma is gamma times extra tau terms from normal distribution. I may misunderstand how BUGS works, tho.

Anyway, as you can see from Wikipedia Normal-Gamma, we can factorize Normal-Gamma with

NG(x,tau|alpha,beta,mu,lambda) = gamma(tau|alpha,beta) normal(x|mu,tau,lambda)

So, I reparameterized tau*lambda = 1/sigma^2 which is defined in

transformed parameters

{

   real<lower=0> sigma;

   sigma = 1/sqrt(lambda*tau);

}

and I took log resulted in

   increment_log_prob(   gamma_log(lambda,alpha,beta)
                       + normal_log(x,mu,sigma)

but this did not return a correct mean for lambda, but x is fine.

I should have tried to write them out explicitly without using gamma_log and normal_log.

I will check it when I get back to my office...

Thank you for your help and kindness.

Shota

2014年3月13日木曜日 20時33分02秒 UTC+9 Bob Carpenter:

[shota...@gmail.com]

It seems like I was wrong.

The model returned correct mode, median for precision and correct mean, median, and mode for x.

With different combination of alpha and beta, the mean is correct, so the model is just diverging with particular combinations of alpha and beta...

For those who are interested in normal-gamma, the attachment is my normal-gamma model in stan without reparameterization and with reparameterization(commented out).

Thank you for your help and kindness.

PS.

By the way, I tried to sample 1E8 and print.exe crached with std::bad_alloc(A possible bug in print.exe)

2014年3月13日木曜日 22時02分57秒 UTC+9 shota...@gmail.com:

Hi Bob,

Edit: fixed some typos.

   Thank you for your quick response.

The following uses Wikipedia notation for convenience (sorry for this confusion).

Here are links to wikipedia.

http://en.wikipedia.org/wiki/Normal-gamma_distribution

http://en.wikipedia.org/wiki/Normal-inverse-gamma_distribution

a) the model's converging properly, 

I have run the model several times and I got an almost same answer and R_hat = 1.000, so I thought it converged well.

b) you're using the same parameterization of gamma (or should that be inverse gamma?), and 

For normal distribution with unknown mean and unknown precision, normal-gamma is a conjugate prior. For unknown mean and unknown variance, normal-inverse-gamma is a conjugate prior. They are actually just parameter transformations.

I am interested in normal-gamma, but I can switch to normal-inverse-gamma for covenience. 

c) whatever you're checking against is giving you the posterior mean, not the posterior mode. 

Yes, I checked againt mean not mode. I should have checked the mode, too. I will check it when I get back to my office.

Also, it's more efficient to use the sampling statements where 
possible, replacing 

>    increment_log_prob(   gamma_log(tau,alpha,beta) 
>                        + normal_log(x,mu,sigma) 
>                      ); 
with 

  tau ~ gamma(alpha,beta); 
  x ~ normal(mu,sigma); 

I understand your point, but I fully do not understand 

  tau ~ gamma(alpha,beta); 
  x ~ normal(mu,sigma); 

is correct or not (but both of them returned almost same answers). I saw many BUGS people said that the same approach works with normal-inverse-gamma in BUGS, and it returned the expected answer in BUGS. However, I still not sure it is really correct way in the Gibbs sampler because given x, Normal-Gamma (or Normal-Inverse-Gamma) is not Gamma (or Inverse-Gamma). Given x, Normal-Gamma is gamma times extra tau terms from normal distribution. I may misunderstand how BUGS works, tho.

Anyway, as you can see from Wikipedia Normal-Gamma, we can factorize Normal-Gamma with

NG(x,tau|alpha,beta,mu,lambda) = gamma(tau|alpha,beta) normal(x|mu,tau,lambda)

So, I reparameterized tau*lambda = 1/sigma^2 which is defined in

transformed parameters

{

   real<lower=0> sigma;

   sigma = 1/sqrt(lambda*tau);

}

and I took log resulted in

   increment_log_prob(   gamma_log(lambda,alpha,beta)
                       + normal_log(x,mu,sigma)

but this did not return a correct mean for lambda, but x is fine.

I should have tried to write them out explicitly without using gamma_log and normal_log.

I will check it when I get back to my office...

Thank you for your help and kindness.

Shota

2014年3月13日木曜日 20時33分02秒 UTC+9 Bob Carpenter:

I would check that

[Daniel Lee]

Can you send me the file that caused print to crash?

[shota...@gmail.com]

Hi Daniel,

   It is 5GB! My virtual machine only has 4GB of memory, so print.exe attempted to allocate memory greater than available one.

print.exe should warn users this instead of throwing std::bad_alloc.

Shota

2014年3月14日金曜日 11時08分15秒 UTC+9 Daniel Lee:

[Bob Carpenter]

>
> On Thu, Mar 13, 2014 at 9:30 PM, <shota...@gmail.com> wrote:
> It seems like I was wrong.
>
> The model returned correct mode, median for precision and correct mean, median, and mode for x.
>
> With different combination of alpha and beta, the mean is correct, so the model is just diverging with particular combinations of alpha and beta...
>

Diverging as in not converging? If so, could you send us
examples of values?

> For those who are interested in normal-gamma, the attachment is my normal-gamma model in stan without reparameterization and with reparameterization(commented out).

Thanks for sharing.

> PS.
> By the way, I tried to sample 1E8 and print.exe crached with std::bad_alloc(A possible bug in print.exe)

The print methods should probably be beefed up so as not to
require everything to be read into memory. Otherwise, we're going
to hit memory limits. We should probably catch those and print
out-of-memory warnings anyway. I just created an issue to try to
print better warnings:

https://github.com/stan-dev/stan/issues/604

- Bob

[shota...@gmail.com]

Hi Bob,

   I have attached a simplest stan-code (alpha = 1.1 and beta = 0.5).

According to Wikipedia (http://en.wikipedia.org/wiki/Normal-gamma_distribution),  an analytic mean, mode and variance are given by

x:

  Mean:x

  Mode:x

  StdDev:sqrt(beta/(lambda*(alpha-1)))

tau:

  Mean:alpha/beta

  Mode:(alpha-0.5)/beta

  StdDev:sqrt(alpha/beta^2)

Case 1: mu = -4.3, alpha = 1.1, beta = 0.5

Expected Values:

x:

  Mean:-4.3

  Mode:-4.3

  StdDev:sqrt(50/31) ~ 1.270...

tau:

  Mean:2.2

  Mode:1.2

  Variance:4.4

Case 1 Results:

x:

  Mean: -4.301 (OK)

  Mode:-4.300 (OK)

  StdDev:1.033 (slightly off due to tau)

tau:

  Mean: 2.246 (2% off. considering MCSE = 5.262e-3, this error is a bit large)

  Mode: 1.200 (Exact)

  StdDev:2.133 (off a half an expected value)

Case 2:mu = -4.3, alpha = 9.0, beta = 4.0 (the code is not attached.)

Expected Values:

x:

  Mean:-4.3

  Mode:-4.3

  StdDev:sqrt(155)/31 ~ 0.401609

tau:

  Mean:2.25

  Mode:2.125

  StdDev:0.5625

Case 2 Results:

x:

  Mean:-4.300 (OK)

  Mode:-4.300 (OK)

  StdDev:0.3991 (close enough)

tau:

  Mean:2.251 (close enough)

  Mode:2.125 (OK)

  StdDev:0.7438 (a bit off)

To reproduce the results, I used the following arguments:

./normal-gamma-version1.exe sample num_warmup=10000 num_samples=1000000 random seed=324393 init=normal-gamma-version1.init.R

As a note, a gamma distribution with alpha = 1.1 and beta = 0.5 is a heavy-tale distribution (this could be a reason). 2% error seems to be small, but the normal-gamma is a conjugate prior to the lognormal distribution and 2% error is deadly significant to my case... 

[Michael Betancourt]

Raw numbers are meaningless by themselves; at the very least

you should be looking at (mean - expected) / MCSE.  And if you

expect that a heavy tail is possibly to blame then increase the

target acceptance probability.

--
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.
For more options, visit https://groups.google.com/d/optout.

<normal-gamma-version1.init.R><normal-gamma-version1.stan>

[shota...@gmail.com]

Hi Michael,

   Increasing adaptation target acceptance statistic (delta=0.90) gave me a better answer than the default delta.

Thank you for your help.