hdarjus / sv-comparison Goto Github PK

  for (i in 1:N) {
    if (i == 1) {
      h_i_mean[i] = 0;
      h_i_sigma[i] = pow(1-square(phi), -0.5);
    } else {
      h_i_mean[i] = fma(phi, h[i-1], 0);
      h_i_sigma[i] = 1;
    }
    exph[i] = exp(fma(h[i], 0.5*sigma, 0.5*mu));
  }

data {
  int<lower=0> N;
  real y[N];
}

transformed data {
  real<lower=0> phi_prior_a = 1;
  real<lower=0> phi_prior_b = 1;
  real mu_prior_mu = 0;
  real<lower=0> mu_prior_sigma = 100;
  real<lower=0> sigma_prior_shape = 0.5;
  real<lower=0> sigma_prior_rate = 0.5;
}

parameters {
  vector[N] h;
  real<lower=0,upper=1> phi_beta;
  real mu;
  real<lower=0> sigma2;
}

transformed parameters {
  real<lower=-1,upper=1> phi = fma(phi_beta, 2, -1);
  real<lower=0> sigma = sqrt(sigma2);
  vector[N] h_i_mean;
  vector<lower=0>[N] h_i_sigma;
  h_i_mean[1] = 0;
  h_i_sigma[1] = pow(1 - square(phi), -0.5);
  h_i_sigma[2:N] = rep_vector(1, (N - 1));
  h_i_mean[2:N] = phi * h[1:(N-1)]; // This added 0 before, but idt that's needed?
}

model {
  mu ~ normal(mu_prior_mu, mu_prior_sigma);
  phi_beta ~ beta(phi_prior_a, phi_prior_b);
  sigma2 ~ gamma(sigma_prior_shape, sigma_prior_rate);

  h ~ normal(h_i_mean, h_i_sigma);
  y ~ normal(0, exp(h * (0.5 * sigma) + (0.5 * mu)));
}