Trouble fitting random walk plus exponential error

Modeling
,
1

I’m having trouble fitting a normal random walk with exponential errors, even on simulated data. I think (hope) I’m misunderstanding something about how to parameterize this so Stan likes it.
The motivation is, there’s a manufacturing line that produces and numbers 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 eventual, larger goal is that there are actually two production lines making different things, and when A finishes it will switch to making B, and the forecast for the run of B is what I’m really interested in, also the rate for line B changes beyond what a random walk would cover, but I’m starting here to not overwhelm myself)

I’m modeling the underlying, unobserved production date for each relic number i as a random walk starting at some initial time (production delay from announcement), each product is made after the previous one is finished and takes some number of days that varies around a mean rate, and the delivery times are production times plus shipping, which is bounded below at 0 and has a long upper tail, so in my fake data I said exponential.

\text{prod}_1 = \text{prod delay} \\ \text{prod}_{i+1} \sim N(\mu_\text{prod} + \text{prod}_i, \sigma_\text{prod}) \\ \text{deliv}_i = \text{prod}_i + \text{ship time}_i \\ \text{ship time}_i \sim \text{exponential}(\lambda_\text{ship})

I understand that there are some other ways I could approach this that might avoid some issues but I think I have OK reasons for preferring this over for example a simple linear model (production is close to sequential maybe more accurately small batches, this handles small changes in rate that happen in production and we can observe in the data without needing change points, also now that I’m in it I’d just like to learn where I’m going wrong fitting this model. But, if you have a strong alternative feel free to say that too.

I started with Stan’s exponentially modified normal since it’s exactly what I’m describing and made some good progress, I could recover parameters for simulated data and the posterior for production time looked reasonable at first, but it had problems: 1) I did something kind of goofy pulling out the production random walk from the exp_mod_normal part and making its sd 0.1, but every time I tried to put the random walk into that mean parameter it had tons of divergences and stopped working 2) when I increased the number of simulated products, something changed that I don’t understand and the random walk overfits the data and blows up very negative in the forecast. 3) I would like to understand how to do this for other error distributions (where you can’t just plug it into a mean shift parameter), so I’d prefer not to use exp_mod_normal
But here’s the first Stan model that tricked me into thinking this was going to be quick work

data {
  int<lower=0> N;
  int<lower=0> N_obs;
  int<lower=1> relic_obs[N_obs];
  int<lower=1> relic_unobs[N - N_obs];
  vector[N_obs] delivery_obs;
}

parameters {
  real<lower=0> prod_delay;
  real<lower=0> prod_mu;
  real<lower=0> prod_sigma;
  real<lower=0> ship_sigma;
  real<lower=0> ship_lambda;
  
  vector[N - N_obs] delivery_unobs;
  vector[N] production;
}

transformed parameters {
  vector[N] delivery;
  delivery[relic_obs] = delivery_obs;
  delivery[relic_unobs] = delivery_unobs;
}

model {
  prod_delay  ~ normal(20, 5);
  prod_mu     ~ scaled_inv_chi_square(10, 0.75);
  prod_sigma  ~ normal(0, 1);
  ship_sigma  ~ normal(0, 0.5);
  ship_lambda ~ normal(0, 0.5);
  
  production[1]   ~ normal(prod_delay, 0.1);
  production[2:N] ~ normal(prod_mu + production[1:(N-1)], prod_sigma);
  delivery ~ exp_mod_normal(production, 0.1, ship_lambda);
}

And a good and bad looking posterior for the production vector from each (the only difference here was more unobserved measurements–a little more sparse and a longer forecast period at the end

image
image445×238 29.8 KB

image
image448×240 44.3 KB

So in service of 3) I tried to separate out the shipping time as the difference of observed delivery and latent production but now I’m just getting divergent transitions. I’m guessing there’s a geometry issue because of a bad parameterization, but I’ve tried a bunch of things, at first I had a parameter vector for unobserved delivery time, but I think that’s not needed, I made this composite shipping time vector from observed delivery times and an unobserved shipping time vector, I tried parameterizing the random walk part with a ‘production standardized’ vector and shift+scale it later. And in a dark moment, I asked AI for help, when its suggestions even compiled, it seemed to just be changing things even more blindly than I am.

One thought that keeps coming to me is, the prior for the production rate (in days per product) was a pain because my domain knowledge says it’s around 11 per week and anything between 5 and 25 would be believable, but I have to be careful because in this ‘days per product’ parameterization anything near zero becomes a large number per week very fast. I’m happy with this prior, but maybe this is causing the problem (it was the same in the exp_mod_normal one though)
I’m out of ideas and could use some advice. Here is the last thing I’ve been tinkering with.

data {
  int<lower=0> N;
  int<lower=0> N_obs;
  int<lower=1> relic_obs[N_obs];
  int<lower=1> relic_unobs[N - N_obs];
  vector[N_obs] delivery_obs;
}

parameters {
  real<lower=0> prod_delay;
  real<lower=0> prod_mu;
  real<lower=0> prod_sigma;
  real<lower=0> ship_lambda;
  
  vector[N-1] production_after_first;
  vector[N-N_obs] delivery_unobs;
  vector<lower=0>[N-N_obs] shipping_unobs;
}

