Bivariate normal distribution

[rtrangucci]

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)
• [x] 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

[bgoodri]

My guess is the Pan paper is the best approach for the _lcdf but we would have to have a heuristic for how many terms of the infinite sum to compute based on Phi(a), Phi(b), asin(rho), etc.

[rtrangucci]

Yeah, that sounds good to me. Here's a reference for the gradient of the CDF: https://blogs.sas.com/content/iml/2013/09/20/gradient-of-the-bivariate-normal-cumulative-distribution.html

[rtrangucci]

FYI, I'm waiting on #978 to work on the nonstandardized distributions.