Random walk with positive differences and missing data

Modeling
,
4

I’d approach as follows. It’s a bit confusing because I’ve left the unobserved and observed vectors separate but you could refactor a lot of this. Notes

  • You need to have a lower bound on the production delay because it doesn’t make sense for that to be negative.
  • A lognormal prior on production_diff doesn’t follow your data generating process coded in R as the difference between productions seemingly can be negative
  • you have to handle the case when the first index is observed vs unobserved
  • I think it makes more sense to model the latent unobserved delivery to get the correct bounds

This samples but has issues. You’ll want to play around with the priors and make sure your data generating process is correct.

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;

  array[N] int<lower=0, upper=1> is_obs;
}
parameters {
  real<lower=(is_obs[1] == 1 ? delivery_obs[1] : 0)> prod_delay;   // prod_delay should be >= 0
  real prod_mu;
  real<lower=0> prod_sigma;
  real<lower=0> ship_lambda;
  
  vector[N - N_obs] delivery_unobs_raw;
  vector[N - 1] production_diff_raw;
}

transformed parameters {
  vector[N] production;
  vector[N - 1] production_diff;
  vector[N - N_obs] delivery_unobs;
  vector[N] delivery;

    production[1] = prod_delay;

  {
    int obs_cnt = 1;
    int unobs_cnt = 1;

    if (is_obs[1] == 1) {
    delivery[1] = delivery_obs[1];
  } else {
    delivery_unobs[1] = lower_bound_jacobian(delivery_unobs_raw[1], prod_delay);
    delivery[1] = delivery_unobs[unobs_cnt];
    unobs_cnt += 1;
  }

  for (i in 2:N) {
    if (is_obs[i] == 1) {
        production[i] = lower_upper_bound_jacobian(production_diff_raw[i - 1], delivery_obs[obs_cnt] - production[i - 1],  delivery_obs[obs_cnt]);
        obs_cnt += 1;
    } else {
        delivery_unobs[unobs_cnt] = lower_bound_jacobian(delivery_unobs_raw[unobs_cnt], production[i - 1]);
        production[i] = lower_upper_bound_jacobian(production_diff_raw[i - 1], delivery_unobs[unobs_cnt] - production[i - 1], delivery_unobs[unobs_cnt]);
        unobs_cnt += 1;
    }
  }
  production_diff = production[2:N] - production[1:N - 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 ~ normal(prod_mu, prod_sigma);
  delivery_obs - production[relic_obs] ~ exponential(ship_lambda);
  delivery_unobs - production[relic_unobs] ~ exponential(ship_lambda);
}