# Priors for hierarchical Gamma - Gamma model (Extended Pareto distribution)

**URL:** <https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418>\
**Category:** Modeling\
**Created:** [April 18, 2020, 11:41am UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418 "2020-04-18T11:41:31Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [April 18, 2020, 11:41am UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/1 "2020-04-18T11:41:31Z")

</div>

Hi everyone,

I hope you are doing well in this context. I have the following code for fitting an extended Pareto distribution, which can be built as a mixture of a Gamma and a rate parameter Gamma distributed.  
The resulting density is denoted as

f(y)=\frac{\Gamma(\alpha+\theta)}{\Gamma(\alpha)\Gamma(\theta)} \frac{1}{\beta} \frac{(y/\beta)^{\theta-1}}{(1+y/\beta)^{\alpha+\theta}}

where \theta is the shape of target Gamma and \alpha and \beta are the shape and rate of the mixing Gamma. Moreover,

E(Y)=\frac{\theta \beta}{\alpha-1}\\ sd(Y)=E(Y) \sqrt{\frac{\alpha+\theta-1}{\theta(\alpha-2)}}

A simple model,

```stan
data {
 int<lower=0> N;
 vector<lower=0>[N] y_i; 
}
parameters {
 real<lower=0> theta[N];
 real<lower=0> alpha;
 real<lower=0> gamma;
 real<lower=0> mu;
}
model { 
// target += normal_lpdf(alpha| 2,0.5);
// target += normal_lpdf(mu| 10,2);
// target += normal_lpdf(gamma| 100,50);
 target += gamma_lpdf(theta| mu,gamma);
 target += gamma_lpdf(y_i| alpha,theta);
}
generated quantities{
  real yrep[N];
  for (n in 1:N){
  yrep[n] = gamma_rng( alpha, gamma_rng(mu,gamma));
  }
}

```

I tested it with the following simulated data: rgamma(shape= 2 , rate = rgamma(n,shape= 10, rate = 100), n=10000)

How would you generally propose appropriate priors for this case?

With these priors  
target += normal\_lpdf(alpha| 2,0.5);  
target += normal\_lpdf(mu| 10,2);  
target += normal\_lpdf(gamma| 100,50);

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/d/d08187ed83cb1dae2b0d99bc8fd4c005561f26a9.png)

But with these  
target += normal\_lpdf(alpha| 0,5);  
target += normal\_lpdf(mu| 0,25);  
target += normal\_lpdf(gamma| 0,250);

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/1ca4a36aae664aca02d4f7521873ffd1dc7b607b.png)

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 18, 2020, 12:26pm UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/2 "2020-04-18T12:26:33Z")

</div>

You have

> [@Juan\_Ignacio\_de\_Oyarbide](#):
>
> parameters { real\<lower=0\> theta[N];  
> real\<lower=0\> alpha;  
> real\<lower=0\> gamma;  
> real\<lower=0\> mu;  
> }

which suggests you’re trying to estimate N + 3 parameters from N observations. I tried letting \theta be a scalar, but got into serious sampling problems. This doesn’t seem to be a very easy model to fit.

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [April 18, 2020, 12:34pm UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/3 "2020-04-18T12:34:45Z")

</div>

Yes. I’ll try to generate samples for lower values of \mu (the shape of the mixing) and see what happens.  
Generally, are the realizations of \theta considered as actual parameters?

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 18, 2020, 12:35pm UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/4 "2020-04-18T12:35:56Z")

</div>

> [@Juan\_Ignacio\_de\_Oyarbide](#):
>
> Generally, are the realizations of θ\theta considered as actual parameters?

“If you don’t see it, it’s a parameter”, Statistician, Bayesian. ;-)

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [April 19, 2020, 4:02pm UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/5 "2020-04-19T16:02:41Z")

</div>

I am trying this now, but the posteriors are strange. I may have some error in the formulae.

