0

I need to extend the following model with a generated quantities block, but I can't figure out how to do it:

import pystan
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
plt.style.use('ggplot')

data = pd.read_csv('unemployment.csv')

stan_code = """data {
  int <lower=0> T;
  vector[T] x;
  vector[T] y;
}

parameters {
  vector[T] u_err; //Slope innovation
  vector[T] v_err; //Level innovation 
  real beta;
  real <lower=0> s_obs;
  real <lower=0> s_slope;
  real <lower=0> s_level;
}

transformed parameters {
  vector[T] u; //Level
  vector[T] v; //Slope
  u[1] = u_err[1];
  v[1] = v_err[1];
  for (t in 2:T) {
    u[t] = u[t-1] + v[t-1] + s_level * u_err[t];
    v[t] = v[t-1] + s_slope * v_err[t];
  }
}

model {
  u_err ~ normal(0,1);
  v_err ~ normal(0,1);
  y ~ normal (u + beta*x, s_obs);
}"""

data_feed = {'y': data['unemployment.office'].values, 'x': np.zeros((data.shape[0], )), 'T': data.shape[0]}
sm = pystan.StanModel(model_code=stan_code)
fit = sm.sampling(data=data_feed, iter=1000)

I have made the following changes to the stan_code:

stan_code = """data {
  int <lower=0> T;
  vector[T] x;
  vector[T] y;
  int <lower=0> n_pred;
}

parameters {
  vector[T] u_err; //Slope innovation
  vector[T] v_err; //Level innovation
  real beta;
  real <lower=0> s_obs;
  real <lower=0> s_slope;
  real <lower=0> s_level;
}

transformed parameters {
  vector[T] u; //Level
  vector[T] v; //Slope
  u[1] = u_err[1];
  v[1] = v_err[1];
  for (t in 2:T) {
    u[t] = u[t-1] + v[t-1] + s_level * u_err[t];
    v[t] = v[t-1] + s_slope * v_err[t];
  }
}

model {
  u_err ~ normal(0,1);
  v_err ~ normal(0,1);
  y ~ normal (u + beta*x, s_obs);
}

generated quantities {
    vector[T+n_pred] u_pred;
    vector[T+n_pred] u_pred_err;
    vector[T+n_pred] v_pred;
    vector[T+n_pred] v_pred_err;
    vector[T+n_pred] y_pred;

    u_pred[1:T] = u;
    v_pred[1:T] = v;
    u_pred_err[1:T] = u_err;
    v_pred_err[1:T] = v_err;
    y_pred[1:T] = y;

    for (t in (T+1):(T+n_pred)) {
        u_pred_err[t] = normal_rng(0, 1);
        v_pred_err[t] = normal_rng(0, 1);
        u_pred[t] = u_pred[t-1] + v_pred[t-1] + s_level * u_err[t];
        v_pred[t] = v_pred[t-1] + s_slope * v_err[t];
        y_pred[t] = normal_rng(u_pred[t] + beta*x, s_obs);
    }
}

"""

But I keep on getting errors when compiling, let alone running the model.

Tom Kealy
  • 2,537
  • 2
  • 27
  • 45

1 Answers1

1

You should post the error message, which was

SYNTAX ERROR, MESSAGE(S) FROM PARSER:
Dimension mismatch in assignment; variable name = y_pred, type = real; right-hand side type = real[ ].
Illegal statement beginning with non-void expression parsed as
  y_pred[t]
Not a legal assignment, sampling, or function statement.  Note that
  * Assignment statements only allow variables (with optional indexes) on the left;
  * Sampling statements allow arbitrary value-denoting expressions on the left.
  * Functions used as statements must be declared to have void returns

 error in 'model50a15d294c93_gc' at line 52, column 8
  -------------------------------------------------
    50:         u_pred[t] = u_pred[t-1] + v_pred[t-1] + s_level * u_err[t];
    51:         v_pred[t] = v_pred[t-1] + s_slope * v_err[t];
    52:         y_pred[t] = normal_rng(u_pred[t] + beta*x, s_obs);
               ^
    53:     }
  -------------------------------------------------

This says explicitly that for the variable whose name is y_pred and whose type is real (scalar), you are trying to assign to it a right-hand side whose type is real[] (array). I think you meant to index the vector x:

y_pred[t] = normal_rng(u_pred[t] + beta * x[t], s_obs);
Ben Goodrich
  • 4,870
  • 1
  • 19
  • 18
  • I suppose the question I really have is - is this approach correct, creating a generated quantities with a transformed parameters block? – Tom Kealy Mar 27 '20 at 09:17
  • You can have both but any calls to `_rng` must be in the generated quantities block. Once you make `y_pred[t] = normal_rng(u_pred[t] + beta * x[t], s_obs);` you code uses legal syntax, but I have no idea if it does what you intended. – Ben Goodrich Mar 28 '20 at 01:32