Hi,
As some of you may have noticed I keep summing gamma distributed random variables… :)
Using the saddle point approximation (pinging expert @martinmodrak ) I have managed to create the likelihood of a sum of three gamma distribution with a common coefficient of variation, also shifted to the right. (I’m using the mean/coefficient of variation (mu/tau) parametrization as detailed here.)
Here, I have a minimal example that shows my latest problem. If the shift is a free parameter, I cannot recover the rest of the parameters (even if I fix most of them).
In this minimal example, I fixed two of the means of the gammas, and so I have left the coefficient of variation, one mean of a gamma, and the shift:
functions {
real t_trans(real t_raw, vector mu, real tau){
//transforms t from unconstrained space to be < 1/(mu *tau^2)
real mint = 1/(max(mu)* tau^2);
return - exp(t_raw) + mint;
}
vector Kpmx(vector y, // unknown: t
vector theta, // parameters
data real[] x_r, // data (real)
data int[] x_i) { // ignored
// equation: K' - x = 0
vector[3] mu = theta[1:3];
real tau = theta[4];
real x = x_r[1] - theta[5];
vector[1] z;
real t = t_trans(y[1], mu, tau); // transformed t
z[1] = -log(x) + log(-(3*mu[1]*mu[2]*mu[3]*t^2*tau^4 - 2*mu[1]*mu[2]*t*tau^2 - 2*mu[1]*mu[3]*t*tau^2 + mu[1] - 2*mu[2]*mu[3]*t*tau^2 + mu[2] + mu[3])/((mu[1]*t*tau^2 - 1)*(mu[2]*t*tau^2 - 1)*(mu[3]*t*tau^2 - 1)));
return z;
}
real gamma_3sum_lpdf(data real x_, real mu_1, real mu_2, real mu_3, real tau, real shift, data int[] debug) {
// saddle point approximation
// sum of 3 gamma dists with the same tau/cv
vector[5] theta = [mu_1, mu_2, mu_3, tau, shift]';
real t; // saddle point
// best guess for t based on the analytical solution
real guess = 0;// t_a(x, theta[1:3],theta[4]);
real t_raw; // unbounded saddlepoint
real x = x_ - shift; //unshifted distribuion
t_raw = algebra_solver(Kpmx,
[guess]', //guess
theta, //params
{x_},
debug)[1];
// gets t in the constrained space:
t = t_trans(t_raw, theta[1:3], theta[4]);
//debug:
if(debug[1]){
print("x=",x_,"; mu_1 = ", mu_1, "; mu_2 =", mu_2,"; mu_3 = " ,mu_3,"; tau = ",tau,"; t=,",t ,"shift=",shift );
}
return -t*x
- log(tau)
+ log1p(-mu_1*t*tau^2)*(-1/tau^2)
+ log1p(-mu_2*t*tau^2)*(-1/tau^2)
+ log1p(-mu_3*t*tau^2)*(-1/tau^2)
- log(.5*(3*mu_1^2*mu_2^2*mu_3^2*t^4*tau^8 + mu_1^2 + mu_2^2 + mu_3^2 + t^3*(-4*mu_1^2*mu_2^2*mu_3*tau^6 - 4*mu_1^2*mu_2*mu_3^2*tau^6 - 4*mu_1*mu_2^2*mu_3^2*tau^6) + t^2*(2*mu_1^2*mu_2^2*tau^4 + 4*mu_1^2*mu_2*mu_3*tau^4 + 2*mu_1^2*mu_3^2*tau^4 + 4*mu_1*mu_2^2*mu_3*tau^4 + 4*mu_1*mu_2*mu_3^2*tau^4 + 2*mu_2^2*mu_3^2*tau^4) + t*(-2*mu_1^2*mu_2*tau^2 - 2*mu_1^2*mu_3*tau^2 - 2*mu_1*mu_2^2*tau^2 - 2*mu_1*mu_3^2*tau^2 - 2*mu_2^2*mu_3*tau^2 - 2*mu_2*mu_3^2*tau^2)))
+ log(fabs(mu_1*t*tau^2 - 1))*2
+ log(fabs(mu_2*t*tau^2 - 1))*2
+ log(fabs(mu_3*t*tau^2 - 1))*2
- log(pi())/2
- log(2)/2;
}
}
data {
int N_obs;
real RT[N_obs];
}
transformed data {
int debug[1] = {0};
// real shift = .075;
}
parameters {
real<lower=0> mu;
real<lower=0, upper = min(RT)> shift;
real<lower=0> tau;
}
model {
target += normal_lpdf(mu|.1,.1);
target += normal_lpdf(tau|.22,.1);
target += normal_lpdf(shift|.075,.01);
for(n in 1:N_obs)
target += gamma_3sum_lpdf(RT[n]| mu,.122, .025,tau,shift, debug);
}
And I cannot recover the parameters at all.
These are my true values and generative process:
N_obs <- 5000
mu <- .1
sigma <- .22
shift <- .075
RT3 <- shift + EnvStats::rgammaAlt(N_obs, mu, sigma) + EnvStats::rgammaAlt(N_obs, .025, sigma)+
EnvStats::rgammaAlt(N_obs, .122, sigma)
And I get this
mgamma3 <- cmdstan_model("gamma3_shifted.stan") # I know 500 is low, but it doesn't matter here
f_gamma3 <- mgamma3$sample(data = list(RT = RT3, N_obs = N_obs), seed = 42, chains =3, parallel_chains = 3, iter_sampling = 500, max_treedepth = 15, adapt_delta = .9)
f_gamma3
> f_gamma3
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ 24020.19 24020.50 1.21 0.89 24017.99 24021.50 1.01 640 825
mu 0.01 0.01 0.00 0.00 0.01 0.01 1.00 627 743
shift 0.20 0.20 0.00 0.00 0.20 0.20 1.01 684 861
sigma 0.56 0.56 0.01 0.01 0.54 0.58 1.00 649 839
As you see everything is quite bad, with mu going to zero.
If I fix the shift parameter everything goes fine though and I can recover the exact values of the other two parameters. So I’m pretty sure it’s not a problem with my likelihood.
If I fix shift to 0.75, I get the exact values, (but weird NAs for ess_tail)
> f_gamma3
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
lp__ 22320.69 22321.00 1.00 0.74 22318.80 22321.70 1.01 715 NA
mu 0.10 0.10 0.00 0.00 0.10 0.10 1.00 1241 868
sigma 0.22 0.22 0.00 0.00 0.22 0.23 1.00 1099 1033
The most obvious thing to do, I think, would be to not let mu to get too close to 0. But how can I do that without changing lower to lower=.05? Any good prior that I could use? Or any other ideas?