transformed parameters {
  vector[N] production;
  production[1] = prod_delay;
  production[2:N] = production_after_first;
  
  vector<lower=0>[N] shipping_time;
  shipping_time[relic_obs] = delivery_obs - production[relic_obs];
  shipping_time[relic_unobs] = shipping_unobs;
}

model {
  prod_delay  ~ normal(20, 5);
  prod_mu     ~ scaled_inv_chi_square(10, 0.75);
  prod_sigma  ~ normal(0, 1);
  ship_lambda ~ normal(0, 0.5);
  
  production[2:N] ~ normal(prod_mu + production[1:(N-1)], prod_sigma);
  shipping_time ~ exponential(ship_lambda);
}

And the core parameters and log prob (likely unhelpful since it’s divergent) plus the posterior for the production vector

image
image672×459 12.3 KB

image
image567×388 5.12 KB

Here’s some R code to generate some simulated data if you’d like to try it yourself

ex_data <- 
  tibble(
    relic = 1:150
  ) %>%
  mutate(prod_delay = 20,
         prod_error = rnorm(n = n(), 0.7, 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:100, size = 50, 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 <= 100),
            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]
)
2

I’m not sure I totally understand the data generating process. From what I gather, the production time is bounded below by zero and above by the delivery time. We have some observed production times and their corresponding delivery times. For the unobserved production and corresponding delivery times we model them through latent parameters. I’ve explicitly set those bounds below. Since your first observation is actually observed we can use that to pin down prod_mu and then model the first order random walk by the successive differences in production times. Noting that we want the differences between delivery and production times to follow the exponential distribution with rate parameter ship_lambda.

I did need to set adapt_delta to 0.9 to avoid one divergence.

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<lower=0> prod_mu;
  real<lower=0> prod_sigma;
  real<lower=0> ship_lambda;
  
  vector<lower=prod_mu>[N - N_obs] delivery_unobs;
  vector<lower=0, upper=delivery_obs>[N_obs] production_obs;
  vector<lower=0, upper=delivery_unobs>[N - N_obs] production_unobs;
}

transformed parameters {
  vector[N] production;

  production[relic_obs] = production_obs;
  production[relic_unobs] = production_unobs;
}
model {
  production[1]  ~ normal(20, 5);
  prod_mu     ~ scaled_inv_chi_square(10, 0.75);
  prod_sigma  ~ normal(0, 1);
  ship_lambda ~ normal(0, 0.5);

  delivery_obs[1] - production_obs[1] ~ normal(prod_mu, prod_sigma);
  production[2:N] - production[1:N - 1] ~ normal(prod_mu, prod_sigma);
  delivery_obs - production_obs ~ exponential(ship_lambda);
  delivery_unobs - production_unobs ~ exponential(ship_lambda);
}

fit_vs_observed
fit_vs_observed1350×825 82.2 KB

n.b. no AI was used for this analysis besides helping to plot the data.

3

That worked! And it also worked with a larger N (and smaller N_obs). I really appreciate the help.
I made some small changes to clean things up and reinsert the ‘prod_delay’ parameter cause I won’t always observe it, but in the simulated data sometimes we do. I left production observed in even though I’m not using right now because I think I can find some factory tours and communications where it’s clear what number is being shipped out. Here’s what I ended up using

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<lower=0> prod_delay;
  real<lower=0> prod_mu;
  real<lower=0> prod_sigma;
  real<lower=0> ship_lambda;
  
  vector<lower=0>[N-N_obs] delivery_unobs;

  vector<lower=0, upper=delivery_obs>[N_obs] production_obs;
  vector<lower=0, upper=delivery_unobs>[N-N_obs] production_unobs;
}

transformed parameters {
  vector[N] production;
  production[relic_obs] = production_obs;
  production[relic_unobs] = production_unobs;
  
  vector[N] delivery;
  delivery[relic_obs] = delivery_obs;
  delivery[relic_unobs] = delivery_unobs;
}
model {
  prod_delay  ~ normal(20, 5);
  prod_mu     ~ scaled_inv_chi_square(10, 0.75);
  prod_sigma  ~ normal(0, 2);
  ship_lambda ~ normal(0, 0.5);

  production[1] ~ normal(prod_delay, prod_sigma);
  production[2:N] - production[1:(N-1)] ~ normal(prod_mu, prod_sigma);
  delivery - production ~ exponential(ship_lambda);
}

I don’t fully understand why my previous one didn’t work, but I’ll keep studying the differences.
For fun, here’s that model with real data (it’s the Dune pinball machine from Barrels of Fun) Still some work to do, but it doesn’t look too bad (I can potentially use shipping location and known delays to address some of those very high shipping times, but I’m happy so far)

image
image878×536 8.87 KB

And the forecast of last production and shipping date also looks great (but the run ended early at 675 last week so it’s only like 2 weeks out. I also tried cutting it off before relic 600 but there was a production delay of about a month that shows up in the data here)