```stan
data {
 int<lower=0> N; //number of observations
 vector<lower=0>[N] y_i; //T 
}
parameters {
 real<lower=0> theta;
 real<lower=0> alpha;
 real<lower=0> gamma;
}
model { 
 target += normal_lpdf(alpha| 2,0.5);
 target += normal_lpdf(gamma| 10,5);
 target += normal_lpdf(theta| 100,30);
 
 target += log(lbeta(alpha,gamma))-log(theta)+(alpha-1)*log(y_i/theta)-(alpha+gamma)*log(1+y_i/theta);
 
}
generated quantities{
  //real yrep[N];
}

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/e9e195b09169c9446f25d3e840c0c879c54868b2.png)

f\_{X\_i}(x)=\int\_{\beta} f\_{X\_i}(x|\beta)f(\beta) d\beta\\ f\_{X\_i}(x)=\int\_{\beta} \frac{\beta ^ \alpha}{\Gamma(\alpha)} x^{(\alpha -1)} e^{- \beta x} \frac{\theta ^ \gamma}{\Gamma(\gamma)} \beta ^{(\gamma -1)} e^{- \theta \beta} d\beta\\ f\_{X\_i}(x)=\frac{x^{(\alpha -1) \theta ^ \gamma}}{\Gamma(\alpha) \Gamma(\gamma)} \int\_{\beta} \beta ^ \alpha e^{-\beta x} \beta ^{(\gamma -1)} e^{- \theta \beta} d\beta\\ f\_{X\_i}(x)=\frac{x^{(\alpha -1) \theta ^ \gamma}}{\Gamma(\alpha) \Gamma(\gamma)} \int\_{\beta} \beta ^ {\alpha+\gamma -1}e^{-\beta (x+\theta)} d\beta\\ f\_{X\_i}(x)=\frac{x^{(\alpha -1)} \theta ^ \gamma}{\Gamma(\alpha) \Gamma(\gamma)} \frac{\Gamma(\alpha + \gamma)}{(x+\theta)^{\alpha + \gamma}} \int\_{\beta} \beta ^ {\alpha+\gamma -1}e^{-\beta (x+\theta)} \frac{(x+\theta)^{\alpha + \gamma}}{\Gamma(\alpha + \gamma)} d\beta\\ f\_{X\_i}(x)=\frac{\Gamma(\alpha + \gamma)}{\Gamma(\alpha) \Gamma(\gamma)} \frac{x^{(\alpha -1)} \theta ^ \gamma}{(x+\theta)^{\alpha + \gamma}}\\ f\_{X\_i}(x)=\frac{\Gamma(\alpha + \gamma)}{\Gamma(\alpha) \Gamma(\gamma)} \frac{x^{(\alpha -1)} \theta ^ \gamma \theta ^ \alpha}{(x+\theta)^{\alpha + \gamma} \theta ^ \alpha}\\ f\_{X\_i}(x)=\frac{\Gamma(\alpha + \gamma)}{\Gamma(\alpha) \Gamma(\gamma)} \frac{x^{(\alpha -1)}}{\frac{(x+\theta)^{\alpha + \gamma} \theta ^ \alpha}{\theta ^ {\alpha + \gamma}}}

f\_{X\_i}(x)=\frac{\Gamma(\alpha+\gamma)}{\Gamma(\alpha)\Gamma(\gamma)} \frac{1}{\theta} \frac{(x/\theta)^{\alpha-1}}{(1+x/\theta)^{\alpha+\gamma}}

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 19, 2020, 4:43pm UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/6 "2020-04-19T16:43:59Z")

</div>

> [@Juan\_Ignacio\_de\_Oyarbide](#):
>
> log(lbeta(alpha,gamma))

Shouldn’t this be just `lbeta(...)`? You may also be missing a minus sign there.

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [April 26, 2020, 8:44am UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/7 "2020-04-26T08:44:50Z")

</div>

Hello, I got back to this. I missed the minus, but isn’t it log(lbeta)? Should I apply any kind of transformation to improve the sampling?

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 26, 2020, 10:25am UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/8 "2020-04-26T10:25:56Z")

</div>

`lbeta` returns the natural logarithm of the beta function.

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [April 29, 2020, 5:03pm UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/9 "2020-04-29T17:03:12Z")

</div>

I corrected those mistakes but still having problems with the sampling. Is it related to the tunning?

The log-increments are defined as

log(f(x\_i)) = -log(B(\alpha,\gamma))+log(\theta)-(\alpha-1)\times log(x\_i/\theta)+(\alpha+\gamma)\times log(1+x\_i/\theta)

```stan
data {
 int<lower=0> N; //number of observations
 vector<lower=0>[N] X; // 
}
parameters {
 real<lower=0> alpha;
 real<lower=0> gamma;
 real<lower=0> theta;
}
model { 
 target += normal_lpdf(alpha| 2,0.1);
 target += normal_lpdf(gamma| 10,1);
 target += normal_lpdf(theta| 100,5);
 
 target += -lbeta(alpha,gamma)+log(theta)-(alpha-1)*log(X/theta)+(alpha+gamma)*log(1+X/theta);
 
}
generated quantities{
  //real xrep[N];
}

