Gamma distribution - parameterization -- gamma distribution

General
1

Hello, I’m new to stan. I’m wondering how gamma distribution was parameterized. According to the manual, it is (alpha, beta) .
However, my simulation results seem to suggest otherwise. Please focus on the stddev estimates. Please bear with me as beginners, and it is very likely I
Stan Model:

data {
  int n;                   // sample size
  real xbar;               // sample mean
  real <lower=0> ss;       // sample sum-of-square
  int mID;                 // gamma parameterization: 1 - assume alpha, beta
                           //                         2 - assumed k, theta
                           //
}

parameters {
  real mu;                 // ped mean
  real <lower=0> stddev;   // std dev: ped
}

model {
  xbar  ~ normal(mu, stddev/sqrt(n));
  if (mID == 1)
    ss ~ gamma((n-1.)/2. , (2*stddev^2));
  else
  ss ~ gamma((n-1.)/2. , 1/(2*stddev^2));
}

####################  simulation results  ######################
## Assumed parameter values
N <- 100
MU <- 0
SIGMA <- 0.50

## sample
set.seed(1)
xbar <- rnorm(1, MU, SIGMA/sqrt(N))
ss <- rchisq(1, df=N-1)*(SIGMA^2)

## stan fits (assumed alpha-beta parameterization)
pCall <- list(n=N, ss=ss, xbar=xbar, mID=1)smp <- sampling(mdl, data=pCall)
smp
## Inference for Stan model: ex.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## mu      -0.03    0.00 0.10  -0.23  -0.10  -0.03   0.04   0.16  4000    1
## stddev   0.99    0.00 0.07   0.86   0.94   0.99   1.04   1.13  3544    1
## lp__   -14.68    0.02 0.98 -17.21 -15.09 -14.38 -13.97 -13.70  2032    1

pCall2 <- list(n=N, ss=ss, xbar=xbar, mID=2)
smp2 <- sampling(mdl, data=pCall2)

## stan fits (assumed k-theta parameterization)

## Inference for Stan model: ex.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
##          mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## mu      -0.03    0.00 0.05  -0.13  -0.06  -0.03   0.00   0.07  3077    1
## stddev   0.51    0.00 0.04   0.45   0.48   0.51   0.53   0.59  3102    1
## lp__   -14.67    0.03 1.01 -17.45 -15.04 -14.37 -13.96 -13.70  1526    1

[edit: escaped model code]

1 Like
2

Sorry for the late reply.

The manual defines the formulas for the parameterizations we use. In this case, alpha is the shape parameter and beta the inverse scale. The kernel of the log density for constant parameters is

log Gamma(y | a, b) = (a - 1) * y - b * y;

P.S. I escaped your code so the post doesn’t look like a ransom note (blame markdown—they like to guess what you mean and apply seemingly random formatting).

1 Like