Bivariate normal distribution

Description

Add the bivariate normal distribution with scalar inputs. Specifically, add:

• binormal_lpdf(reals y_1, reals y_2 | reals mu_1, reals mu_2, reals sigma_1, reals sigma_2, reals rho)
• binormal_cdf(reals y_1, reals y_2 | reals mu_1, reals mu_2, reals sigma_1, reals sigma_2, reals rho)
• binormal_lcdf(reals y_1, reals y_2 | reals mu_1, reals mu_2, reals sigma_1, reals sigma_2, reals rho)
• binormal_lccdf(reals y_1, reals y_2 | reals mu_1, reals mu_2, reals sigma_1, reals sigma_2, reals rho)
• binormal_rng(reals mu_1, reals mu_2, reals sigma_1, reals sigma_2, reals rho)
• std_binormal_lcdf(reals y_1, reals y_2 | reals rho)
• std_binormal_lpdf(reals y_1, reals y_2 | reals rho)
• std_binormal_cdf(reals y_1, reals y_2 | reals rho)

Reference: #2356

Example

Bivariate probit model

functions {
  int count_cases(int[] y1, int[] y2, int case_1, int case_2) {
    int N = size(y1);
    int N_cases = 0;
    for (n in 1:N) {
      if (y1[n] == case_1 && y2[n] == case_2)
        N_cases += 1;
    }
    return N_cases;
  }
}
data {
  int<lower=1> N;
  int<lower=1> K1;
  int<lower=1> K2;
  int<lower=0,upper=1> y1[N];
  int<lower=0,upper=1> y2[N];
  matrix[N, K1] X1;
  matrix[N, K2] X2;
}
transformed data {
  int N_00 = count_cases(y1, y2, 0, 0);
  int N_10 = count_cases(y1, y2, 1, 0);
  int N_01 = count_cases(y1, y2, 0, 1);
  int N_11 = count_cases(y1, y2, 1, 1);
  int n_00 = 1;
  int n_10 = 1;
  int n_01 = 1;
  int n_11 = 1;
  int idx_00[N_00];
  int idx_10[N_10];
  int idx_01[N_01];
  int idx_11[N_11];
  for (n in 1:N) {
    if(y1[n] == 0 && y2[n] == 0) {
      idx_00[n_00] = n;
      n_00 += 1;
    }
    if(y1[n] == 1 && y2[n] == 0) {
      idx_10[n_10] = n;
      n_10 += 1;
    }
    if(y1[n] == 0 && y2[n] == 1) {
      idx_01[n_01] = n;
      n_01 += 1;
    }
    if(y1[n] == 1 && y2[n] == 1) {
      idx_11[n_11] = n;
      n_11 += 1;
    }
  }
}
parameters {
  real<lower=0,upper=1> rho_raw;
  vector[K1] beta_1;
  vector[K2] beta_2;
  real alpha_1;
  real alpha_2;
}
transformed parameters {
  real rho = -1 + 2 * rho_raw;
}
model {
  vector[N] mu_1 = alpha_1 + X1 * beta_1;
  vector[N] mu_2 = alpha_2 + X2 * beta_2;
  rho_raw ~ beta(4,4);
  alpha_1 ~ normal(0, 2);
  alpha_2 ~ normal(0, 2);
  beta_1 ~ normal(0, 2);
  beta_2 ~ normal(0, 2);
  target += std_binormal_lcdf(-mu_1[idx_00], -mu_2[idx_00] | rho);
  target += std_binormal_lcdf(mu_1[idx_10], -mu_2[idx_10] | -rho);
  target += std_binormal_lcdf(-mu_1[idx_01], mu_2[idx_01] | -rho);
  target += std_binormal_lcdf(mu_1[idx_11], mu_2[idx_11] | rho);
}

Additional Information

binormal_cdf and binormal_lcdf can implement the following algorithms:

• http://www.math.wsu.edu/faculty/genz/papers/bvnt.pdf
• https://arxiv.org/pdf/1004.3616.pdf
• https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2924071
• https://projecteuclid.org/euclid.aos/1176344739

Current Math Version

v2.18.0