I would like to model a continuous-time (partly) stochastic process X(t) that has a jump at a random time point \tau and that is not directly observed.
The process is deterministic up to \tau
X(t) \equiv 0 \qquad \text{ for } t < \tau,
at \tau, it jumps to a random value m_0
X(\tau) = m_0,
and after \tau, it evolves as a Markov process with a transition probability that depends on a parameter \theta:
\left\lbrace X(t)\, {\large\vert}\, X(s) = x_s\right\rbrace \sim \mathcal{N}_{>0} \left(x_s - \theta \cdot x_s\cdot (t-s), \,\sqrt{\theta \cdot x_s\cdot (t-s)}\right)\qquad \text{ for } \tau\leq s < t.
I assume multiplicative measurement errors for the observations y_1,\dots,y_M at time points t_1,\dots, t_M:
\log(y_k)\sim\mathcal{N}(\log(X(t_k) + 1),\, \sigma)\qquad \text{ for } k=1,\dots M.
Moreover, I assume the following priors
\tau \sim \mathcal{Unif}(0,10), \quad m_0 \sim \mathcal{N}_{>0}(300,300), \quad \theta \sim \mathcal{N}_{>0}(0,5), \quad \sigma\sim\mathcal{Unif}(0.001,2),
where \mathcal{N}_{>0} denotes a truncated normal distribution. Note that the random time point \tau does not need to be one of the observed time point.
I am able to simulate the model in rstan with simulate_stoch_exp_decay_with_random_release_time_fix_param.stan (2.1 KB) :
library(rstan)
stansim <- stan(file = "simulate_stoch_exp_decay_with_random_release_time_fix_param.stan",
data = list(M=181, time_points=seq(from=0, to=30, length.out=181),
tau=1.1, m0=200, theta=.11, sigma=.01),
algorithm = "Fixed_param", seed = 4242, iter = 1, chains = 1)
y_obs <- extract(stansim, pars = "y_obs")[[1]][1,]
write.table(y_obs, file = "y_obs.txt", col.names=FALSE, row.names=FALSE)
Now, I would like to infer the four parameters \tau, m_0, \theta, \sigma from the simulated dataset y_obs.txt (3.0 KB).
Among others, Iāve tried the following Stan Code (stanmodel_stoch_exp_decay_with_random_release_time_v1.stan (2.1 KB) ):
data {
int<lower=0> M;
vector<lower=0>[M] y_obs; // observations
vector<lower=0>[M] time_points; // time points of the observations
}
// we assume log-normally distributed noise and therefore log-transform the data
transformed data{
vector[M] log_y_obs;
log_y_obs = log(y_obs);
}
parameters {
real<lower=0, upper=10> tau; // time point at which initial amount m0 is released
real<lower=0> m0; // initial amount released at time tau
real<lower=0> theta; // decay rate
real<lower=0.001, upper=2> sigma; // measurement noise parameter
real<lower=0> x_k[M];// random process states
}
transformed parameters{
real<lower=0> x[M]; // process states at observed time points
for (t in 1:M) {
if (time_points[t] < tau ){
x[t] = 0; // process states before tau are deterministic
}else{
if(time_points[t] == tau ){// tau is an observed time point
x[t] = m0;
}else{ // time_points[t] > tau
x[t] = x_k[t];
}
}
}
}
model {
// log-likelihood part 1: prod. of transition probabilities of the diffusion process
int first; // variable to check whether the first transition after tau is yet to come
first = 1;
for (t in 1:M) {
if(time_points[t] >= tau ){
if(time_points[t] == tau ){ // tau is an observed time point
first = 0;
}else{// time_points[t] > tau
if(first){// transtition prob. between tau and the first observed time point after tau
x[t] ~ normal(m0 - theta * m0 * (time_points[t] - tau),
sqrt(theta * m0 * (time_points[t] - tau)));
first = 0;
}else{// transitions prob. between observed time points after tau
x[t] ~ normal(x[t-1] - theta * x[t-1] * (time_points[t] - time_points[t-1]),
sqrt(theta * x[t-1] * (time_points[t] - time_points[t-1])));
}
}
}
}
// log-likelihood part 2: multiplicative normal noise
for (i in 1:M){
log_y_obs[i] ~ normal(log(x[i] + 1), sigma);
}
// priors
theta ~ normal(0, 5);
m0 ~ normal(300,300);
}
and in R:
y_obs <- read.table("y_obs.txt")[[1]]
stanfit1 <- stan(file = "stanmodel_stoch_exp_decay_with_random_release_time_v1.stan",
data = list(M=181, time_points=seq(from=0, to=30, length.out=181),
y_obs=y_obs),
seed = 42424)
This results in mostly divergent transitions (increasing adapt_delta did not help) and the chains just get stuck:
Iām suspecting that sampling with Stan is simply not possible for this kind of model due to the fact that the components of X would have to switch between a deterministic and a stochastic state depending on the value of \tau. But I am not completely sure. Maybe there is some technique that I have not come across yet.
Does anyone have an idea on how to resolves this issue or can confirm my conjecture that sampling with Stan is not possible?
