regression.stan

Go to the documentation of this file.

 
real intercept_centring_offset(vector means_fixed, vector beta,
                               vector beta_sd, matrix rdesign, int neffs) {
  if (neffs == 0) return 0.0;
  vector[1 + num_elements(beta_sd)] ext_beta_sd = append_row(1.0, beta_sd);
  vector[neffs] scaled_beta = beta .* (rdesign * ext_beta_sd);
  return dot_product(means_fixed, scaled_beta);
}
 
tuple(real, real) centring_offsets(
  vector means_fixed, vector beta, vector beta_sd, matrix rdesign, int neffs,
  int centre,
  int arima_present, int arima_T, int arima_G,
  int arima_p, int arima_d, int arima_q,
  matrix arima_z, vector arima_pacf, vector arima_theta,
  array[] real arima_sigma,
  int gp_present, int gp_T, int gp_G, int gp_M, real gp_L,
  int gp_type, real gp_nu, int gp_d, matrix gp_PHI, matrix gp_eta,
  array[] real gp_rho, array[] real gp_alpha
) {
  real design = intercept_centring_offset(
    means_fixed, beta, beta_sd, rdesign, neffs
  );
  real latent = arima_latent_mean_offset(
    arima_present, centre, arima_T, arima_G, arima_p, arima_d, arima_q,
    arima_z, arima_pacf, arima_theta, arima_sigma
  ) + gp_latent_mean_offset(
    gp_present, centre, gp_T, gp_G, gp_M, gp_L, gp_type, gp_nu, gp_d,
    gp_PHI, gp_eta, gp_rho, gp_alpha
  );
  return (design, latent);
}
 
vector regression_predictor(
  array[] real intercept, vector beta, int nobs, int neffs,
  matrix fdesign,
  tuple(vector, array[] int, array[] int) sparse,
  vector beta_sd, matrix rdesign,
  int add_intercept, int sparse_design, int centre,
  int arima_present, int arima_T, int arima_G,
  int arima_p, int arima_d, int arima_q, int arima_n_obs,
  matrix arima_z, vector arima_pacf, vector arima_theta,
  array[] real arima_sigma,
  array[] int arima_flat_idx,
  int gp_present, int gp_T, int gp_G, int gp_M, real gp_L,
  int gp_type, real gp_nu, int gp_d, matrix gp_PHI, matrix gp_eta,
  array[] real gp_rho, array[] real gp_alpha,
  array[] int gp_flat_idx
) {
  // Short-circuit when neither latent term is present: return
  // combine_effects() directly so we skip materialising the
  // intermediate `base` vector and the apply_* call frames on every
  // gradient evaluation. Behaviour is identical when a term is present.
  // Design centring (the `- offset` shift) is a constant applied by the
  // caller, so it does not enter here.
  if (!arima_present && !gp_present) {
    return combine_effects(
      intercept, beta, nobs, neffs, fdesign, sparse, beta_sd, rdesign,
      add_intercept, sparse_design
    );
  }
  vector[nobs] base = combine_effects(
    intercept, beta, nobs, neffs, fdesign, sparse, beta_sd, rdesign,
    add_intercept, sparse_design
  );
  // `centre` is the free-intercept flag: an integrated residual is only
  // mean-centred when a free intercept exists to absorb the level it gives
  // up. ARIMA and GP compose additively; each `apply_*` is inert on its
  // own presence flag, so passing both through is safe when only one is set.
  base = apply_arima_residual(
    base, arima_n_obs, arima_present, centre, arima_T, arima_G,
    arima_p, arima_d, arima_q,
    arima_z, arima_pacf, arima_theta, arima_sigma,
    arima_flat_idx
  );
  return apply_gp_term(
    base, gp_present, centre, gp_T, gp_G, gp_M, gp_L, gp_type, gp_nu, gp_d,
    gp_PHI, gp_eta, gp_rho, gp_alpha, gp_flat_idx
  );
}
 
void regression_priors_lp(
  vector beta, vector beta_sd, array[,] real beta_sd_p,
  int fixed, int random,
  int arima_present, int arima_p, int arima_q,
  matrix arima_z, vector arima_pacf, vector arima_theta,
  array[] real arima_sigma, array[,] real arima_sigma_p,
  array[,] real arima_pacf_p
) {
  effect_priors_lp(beta, beta_sd, beta_sd_p, fixed, random);
  if (arima_present) {
    to_vector(arima_z) ~ std_normal();
    // Partial autocorrelations are Uniform(-1, 1) by default via their
    // parameter bounds. A positive prior sd switches to a Normal(mean, sd)
    // truncated to (-1, 1); the truncation constant is fixed by the bounds
    // and so is dropped. A non-positive sd leaves the Uniform default.
    if (arima_p > 0 && arima_pacf_p[2, 1] > 0) {
      arima_pacf ~ normal(arima_pacf_p[1, 1], arima_pacf_p[2, 1]);
    }
    if (arima_q > 0) {
      arima_theta ~ std_normal();
    }
    arima_sigma[1] ~ normal(
      arima_sigma_p[1, 1], arima_sigma_p[2, 1]
    ) T[0, ];
  }
}
 
void gp_priors_lp(
  int gp_present, matrix gp_eta,
  array[] real gp_rho, array[] real gp_alpha,
  array[,] real gp_rho_p, array[,] real gp_alpha_p
) {
  if (gp_present) {
    to_vector(gp_eta) ~ std_normal();
    gp_rho[1] ~ lognormal(gp_rho_p[1, 1], gp_rho_p[2, 1]);
    gp_alpha[1] ~ normal(gp_alpha_p[1, 1], gp_alpha_p[2, 1]) T[0, ];
  }
}