Vectorized model yielding different results to (original) non-vectorized one

Modeling
1

Hi! I have a working model that I’m trying to vectorize. This is the original (working) model:

data {
    int<lower=0> N_w16_obs;
    int<lower=0> N_w16_mis;
    vector[N_w16_obs] y_w16_obs;
    int<lower=0> N_oct_obs;
    int<lower=0> N_oct_mis1;
    int<lower=0> N_oct_mis2;
    vector[N_oct_obs] y_oct_obs;
    int<lower=0> N_nov_obs;
    int<lower=0> N_nov_mis1;
    int<lower=0> N_nov_mis2;
    vector[N_nov_obs] y_nov_obs;
    int<lower=0> N_dec_obs;
    int<lower=0> N_dec_mis1;
    int<lower=0> N_dec_mis2;
    vector[N_dec_obs] y_dec_obs;
    int<lower=0> N_jan_obs;
    int<lower=0> N_jan_mis1;
    int<lower=0> N_jan_mis2;
    vector[N_jan_obs] y_jan_obs;
    int<lower=0> N_feb_obs;
    int<lower=0> N_feb_mis1;
    int<lower=0> N_feb_mis2;
    vector[N_feb_obs] y_feb_obs;
    int<lower=0> N_mar_obs;
    int<lower=0> N_mar_mis;
    vector[N_mar_obs] y_mar_obs;
}
parameters {
    real rho0;
    real rho1;
    real<lower=0> sigma_w16;
    real<lower=0> sigma_diff;
    real<lower=0> sigma_y;
    real<lower=0> mu_diff;
    vector[N_w16_mis] y_w16_mis;
    vector[N_oct_mis1] y_oct_mis1;
    vector[N_oct_mis2] y_oct_mis2;
    vector[N_nov_mis1] y_nov_mis1;
    vector[N_nov_mis2] y_nov_mis2;
    vector[N_dec_mis1] y_dec_mis1;
    vector[N_dec_mis2] y_dec_mis2;
    vector[N_jan_mis1] y_jan_mis1;
    vector[N_jan_mis2] y_jan_mis2;
    vector[N_feb_mis1] y_feb_mis1;
    vector[N_feb_mis2] y_feb_mis2;
    vector[N_mar_mis] y_mar_mis;
}
transformed parameters {
    vector[N_w16_obs + N_w16_mis] y_w16 = append_row(y_w16_obs, y_w16_mis);
    vector[N_w16_obs + N_w16_mis] y_oct = append_row(append_row(y_oct_mis1, y_oct_obs), y_oct_mis2);
    vector[N_w16_obs + N_w16_mis] diff_oct = y_w16 - y_oct;
    vector[N_w16_obs + N_w16_mis] y_nov = append_row(append_row(y_nov_mis1, y_nov_obs), y_nov_mis2);
    vector[N_w16_obs + N_w16_mis] diff_nov = y_w16 - y_nov;
    vector[N_w16_obs + N_w16_mis] y_dec = append_row(append_row(y_dec_mis1, y_dec_obs), y_dec_mis2);
    vector[N_w16_obs + N_w16_mis] diff_dec = y_w16 - y_dec;
    vector[N_w16_obs + N_w16_mis] y_jan = append_row(append_row(y_jan_mis1, y_jan_obs), y_jan_mis2);
    vector[N_w16_obs + N_w16_mis] diff_jan = y_w16 - y_jan;
    vector[N_w16_obs + N_w16_mis] y_feb = append_row(append_row(y_feb_mis1, y_feb_obs), y_feb_mis2);
    vector[N_w16_obs + N_w16_mis] diff_feb = y_w16 - y_feb;
    vector[N_w16_obs + N_w16_mis] y_mar = append_row(y_mar_mis, y_mar_obs);
    vector[N_w16_obs + N_w16_mis] diff_mar = y_w16 - y_mar;
}
model {
    rho0 ~ normal(0, 1);
    rho1  ~ normal(0, 1);
    sigma_w16 ~ gamma(0.01, 0.01);
    sigma_diff ~ gamma(0.01, 0.01);
    sigma_y ~ gamma(0.01, 0.01);
    mu_diff ~ gamma(0.01, 0.01);
    for (i in 2:(N_w16_obs+N_w16_mis))
        y_w16[i] ~ normal(rho0 + rho1 * y_w16[i-1], sigma_w16);
    diff_oct ~ normal(mu_diff, sigma_diff);
    y_oct ~ normal(y_w16 + diff_oct, sigma_y);
    diff_nov ~ normal(mu_diff, sigma_diff);
    y_nov ~ normal(y_w16 + diff_nov, sigma_y);
    diff_dec ~ normal(mu_diff, sigma_diff);
    y_dec ~ normal(y_w16 + diff_dec, sigma_y);
    diff_jan ~ normal(mu_diff, sigma_diff);
    y_jan ~ normal(y_w16 + diff_jan, sigma_y);
    diff_feb ~ normal(mu_diff, sigma_diff);
    y_feb ~ normal(y_w16 + diff_feb, sigma_y);
    diff_mar ~ normal(mu_diff, sigma_diff);
    y_mar ~ normal(y_w16 + diff_mar, sigma_y);
}

