The generalized gamma is an attractive likelihood for modelling survival data because other common survival likelihoods (i.e. the lognormal, the exponential, the gamma, and the weibull) can be considered as special cases of the generalized gamma.

As such, I’ve implemented the `_lpdf`

and `_lcdf`

functions based on the `{flexsurv}`

implementation found here. Those functions are shown below.

```
real generalized_gamma_lpdf(real x, real mu, real sigma, real Q) {
real lpdf;
if (Q!=0) {
real y = log(x);
real w = (y-mu) / sigma;
real abs_q = abs(Q);
real qw = Q*w;
real qi = 1.0 / (Q*Q);
lpdf = -log(sigma * x) +
log(abs_q) * (1 - 2*qi) +
qi * (qw - exp(qw)) -
lgamma(qi);
} else {
lpdf = lognormal_lpdf(x | mu, sigma);
}
return lpdf;
}
real generalized_gamma_lcdf(real x, real mu, real sigma, real Q) {
real y = log(x);
real w = (y-mu) / sigma;
real lcdf = 0;
if (Q != 0) {
real qq = 1.0 / (Q*Q);
real expnu = exp(Q*w)*qq;
if (Q > 0) {
lcdf = gamma_lcdf(expnu | qq, 1);
} else {
lcdf = gamma_lccdf(expnu | qq, 1);
}
} else {
lcdf = lognormal_lcdf(x | mu, sigma);
}
return lcdf;
}
```

My implementation in agrees with the flexsurv well I have a small model, shown below, for a single set of parameters \mu, \sigma, and Q. I can recover these parameters fairly well with good sampler diagnostics. But every once in a while, the model will experience ~2-4 divergences.

```
functions {
#include generalized_gamma_functions.stan
}
data {
int n;
array[n] real x;
}
parameters {
real mu;
real<lower=0> sigma;
real Q;
}
model {
mu ~ student_t(30, 0, 1);
sigma ~ gamma(3, 3);
Q ~ student_t(30, 0, 1);
for(i in 1:n){
target += generalized_gamma_lpdf(x[i]|mu, sigma, Q);
}
}
generated quantities {
array[n] real x_ppc;
for(i in 1:n) {
x_ppc[i] = generalized_gamma_rng(mu, sigma, Q);
}
}
```

Increasing `adapt_delta`

to 0.99 reduces the frequency of divergent transitions, but I still get divergences from time to time. I simulated the data above a 250 times and plotted what parameter combinations lead to divergences. Here is that plot.

There seems to be an issue near Q=0. Looking at the density(below), I imagine the \vert Q \vert is the culpriate.

$$ f(x \mid \mu, \sigma, Q)=\frac{Q Q\left(Q^{-2}\right)^{Q^{-2}}}{\sigma x \Gamma\left(Q^{-2}\right)} \exp \left(Q^{-2}(Q w-\exp (Q w))\right) $$

## Question

What are some strategies to prevent the divergences?

One pragmatic solution would be to just use the log normal if I suspect Q\approx 0 and bound Q above/below in the `parameter`

block otherwise.

Are there any other approaches? Ideally I could use one density, but model stacking is something my team regularly does, so we could always fit a log normal survival model *and* a generalized gamma with bounded Q.

Open to suggestions.