Bad estimation of shift and location in a complex likelihood

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) {
// sum of 3 gamma dists  with the same tau/cv
vector[5] theta = [mu_1, mu_2, mu_3, tau, shift]';
// best guess for t based on the analytical solution
real guess = 0;// t_a(x, theta[1:3],theta[4]);
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?

1 Like

It also worries me that the model that doesn’t recover the true values has biased posteriors with CrI that are too precise and wrong; see my other general question here.

I looked into this, can confirm that I recover params well with `shift` fixed, but there is weird bias when `shift` is taken as a parameter. I have no idea what is the cause…

2 Likes