Error-prone measurement of an exponential - posterior in the wrong place

Modeling
1

I am trying to simulate and fit a model to a model where: T \sim Exp(\lambda), T^{*} = Tv where v \sim Gamma(c,d), where c and d are known and T^{*} but not T is observed. My Stan model is:

data {
  int<lower=0> N;
  real<lower=0> tstar[N];
  real<lower=0> a;
  real<lower=0> b;
  real<lower=0> c;
  real<lower=0> d;
}

parameters {
  real<lower=0> lambda;   
  real<lower=0> v[N];      
}

model {
  lambda ~ gamma(a, b);
  for (n in 1:N) {
    real t = tstar[n]/v[n];
    t ~ exponential(lambda);  
  }
  for (n in 1:N) {
    v[n] ~ gamma(c,d);
  }
}

I simulate the data and pass it to Stan using:

library(rstan)

# sample size
N <- 1000

# rate parameter for underlying exponential
lambda_true <- 2
T <- rexp(N, rate = lambda_true)
tstar <- T*rgamma(N,10,10)

# Prepare the data for Stan
stan_data <- list(
  N = N,
  tstar = tstar,
  a = 1,
  b = 1,
  c = 10,
  d = 10
)

# Fit the updated model
fit <- stan(
  file = "mymodel.stan",
  data = stan_data,
  iter = 2000,
  chains = 4,
  warmup = 1000,
)

This runs fine, with no warnings from Stan. With the above N and parameter values, the posterior for lambda looks close to normal, but the mean (and most of the posterior) is always too high, and not centered at the true value (2 in the above). Despite trying for hours to figure out what is going wrong I’m stumped. If I make the prior for lambda very very strong, centered at 2, then the posterior eventually centers at 2, but I think it ought, with N=1000 to be pretty close even with a Gamma(1,1) prior on lambda.

I am almost certainly doing something very stupid, but I really can’t figure out what. If anyone has any insights, I’d be very grateful.

Thanks
Jonathan

2

If T\sim\mathrm{Exp}\left(\lambda\right) and T^{\star}=Tv then (treating v as constant) T^{\star}\sim\mathrm{Exp}\left(\frac{\lambda}{v}\right).
The model seems to work if you write the likelihood as

  for (n in 1:N) {
    tstar[n] ~ exponential(lambda/v[n]);
  }

The difference between these models is equal to target += log(v[i]); (it’s a Jacobian adjustment)

1 Like