Simulation and inference for a stochastic process with a random time of event in rstan

[SuPie]

Hi Martin,
Thank you for your reply!

@martinmodrak

Hi, the main problem is IMHO that conditioning on a parameter in the way you do on tau makes the posterior discontinuous and this breaks Stan. When you need to handle a discrete parameter (and τ\tau is a discrete parameter here) in Stan, you need to marginalize it out, this amounts to (very roughly):

Actually \tau is not supposed to be a discrete parameter. By \tau\sim\mathcal{U}inf(0,10), I meant to refer to the continuous uniform distribution. Sorry for not making that clear! (I did try to emphasise that \tau does not necessarily coincide with one of the time points t_k of the observations though) However, I do realize that the conditioning on \tau is weird, I just don’t know how to prevent it.
Do you think it would make sense to “split” \tau into two parameters by setting \tau = k \cdot\Delta t - s with discrete parameter k being the index of the first observed time point after \tau, continuous parameter s being the time interval between \tau and the observed time point t_k, and \Delta t being the known time step between observations; and then treating k in the way you have suggested?

@martinmodrak

Nevertheless even the rest of the model seems weird, so I would start by treating tau as known and making sure the model works in this case and only then move to marginalizing tau out.

In particular, the hard bounds on sigma are likely to be problematic, and you will likely get better results with a soft bound (e.g. real<lower=0> sigma and sigma ~ normal(0,1) or sigma ~ inv_gamma(something suitable) , which has low prior probability of both sigma < 0.01 and sigma > 2 but will not hinder sampling.

I had already tried whether inference for the remaining parameters works when \tau is assumed to be known. In that case, things are fine:
E.g. for the simulated data from my first post (where the parameters used for simulating were m_0=200, \theta=0.11, and \sigma=0.01), the Stan code stanmodel_stoch_exp_decay_with_known_release_time.stan (1.9 KB) and in R:

y_obs <- read.table("y_obs.txt")[[1]]
stanfit <- stan(file = "stanmodel_stoch_exp_decay_with_known_release_time.stan", 
                data = list(M=181, time_points=seq(from=0, to=30, length.out=181), 
                            y_obs=y_obs, tau=1.1, index_before_tau=7), 
                seed = 42424)

I obtained the following output:

##          mean se_mean   sd   2.5%    50%  97.5% n_eff  Rhat
## m0    201.439    0.03 1.93 197.73 201.45 205.29  3708 1.001
## theta   0.111    0.00 0.01   0.10   0.11   0.13  3638 1.000
## sigma   0.009    0.00 0.00   0.01   0.01   0.01   485 1.003

@martinmodrak

The way the process is described is also a bit weird - if your values are bounded by 0 from the bottom but have varying means, maybe the process can be better described by a log normal or gamma distribution?

Describing the transition probabilities of the process by the normal distribution comes from the Euler-Maruyama approximation for the solution of a stochastic differential equation which is the kind of process that I am trying to model, therefore, I would definitely like to keep that part the way it is.