Hi! Thank you for your suggestions! I have implemented all three of them.
This is my new Stan model code:
// This Stan program defines a simple stochastic model of exponential decay where the
// initial time point tau is random.
// in this model, the index of the observation after tau is included as a parameter
// tau = index_after_tau_p * time_step - s
// and we have one set of parameters per possible value of index_after_tau_p
// non-centered parametrization
data {
int<lower=0> M; // number of observed time points
vector<lower=0>[M] y_obs; // observations
vector<lower=0>[M] time_points; // time points of the observations
real<lower=0> time_step; // time step between observations
int<lower=2> max_index; // upper bound for index of the observation after tau
int<lower=2, upper=max_index> min_index; // lower bound for index of the observation after tau
// the index after tau cannot be one (i.e. that of the first observation)
}
transformed data{
vector[M] log_y_obs;
real log_unif;
int n_poss_ind;
int n_random_states;
int x_start_index[max_index - min_index + 1];
log_y_obs = log(y_obs); // we assume log-normally distributed noise and therefore log-transform the data
n_poss_ind = max_index - min_index + 1; // number of possible values that the index after tau can take
log_unif = - log(n_poss_ind); // we assume a uniform prior for the index after tau;
n_random_states = (M - max_index + 1) * n_poss_ind + (n_poss_ind - 1) *
n_poss_ind / 2; // the number of random process states = length of x_flat
// calculate the index within x_flat where the values of the random process states
// start for the i-th possible value of index_after_tau
// i = index_after_tau - min_index + 1
for (i in 1:n_poss_ind) {
x_start_index[i] = (i - 1) * (M - i - min_index + 3) + (i - 1) * (i - 2) / 2 + 1;
}
}
parameters {
vector<lower=0, upper=time_step>[n_poss_ind] s; // time step from tau (at which
// initial amount m0 is released) to the next observed time point
vector<lower=0>[n_poss_ind] m0; // initial amount released at time tau
vector<lower=0>[n_poss_ind] theta; // decay rate
vector<lower=0>[n_poss_ind] sigma; // measurement noise parameter
vector[n_random_states] x_flat; // random process states
simplex[n_poss_ind] index_after_tau_p; //index_after_tau_p[i] == P(index_after_tau = i + min_index - 1)
}
transformed parameters {
vector[n_poss_ind] partial_lp;
vector<lower=0>[n_poss_ind] x[M]; // all process states
partial_lp = rep_vector(log_unif, n_poss_ind);
for (i in 1:n_poss_ind) { // again: i = index_after_tau - min_index + 1
//Compute the likelihood assuming a single given value of the index after tau
for (t in 1:(i + min_index - 2)) { // before tau x[t] = 0
x[t, i] = 0;
// observation model
partial_lp[i] += normal_lpdf(log_y_obs[t] | 0, sigma[i]);
}
// factor of the likelihood of the process for the time step from the state at
// time point tau to the state at the first observed time point
partial_lp[i] += normal_lpdf(x_flat[x_start_index[i]] | 0, 1);
x[i + min_index - 1, i] = m0[i] - theta[i] * m0[i] * s[i] +
sqrt(theta[i] * m0[i] * s[i]) * x_flat[x_start_index[i]];
// observation model for the first observation after tau
partial_lp[i] += normal_lpdf(log_y_obs[i + min_index - 1] |
log(x[i + min_index - 1, i] + 1), sigma[i]);
for (t in (i + min_index):M) { // (index_after_tau + 1):M
// factor of the likelihood of the process
partial_lp[i] += normal_lpdf(x_flat[x_start_index[i] + t - i - min_index + 1] | 0, 1);
x[t, i] = x[t-1, i] - theta[i] * x[t-1, i] * (time_points[t] - time_points[t-1]) +
sqrt(theta[i] * x[t-1, i] * (time_points[t] - time_points[t-1])) *
x_flat[x_start_index[i] + t - i - min_index + 1];
// observation model for the remaining observations after tau
partial_lp[i] += normal_lpdf(log_y_obs[t] | log(x[t, i] + 1), sigma[i]);
}
}
}
model {
//Weight the likelihood by the probability
target += log_sum_exp(log(index_after_tau_p) + partial_lp);
// priors
theta ~ normal(0, 5);
m0 ~ normal(300, 300);
sigma ~ normal(0,1);
}
When I now run (still providing the correct value of the index after \tau via min_index and max_index)
y_obs <- read.table("y_obs.txt")[[1]]
init_values <- c(rep(list(list(x_flat=as.array(rep(1, times= 174)))), times = 4))
stanfit <- stan(file = "stanmodel_stoch_exp_decay_with_random_release_time_incl_index_with_sep_param_x_non_centered.stan",
data = list(M=181, time_points=seq(from=0, to=30, length.out=181),
y_obs=y_obs, time_step=1/6, max_index=8, min_index=8),
seed = 42424,
control = list(max_treedepth = 15),
init = init_values)
there are no divergent transitions and the estimates look reasonable:
## mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
## m0[1] 201.798 0.03 2.31 197.71 201.68 206.61 4667 1.000
## theta[1] 0.111 0.00 0.01 0.10 0.11 0.13 1392 1.004
## sigma[1] 0.009 0.00 0.00 0.01 0.01 0.01 448 1.015
## s[1] 0.084 0.00 0.05 0.00 0.08 0.16 6319 1.000
(I’m initializing x_flat to ensure that the sampling starts in a region where the transformed parameters x are non-negative. I couldn’t find a better way to deal with the lower bound on x.)
Still, when I try to allow for several possible values for the index after \tau , e.g. max_index=9, min_index=7 , then this results in mostly divergent transitions and the chains looking messy:
y_obs <- read.table("y_obs.txt")[[1]]
M <- 181
max_index <- 9
min_index <- 7
length_x_flat <- (M - max_index + 1) * (max_index - min_index + 1) +
(max_index - min_index) * (max_index - min_index + 1) / 2
init_values <- c(rep(list(list(x_flat=as.array(rep(1, times = length_x_flat)))),
times = 4))
stanfit21 <- stan(file = "stanmodel_stoch_exp_decay_with_random_release_time_incl_index_with_sep_param_x_non_centered.stan",
data = list(M=M, time_points=seq(from=0, to=30, length.out=181),
y_obs=y_obs, time_step=1/6, max_index=max_index, min_index=min_index),
seed = 42424,
control = list(max_treedepth = 15),
init = init_values)
I am pretty sure that there is no error in the indexing of x or x_flat, as I’ve looked at the values of x at different iterations after warm up and those look reasonable for those chains that sample from a region around the true value (chains 1,2, and 4 for index_after_tau=8, and chain 3 for index_after_tau=9):
But now I don’t know what else to look at. Do have another idea what could still be causing these problems?