gaussian_process.stan

Go to the documentation of this file.

 
vector diagSPD_EQ(real alpha, real rho, real L, int M) {
  vector[M] indices = linspaced_vector(M, 1, M);
  real factor = alpha * sqrt(sqrt(2 * pi()) * rho);
  real exponent = -0.25 * (rho * pi() / 2 / L)^2;
  return factor * exp(exponent * square(indices));
}
 
vector matern_indices(int M, real L) {
  vector[M] indices = linspaced_vector(M, 1, M);
  return square(pi() / (2 * L) * indices);
}
 
vector diagSPD_Matern12(real alpha, real rho, real L, int M) {
  vector[M] denom = 1 / rho + rho * matern_indices(M, L);
  return alpha * sqrt(2 ./ denom);
}
 
vector diagSPD_Matern32(real alpha, real rho, real L, int M) {
  real factor = 2 * alpha * (sqrt(3) / rho)^1.5;
  vector[M] denom = 3 / square(rho) + matern_indices(M, L);
  return factor ./ denom;
}
 
vector diagSPD_Matern52(real alpha, real rho, real L, int M) {
  real factor = 16 * pow(sqrt(5) / rho, 5);
  vector[M] denom = 3 * pow(5 / square(rho) + matern_indices(M, L), 3);
  return alpha * sqrt(factor ./ denom);
}
 
vector diagSPD_Periodic(real alpha, real rho, int M) {
  real a = inv_square(rho);
  vector[M] indices = linspaced_vector(M, 1, M);
  vector[M] q = exp(
    log(alpha) + 0.5 *
      (log(2) - a + to_vector(log_modified_bessel_first_kind(indices, a)))
  );
  return append_row(q, q);
}
 
vector update_gp(matrix PHI, int M, real L, real alpha,
                 real rho, vector eta, int type, real nu) {
  vector[type == 1 ? 2 * M : M] diagSPD;
  if (type == 0) {
    diagSPD = diagSPD_EQ(alpha, rho, L, M);
  } else if (type == 1) {
    diagSPD = diagSPD_Periodic(alpha, rho, M);
  } else if (type == 2) {
    if (nu == 0.5) {
      diagSPD = diagSPD_Matern12(alpha, rho, L, M);
    } else if (nu == 1.5) {
      diagSPD = diagSPD_Matern32(alpha, rho, L, M);
    } else if (nu == 2.5) {
      diagSPD = diagSPD_Matern52(alpha, rho, L, M);
    } else {
      reject("nu must be one of 0.5, 1.5, or 2.5; found nu=", nu);
    }
  }
  return PHI * (diagSPD .* eta);
}
 
vector apply_gp_term(vector base, int present, int T, int G,
                     int M, real L, int type, real nu, int d,
                     matrix PHI, matrix eta,
                     array[] real rho, array[] real alpha,
                     array[] int flat_idx) {
  if (!present) return base;
  matrix[T, G] gp_eps;
  for (g in 1:G) {
    // Free GP values (length T for d = 0, T - d for d >= 1).
    vector[rows(PHI)] f = update_gp(
      PHI, M, L, alpha[1], rho[1], eta[, g], type, nu
    );
    if (d == 0) {
      gp_eps[, g] = f;
    } else {
      // Anchor the first d entries to zero and place the free values
      // from (d + 1):T, then integrate d times. The leading zeros are
      // preserved by each cumulative_sum pass.
      vector[T] col = rep_vector(0.0, T);
      col[(d + 1):T] = f;
      for (i in 1:d) {
        col = cumulative_sum(col);
      }
      gp_eps[, g] = col;
    }
  }
  return base + to_vector(gp_eps)[flat_idx];
}