I’m trying to modify a random walk model with exponential errors and missing observations to have differences that are also constrained to be positive. I’m struggling with the model because there’s something that feels circular about the combination of missing data and the constraints, every time I try to push down one part, another part pops up. I think there’s some other way to think about it that I’m missing.
I had previously asked for help with a normal random walk in a question here Trouble fitting random walk plus exponential error . The motivation is the same, there are two manufacturing lines that produce products sequentially (the sequential IDs are called ‘relic numbers’ in their terminology, so that’s what’s in the code below). I want to estimate the production rate and how it changes over time to forecast the end of the production run, but I can only observe some delivery dates and no production dates. The errors are positive and have a long tail (since they’re essentially shipping times), and the previous normal + exponential random walk worked well for individual lines.
Where it failed was, I’m modeling the production rate as time between sequential products, and the larger goal is that lines A and B make different products right now, but when A finishes its run it will switch to the same product as B and I really want to forecast the production time for the end of the run of B. My plan was to add the rates together as 1/(1/A + 1/B) = AB/A+B using sampled differences instead of having to work out the distribution of that (I spent a little time on that but this seems better) but when the differences can be negative as in the previous normal random walk, that falls apart, so I’d like to constrain them to be positive.
One problem may be that the production time has to be below the observed delivery time, but when there are long gaps between observations, the production differences overshoot the next observation. But I don’t know how to determine if that’s causing the issue, and the constraints from the previous model don’t work the same when when I’m declaring a parameter for the differences.
I changed the previous model to make the random walk differences production_diff a parameter directly so it can get the <lower=0> constraint. But then I have to write the transformed production vector so that it uses the difference, but then I don’t know how to make sure those values are below the delivery_obs values since the upper is really delivery - prev production. The next thing I’m going try is splitting up production_diff into observed and unobserved. I tried working with lookups for whether each index is observed or unobserved and which index of that respective vector it was from, then building the production vector as either prev+diff or observed-error, but I couldn’t get that to work either.
Here’s where I’ve landed so far. Lots of ‘rejecting initial value’ errors like this, which makes me I think I’m looking in the right place for the problem, but I don’t know how to fix it
Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: exponential_lpdf: Random variable[5] is -0.9674, but must be nonnegative! (in 'string', line 29, column 2 to column 66)
And here’s the model
data {
int<lower=0> N;
int<lower=0> N_obs;
array[N_obs] int<lower=1> relic_obs;
array[N-N_obs] int<lower=1> relic_unobs;
vector[N_obs] delivery_obs;
}
parameters {
real prod_delay;
real prod_mu;
real<lower=0> prod_sigma;
real<lower=0> ship_lambda;
vector<lower=0>[N-1] production_diff;
}
transformed parameters {
vector[N] production;
production[1] = prod_delay;
for (i in 2:N) {
production[i] = production[i-1] + production_diff[i-1];
}
}
model {
prod_delay ~ normal(20, 5);
prod_mu ~ normal(0, 1);
prod_sigma ~ normal(0, 2);
ship_lambda ~ normal(0, 0.5);
production_diff ~ lognormal(prod_mu, prod_sigma);
delivery_obs - production[relic_obs] ~ exponential(ship_lambda);
}
And some code for generating simulated to pass to rstan
ex_data <-
tibble(
relic = 1:200
) %>%
mutate(prod_delay = 20,
prod_error = rlnorm(n = n(), meanlog = log(0.5), sdlog = 1),
prod = prod_delay + cumsum(prod_error),
ship_error = rexp(n = n(), rate = 1/7),
deliv = prod + ship_error,
is_observed = row_number() %in% sample(1:150, size = 35, replace = FALSE))
ex_data %>%
ggplot(aes(x = relic)) +
geom_point(data = . %>% filter(!is_observed),
aes(y = deliv, color = is_observed)) +
geom_line(data = . %>% filter(relic <= 150),
aes(y = prod)) +
geom_point(data = . %>% filter(is_observed),
aes(y = deliv, color = is_observed)) +
scale_color_manual(guide = "none",
values = c("TRUE"= "#000000", "FALSE" = "#CCCCCC"))
stan_data <- list(
N = nrow(ex_data),
N_obs = sum(ex_data$is_observed),
relic_obs = ex_data$relic[ex_data$is_observed],
relic_unobs = ex_data$relic[!ex_data$is_observed],
delivery_obs = ex_data$deliv[ex_data$is_observed]
)