```

There were 42 divergent transitions after warmup. Increasing adapt\_delta above 0.8 may help. See  
[http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupThere](http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupThere) were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See  
[http://mc-stan.org/misc/warnings.html#bfmi-lowExamine](http://mc-stan.org/misc/warnings.html#bfmi-lowExamine) the pairs() plot to diagnose sampling problems  
The largest R-hat is 2.85, indicating chains have not mixed.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#r-hatBulk](http://mc-stan.org/misc/warnings.html#r-hatBulk) Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#bulk-essTail](http://mc-stan.org/misc/warnings.html#bulk-essTail) Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.  
Running the chains for more iterations may help. See  
[http://mc-stan.org/misc/warnings.html#tail-ess](http://mc-stan.org/misc/warnings.html#tail-ess)

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/ec2ff6bd8d8460325f9c238896a0f29fc859eb91.png)

---

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [April 29, 2020, 5:52pm UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/10 "2020-04-29T17:52:05Z")

</div>

> [@Juan\_Ignacio\_de\_Oyarbide](#):
>
> target += -lbeta(alpha,gamma)+log(theta)-(alpha-1)\*log(X/theta)+(alpha+gamma)\*log(1+X/theta);

There’s a conspicuous absence of a loop in the line incorporating the likelihood. Also, just in case these are caused by numerical issues, try this implementation:

```stan
for (i in 1:N) target += -lbeta(alpha,gamma) + log(theta)-(alpha-1)*( log(X[i]) - log(theta)) + alpha + gamma)*log1p_exp(log(X[i]) - log(theta));

```

You might also move this line into its own function so you can re-use terms like `log(theta)` and `log(X[i])`.

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [April 30, 2020, 7:02am UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/11 "2020-04-30T07:02:49Z")

</div>

I am not sure a loop is necessary, I was using a vectorized notation. Anyways, I tried that way also with the transformations you suggested, but didn’t succeed.  
Isn’t it the shape of the beta function causing trouble in the likelihood?  
This density is challenging…

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/8/88f0609f407b58688024ccb8b3e3c1cae395ff10.png)

---

<div class="post-metadata">

**Author:** ![Juan\_Ignacio\_de\_Oyarbide](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/juan_ignacio_de_oyarbide/32/7447_2.png) [@Juan\_Ignacio\_de\_Oyarbide](https://discourse.mc-stan.org/u/Juan_Ignacio_de_Oyarbide)\
**Post date:** [May 1, 2020, 10:05am UTC](https://discourse.mc-stan.org/t/priors-for-hierarchical-gamma-gamma-model-extended-pareto-distribution/14418/12 "2020-05-01T10:05:51Z")

</div>

I am using the fact that

X\_i= \frac{\theta \alpha}{\gamma} \frac{G\_{\alpha}^{(i)}}{G\_{\gamma}}

where G\_{\alpha} and G\_{\gamma} are rv’s gamma distributed with mean 1. Then X\_i is an Extended Pareto as we defined before.  
I am defining those two latent variables and the sampling works, but I still have problems with wide priors. This case seems to be similar to the mixture.

```stan
data {
 int<lower=0> N; //number of observations
 vector<lower=0>[N] X; // 
}
parameters {
 real<lower=0> alpha;
 real<lower=0> gamma;
 real<lower=0> theta;
 vector<lower=0>[N] GA;
 }
transformed parameters{
 real<lower=0> GG;
 for (i in 1:N){
   GG= (theta * alpha ) * (GA[i] / X[i]);
 }
}
model { 
 target += normal_lpdf(alpha| 2,0.1);
 target += normal_lpdf(gamma| 10,0.1);
 target += normal_lpdf(theta| 100,5);
 
 target += gamma_lpdf(GA|alpha,alpha);
 target += gamma_lpdf(GG|gamma,1);
}
generated quantities{
  //real xrep[N];
}

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/a/aa574fa2cfcf832357b990461f5e7817615acccf.png)