and my attempt at vectorization:

data {
    int<lower=0> N_w16_obs;
    int<lower=0> N_w16_mis;
    vector[N_w16_obs] y_w16_obs;
    int N[3,6];                         // 1: mis1, 2: obs, 3: mis2
    real y_month_obs[65,6];
}
parameters {
    real rho0;
    real rho1;
    real<lower=0> sigma_w16;
    row_vector<lower=0>[6] sigma_diff;
    row_vector<lower=0>[6] sigma_y;
    row_vector<lower=0>[6] mu_diff;
    vector[N_w16_mis] y_w16_mis;
    real y_month_mis1[167,6];
    real y_month_mis2[104,6];
}
transformed parameters {
    vector[N_w16_obs + N_w16_mis] y_w16 = append_row(y_w16_obs, y_w16_mis);
    matrix[N_w16_obs+N_w16_mis,6] y_month;
    matrix[N_w16_obs+N_w16_mis,6] diff_month;
    for (n in 1:6) {
        int M1 = N[1,n];
        int M2 = N[3,n];
        y_month[:,n] = append_row(append_row(to_vector(y_month_mis1[1:M1,n]), to_vector(y_month_obs[:,n])), to_vector(y_month_mis1[1:M2,n]));
        diff_month[:,n] = y_w16[n] - y_month[:,n];
    }
}
model {
    rho0 ~ normal(0, 1);
    rho1  ~ normal(0, 1);
    sigma_w16 ~ gamma(0.01, 0.01);
    for (n in 1:6) {
        sigma_diff[n] ~ gamma(0.01, 0.01);
        sigma_y[n] ~ gamma(0.01, 0.01);
        mu_diff[n] ~ gamma(0.01, 0.01);
    }
    for (n in 2:(N_w16_obs+N_w16_mis))
        y_w16[n] ~ normal(rho0 + rho1 * y_w16[n-1], sigma_w16);
    for (n in 1:6) {
        diff_month[:,n] ~ normal(mu_diff[n], sigma_diff[n]);
        y_month[:,n] ~ normal(y_w16[n] + diff_month[:,n], sigma_y[n]);
    }
}

Both models run and converge but the graphs that I go on to plot afterwards look very different. Grateful for some advice.

2

I’m still a beginner and I’m probably missing something here, but is there a reason why you don’t include the effect of month as a parameter to be estimated? E.g. something like:

parameters {
   vector[n_months - 1] beta_month;
}

Also, describing your data might go a long way in helping others diagnose your problem and help you out!

3

Thanks for the suggestion @abartonicek. I’ve found a few things wrong in the code. The reason I’m getting different results is down to pooling of the parameters. I also think the way I’ve dealt with a ragged matrix is incorrect in the vectorized version.