Update:
Silly me! It is the log(current_history) when current_history == 0, which is of course -Inf, that was causing the problem. It is “corrected” by the exp (-current_history)=exp(Inf) that is 0, which is why the final numbers were correct (and why R never complained). However, having that negative infinity, even it was exponentiated to 0 later on, apparently, still breaks the gradient computation in some (subtle?) ways.
Original post:
I am new to Stan and I am trying to program a simple history integrator. And I am REALLY STUCK! Any help, hint, suggestion is highly appreciated. I suspect that it is because of local variables that I use in the loop but I have not a slightest clue of what am I doing so terribly wrong. :(
Here is the small data frame I use for the example.
test_df <- tibble(
Duration = c(3.181, 2.553, 2.552, 2.976, 3.094, 3.517, 4.764, 2.894, 2.718, 1.435, 2.105, 2.799, 5.398, 2.84699999999999, 3.48099999999999, 3.058, 4.023, 2.599, 2.552, 2.341),
IsOn = c(0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1),
Tmean = 2.5213
)
The history starts at 0 state. Whenever the state IsOn, it grows as negative exponential to 1, whenever it is not IsOn it drops as negative exponential. The grow/drop speed is controlled via a single parameter norm_tau, so that actual grow scale is tau = norm_tau * Tmean to match the average duration. Here is a simple implementation of it in R for norm_tau=1:
analytic_history <- function(df, norm_tau){
df$History <- NA
tau <- norm_tau * df$Tmean[1]
current_history <- 0
for(iT in 1:nrow(df)){
df$History[iT] <- current_history
if (df$IsOn[iT]){
current_history <- 1 - exp(-(df$Duration[iT] - tau * log(1 - current_history)) / tau)
}
else {
current_history <- exp(-(df$Duration[iT] - tau * log(current_history)) / tau)
}
}
df
}
analytical_df <- analytic_history(test_df, 1)
analytical_df
As a first simple step, I wanted to implement the same algorithm in Stan and, as a test, let it recover the norm_tau parameter. Here is the same code implemented in Stan:
data{
int<lower=1> N;
real<lower=0> Duration[N];
real<lower=0,upper=1> IsOn[N];
real<lower=0,upper=1> History[N];
real Tmean;
}
parameters{
real<lower=0.9, upper=1.1> norm_tau;
real<lower=0> sigma;
}
transformed parameters{
real<lower=0, upper=1> Predictions[N];
{
real tau = norm_tau * Tmean;
real current_history = 0;
for(iT in 1:N){
Predictions[iT] = current_history;
if (IsOn[iT] > 0){
current_history = 1 - exp(-(Duration[iT] - tau * log(1 - current_history)) / tau);
}
else {
current_history = exp(-(Duration[iT] - tau * log(current_history)) / tau);
}
}
}
}
model{
norm_tau ~ normal(1, 0.1);
sigma ~ exponential(1);
History ~ normal(Predictions, sigma);
}
And here is the sampling call
test_data <- list(
N = nrow(analytical_df),
Duration = analytical_df$Duration,
IsOn = analytical_df$IsOn,
History = analytical_df$History,
Tmean = analytical_df$Tmean[1]
)
test_init <- function(){
list(norm_tau = 1, sigma = 0.1)
}
test_fit <- rstan::sampling(single_state_history,
chains=1,
data=test_data,
init=test_init)
However, if I run it, I get
Chain 1: Rejecting initial value:
Chain 1: Gradient evaluated at the initial value is not finite.
Chain 1: Stan can’t start sampling from this initial value.
Chain 1:
Chain 1: Initialization between (-2, 2) failed after 100 attempts.
Chain 1: Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] “Error in sampler$call_sampler(args_list[[i]]) : Initialization failed.”
[1] “error occurred during calling the sampler; sampling not done”
If I fix tau = Tmean in the code (throwing out **norm_tau * ** bit, sampling works and the numbers for Predictions match those of History exactly (so, the code does compute the same thing). It makes no difference whether I init norm_tau to 1 (as I do in the code above) or to “random”. It makes no difference whether I restrict norm_tau (in the code above to be between 0.9 and 1.1) or just state <lower=0> or let it be unrestricted and use exp(norm_tau). I’ve also printed out values during the sampling, all finite, within a reasonable range or matching History values exactly, if initializing to norm_tau == 1.
To recap, I am stuck as I cannot understand what sort of diagnostics should I look at and what could I possibly be doing wrong. Any help is highly appreciated!