From 8f0dea3bc3e39dbe31e62c4c629d51ded74ce4a8 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 9 Mar 2022 16:17:55 +0800 Subject: [PATCH 01/11] First impl reduced form --- .../copula/stan/mixed_copula_example.stan | 74 ++-- functions/copula/mixed_copula.stanfunctions | 381 +++++------------- 2 files changed, 144 insertions(+), 311 deletions(-) diff --git a/examples/copula/stan/mixed_copula_example.stan b/examples/copula/stan/mixed_copula_example.stan index d371a198..7c09136c 100644 --- a/examples/copula/stan/mixed_copula_example.stan +++ b/examples/copula/stan/mixed_copula_example.stan @@ -5,7 +5,7 @@ data { int N; // number of observations array[3] int J; // total number of outcomes int K; // number of covariates for concatenated design matrices (X1, .., XK) - array[N, J[1]] real Yn; // normal outcomes + matrix[N, J[1]] Yn; // normal outcomes array[N, J[2]] int Yb; // bernoulli outcomes array[N, J[3]] int Yp; // poisson outcomes matrix[N, K] X; // concatenated design matrices (X1, ..., XK) @@ -14,52 +14,48 @@ data { } transformed data { int J_all = sum(J); - array[2] matrix[J_all, N] bounds_ind = get_bounds_indicator(Yn, Yb, Yp, N, - J[1], J[2], - J[3], J_all); + // Create separate design matrices for each outcome type + matrix[N, Kj[1]] Xn = X[ : , 1:Kj[1]]; + matrix[N, Kj[2]] Xb = X[ : , (Kj[1] + 1):(Kj[1] + Kj[2])]; + matrix[N, Kj[3]] Xp = X[ : , (Kj[2] + 1):(Kj[2] + Kj[3])]; vector[J_all] mu_zero = rep_vector(0, J_all); } parameters { - vector[K] beta; // long vector of all regression coefficients + vector[Kj[1]] beta_n; // Vector of normal regression coefficients + vector[Kj[2]] beta_b; // Vector of bernoulli regression coefficients + vector[Kj[3]] beta_p; // Vector of poisson regression coefficients cholesky_factor_corr[J_all] L; // Cholesky decomposition of JxJ correlation matrix - array[J[1]] real sigmasq; // Jn-dim array of variances (may be 0-dim) - array[J[2]] vector[N] uraw_b; // latent variables for bernoulli outcomes - array[J[3]] vector[N] uraw_p; // latent variables for Poisson outcomes + vector[J[1]] sigmasq; // Jn-dim vector of variances (may be 0-dim) + matrix[N, J[2]] uraw_b; // latent variables for bernoulli outcomes + matrix[N, J[3]] uraw_p; // latent variables for Poisson outcomes } model { // initialize variables - array[J[1]] real sigma; // stdev - matrix[N, J_all] mu_glm = get_means(beta, X, J, Kj); - - // get standard deviations and declare prior on sigmasq if applicable - if (J[1] > 0) { - sigma = sqrt(sigmasq); - sigmasq ~ inv_gamma(1.0e-4, 1.0e-4); - - for (j in 1 : J[1]) - to_vector(Yn[ : , j]) ~ normal(to_vector(mu_glm[ : , j]), sigma[j]); + vector[J[1]] sigma = sqrt(sigmasq); // stdev + // Calculate the means for each regression + matrix[N, J[1]] mu_n = Xn * rep_matrix(beta_n, J[1]); + matrix[N, J[2]] mu_b = Xb * rep_matrix(beta_b, J[2]); + matrix[N, J[3]] mu_p = Xp * rep_matrix(beta_p, J[3]); + + sigmasq ~ inv_gamma(1.0e-4, 1.0e-4); + // priors for regression coefficients + beta_n ~ normal(0.0, 10.0); + beta_b ~ normal(0.0, 10.0); + beta_p ~ normal(0.0, 10.0); + + for (n in 1 : N) { + // For each individual, construct a matrix containing all needed + // marginal calculations (e.g., latent variables, bounds, and indicators) + matrix[J_all, 5] marginals = append_row( + append_row(normal_marginal(Yn[n], mu_n[n], sigma), + bernoulli_marginal(Yb[n], mu_b[n], uraw_b[n])), + poisson_marginal(Yp[n], mu_p[n], uraw_p[n]) + ); + + // Increment log-likelihoo + marginals ~ mixed_copula_cholesky(mu_zero, L); + Yn[n] ~ normal(mu_n[n], sigma); } - // prior for beta - beta ~ normal(0.0, 10.0); - - // log target increment - if (special == 1) { - target += mixed_cop_sp_lp(Yn, Yb, Yp, - uraw_b, uraw_p, - mu_glm, - sigma, - L, - bounds_ind, - mu_zero, - J, N); - } else target += mixed_cop_lp(Yn, Yb, Yp, - uraw_b, uraw_p, - mu_glm, - sigma, - L, - bounds_ind, - mu_zero, - J, N); } generated quantities { corr_matrix[J_all] Gamma = multiply_lower_tri_self_transpose(L); diff --git a/functions/copula/mixed_copula.stanfunctions b/functions/copula/mixed_copula.stanfunctions index ad0e27ee..a73f44ee 100644 --- a/functions/copula/mixed_copula.stanfunctions +++ b/functions/copula/mixed_copula.stanfunctions @@ -1,8 +1,8 @@ /** @addtogroup mixed_copula Mixed Discrete-Continuous Gaussian Copula Functions * - * The mixed-discrete normal copula. The method is from - * Smith, Michael & Khaled, Mohamad. (2011). Estimation of Copula Models With Discrete Margins via Bayesian Data Augmentation. - * Journal of the American Statistical Association. 107. 10.2139/ssrn.1937983. + * The mixed-discrete normal copula. The method is from + * Smith, Michael & Khaled, Mohamad. (2011). Estimation of Copula Models With Discrete Margins via Bayesian Data Augmentation. + * Journal of the American Statistical Association. 107. 10.2139/ssrn.1937983. * * @include \copula\mixed_copula.stanfunctions * @@ -25,25 +25,33 @@ * @param J Integer * @return Vector in [0, 1] of continuous outcome */ -vector normal_marginal(array[] real y, row_vector mu_glm, array[] real sigma, int J) { - vector[J] u; - + + +matrix normal_marginal(row_vector y, row_vector mu_glm, vector sigma) { + int J = num_elements(y); + matrix[J, 5] rtn; + real lb = negative_infinity(); + real ub = positive_infinity(); + real lb_ind = 0; + real ub_ind = 0; + for (j in 1 : J) { - u[j] = Phi_approx((y[j] - mu_glm[j]) / sigma[j]); + real u = Phi_approx((y[j] - mu_glm[j]) / sigma[j]); + rtn[j] = [u, lb, ub, lb_ind, ub_ind]; } - - return u; + + return rtn; } /** * Bernoulli marginal * * Bernoulli marginal for mixed continuous-discrete Gaussian - * copula. + * copula. * * The lb and ub: * if y == 0, upper bound at inv_Phi( y, mu ) - * if y == 1, lower bound at inv_Phi( y - 1, mu ) + * if y == 1, lower bound at inv_Phi( y - 1, mu ) * * @copyright Ethan Alt, Sean Pinkney 2021 \n * @@ -52,124 +60,114 @@ vector normal_marginal(array[] real y, row_vector mu_glm, array[] real sigma, in * @param J Integer * @return Vector in [0, 1] of continuous outcome */ -array[] vector bernoulli_marginal(array[] int y, row_vector mu_glm, int J) { - array[2] vector[J] b; - +matrix bernoulli_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { + int J = num_elements(y); + matrix[J, 5] rtn; + row_vector[J] mu_glm_logit = inv_logit(mu_glm); + +// u, lb, ub, lb_ind, ub_ind + for (j in 1 : J) { + real u = u_raw[j]; + real lb; + real ub; + real lb_ind = 0; + real ub_ind = 0; if (y[j] == 0) { - b[1, j] = negative_infinity(); - b[2, j] = inv_Phi(bernoulli_cdf(0 | mu_glm[j])); + lb = negative_infinity(); + ub = inv_Phi(bernoulli_cdf(0 | mu_glm_logit[j])); + ub_ind = 1; } else { - b[1, j] = inv_Phi(bernoulli_cdf(0 | mu_glm[j])); - b[2, j] = positive_infinity(); + lb = inv_Phi(bernoulli_cdf(0 | mu_glm_logit[j])); + ub = positive_infinity(); + lb_ind = 1; } + rtn[j] = [u, lb, ub, lb_ind, ub_ind]; } - - return b; + + return rtn; } /** - * Poisson marginal + * Binomial marginal * - * Poisson marginal for mixed continuous-discrete Gaussian - * copula. + * Binomial marginal for mixed continuous-discrete Gaussian + * copula. * * The lb and ub: * Always upper bound at inv_Phi( y, mu ) - * If y != 0, lower bound at inv_Phi(y-1, mu) + * If y != 0, lower bound at inv_Phi(y-1, mu) * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param y int[] Array of integers + * @param n int[] Array of numerator integers + * @param d int[] Array of denominator integers * @param mu_glm Row vector * @param J Integer * @return Vector in [0, 1] of continuous outcome */ -array[] vector poisson_marginal(array[] int y, row_vector mu_glm, int J) { - array[2] vector[J] b; - + +// u, lb, ub, lb_ind, ub_ind +matrix binomial_marginal(array[] int n, array[] int d, row_vector mu_glm, row_vector u_raw) { + int J = num_elements(n); + matrix[J, 5] rtn; + row_vector[J] mu_glm_logit = inv_logit(mu_glm); + for (j in 1 : J) { - b[2, j] = inv_Phi(poisson_cdf(y[j] | mu_glm[j])); - if (y[j] > 0) { - b[1, j] = inv_Phi(poisson_cdf(y[j] - 1 | mu_glm[j])); + real u = u_raw[j]; + real lb = negative_infinity(); + real ub = inv_Phi(binomial_cdf(n[j] | d[j], mu_glm_logit[j])); + real lb_ind = 0; + real ub_ind = 1; + if (n[j] > 0) { + lb_ind = 1; + lb = inv_Phi(binomial_cdf(n[j] - 1 | d[j], mu_glm_logit[j])); } + rtn[j] = [u, lb, ub, lb_ind, ub_ind]; } - - return b; + + return rtn; } /** -* Mixed Copula Special Log-Probability Function -* -* This is potentially faster when there are many margins. However, -* it is slower with few margins. -* -* @copyright Ethan Alt, Sean Pinkney 2021 \n -* -* @param Yn real[,] 2d-Array of Normal depended variates -* @param Yb int[,] 2d-Array of Bernoulli dependent variates -* @param Yp int[,] 2d-Array of Poisson dependent variates -* @param uraw_b Vector of standard uniform variates for Bernoulli margins -* @param uraw_p Vector of standard uniform variates for Poisson margins -* @param mu_glm Matrix of GLM means -* @param sigma Vector of standard deviations for Gaussian variables -* @param L Cholesky Factor Correlation -* @param bounds_ind Matrix giving whether there is an lower([1]) or upper([2]) bound -* @param tmvn_mu Vector of means for the Multivariate (Truncated) Gaussian, almost always 0 -* @param J Integer array specifying the number of margins for each distribution -* @param N Integer number of observations -* @return Real log-probability -*/ -real mixed_cop_sp_lp(array[,] real Yn, array[,] int Yb, array[,] int Yp, - array[] vector uraw_b, array[] vector uraw_p, matrix mu_glm, - array[] real sigma, matrix L, array[] matrix bounds_ind, - vector tmvn_mu, array[] int J, int N) { - real lp = 0; - int J_all = sum(J); - array[2] vector[J[2]] bb; - array[2] vector[J[2]] bp; - - vector[J_all] u; // latent variables - vector[J_all] lb; - vector[J_all] ub; - - // loop through independent observations - for (i in 1 : N) { - // get indicators of upper/lower bound - vector[J_all] lb_ind = col(bounds_ind[1], i); - vector[J_all] ub_ind = col(bounds_ind[2], i); - - int pos = 1; - - if (J[1] > 0) { - u[pos : J[1]] = normal_marginal(Yn[i], head(mu_glm[i], J[1]), sigma, J[1]); - lb[pos : J[1]] = rep_vector(negative_infinity(), J[1]); - ub[pos : J[1]] = rep_vector(positive_infinity(), J[1]); - pos += J[1]; - } - - if (J[2] > 0) { - bb = bernoulli_marginal(Yb[i], segment(mu_glm[i], pos, J[2]), J[2]); - u[pos : J[1] + J[2]] = to_vector(uraw_b[1 : J[2], i]); - lb[pos : J[1] + J[2]] = bb[1]; - ub[pos : J[1] + J[2]] = bb[2]; - - pos += J[2]; - } - - if (J[3] > 0) { - bp = poisson_marginal(Yp[i], segment(mu_glm[i], pos, J[3]), J[3]); - - u[pos : J_all] = to_vector(uraw_p[1 : J[3], i]); - lb[pos : J_all] = bp[1]; - ub[pos : J_all] = bp[2]; + * Poisson marginal + * + * Poisson marginal for mixed continuous-discrete Gaussian + * copula. + * + * The lb and ub: + * Always upper bound at inv_Phi( y, mu ) + * If y != 0, lower bound at inv_Phi(y-1, mu) + * + * @copyright Ethan Alt, Sean Pinkney 2021 \n + * + * @param y int[] Array of integers + * @param mu_glm Row vector + * @param J Integer + * @return Vector in [0, 1] of continuous outcome + */ +matrix poisson_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { + int J = num_elements(y); + matrix[J, 5] rtn; + row_vector[J] mu_glm_exp = exp(mu_glm); + + for (j in 1 : J) { + real u = u_raw[j]; + real lb = negative_infinity(); + real ub = inv_Phi(poisson_cdf(y[j] | mu_glm_exp[j])); + real lb_ind = 0; + real ub_ind = 1; + if (y[j] > 0) { + lb_ind = 1; + lb = inv_Phi(poisson_cdf(y[j] - 1 | mu_glm_exp[j])); } - - lp += multi_normal_cholesky_truncated_lpdf(u | tmvn_mu, L, lb, ub, lb_ind, ub_ind); + rtn[j] = [u, lb, ub, lb_ind, ub_ind]; } - return lp; + + return rtn; } + /** * Mixed Copula Log-Probability Function * @@ -178,187 +176,26 @@ real mixed_cop_sp_lp(array[,] real Yn, array[,] int Yb, array[,] int Yp, * @param Yn real[,] 2d-Array of Normal depended variates * @param Yb int[,] 2d-Array of Bernoulli dependent variates * @param Yp int[,] 2d-Array of Poisson dependent variates - * @param uraw_b Vector of standard uniform variates for Bernoulli margins + * @param uraw_be Vector of standard uniform variates for Bernoulli margins * @param uraw_p Vector of standard uniform variates for Poisson margins * @param mu_glm Matrix of GLM means * @param sigma Vector of standard deviations for Gaussian variables - * @param L Cholesky Factor Correlation + * @param L Cholesky Factor Correlation * @param bounds_ind Matrix giving whether there is an lower([1]) or upper([2]) bound * @param tmvn_mu Vector of means for the Multivariate (Truncated) Gaussian, almost always 0 * @param J Integer array specifying the number of margins for each distribution * @param N Integer number of observations - * @return Real log-probability + * @return Real log-probability */ -real mixed_cop_lp(array[,] real Yn, array[,] int Yb, array[,] int Yp, - array[] vector uraw_b, array[] vector uraw_p, matrix mu_glm, - array[] real sigma, matrix L, array[] matrix bounds_ind, - vector tmvn_mu, array[] int J, int N) { - real lp = 0; - int J_all = sum(J); - - // loop through independent observations - for (i in 1 : N) { - vector[J_all] u; // latent variables - vector[J_all] lb = rep_vector(negative_infinity(), J_all); // initialize lower bounds to -inf - vector[J_all] ub = rep_vector(positive_infinity(), J_all); // initialize upper bounds to +inf - - // get indicators of upper/lower bound - vector[J_all] lb_ind = col(bounds_ind[1], i); - vector[J_all] ub_ind = col(bounds_ind[2], i); - - int j = 0; // counter for 1:J - int jn = 0; // counter for 1:Jn - int jb = 0; // counter for 1:Jb - int jp = 0; // counter for 1:Jp - - // for normal; increment by log likelihood and obtain u|y deterministically - while (jn < J[1]) { - real z_ij; - jn += 1; - j += 1; - z_ij = (Yn[i, jn] - mu_glm[i, j]) / sigma[jn]; - u[j] = Phi_approx(z_ij); - } - - // for Bernoulli, u|y is random and must update lb and ub: - // if y == 0, upper bound at inv_Phi( y, mu ) - // if y == 1, lower bound at inv_Phi(y-1, mu) (y-1 == 0). - while (jb < J[2]) { - jb += 1; - j += 1; - u[j] = uraw_b[jb, i]; - if (Yb[i, jb] == 0) { - ub[j] = inv_Phi(bernoulli_cdf(0 | mu_glm[i, j])); - } else { - lb[j] = inv_Phi(bernoulli_cdf(0 | mu_glm[i, j])); - } - } - - // for Poisson, u|y is random and must update lb and ub: - // Always upper bound at inv_Phi( y, mu ) - // If y != 0, lower bound at inv_Phi(y-1, mu) - while (jp < J[3]) { - int y_ij; - real mu_ij; - jp += 1; - j += 1; - u[j] = uraw_p[jp, i]; - y_ij = Yp[i, jp]; - mu_ij = mu_glm[i, j]; - ub[j] = inv_Phi(poisson_cdf(y_ij | mu_ij)); - if (y_ij > 0) { - lb[j] = inv_Phi(poisson_cdf(y_ij - 1 | mu_ij)); - } - } - lp += multi_normal_cholesky_truncated_lpdf(u | tmvn_mu, L, lb, ub, lb_ind, ub_ind); - } // end i loop - return lp; -} +real mixed_copula_cholesky_lpdf(matrix marginals, vector tmvn_mu, matrix L) { + int J = num_elements(tmvn_mu); + vector[J] u = marginals[ : , 1]; + vector[J] lb = marginals[ : , 2]; + vector[J] ub = marginals[ : , 3]; + vector[J] lb_ind = marginals[ : , 4]; + vector[J] ub_ind = marginals[ : , 5]; -//----------------------------------------------------------------- -// Get GLM means -// -// * @param beta long vector of regression coefs for all outcomes -// * @param X concatenated design matrices -// * @param Jn number of normal outcomes -// * @param Jb number of bernoulli outcomes -// * @param Jp number of poisson outcomes -// * @param J total number of outcomes -// * @param N number of observations -// * @param Kj J-dim integer vector giving number of covariates per variable -// * @param Xindx Xindx Jx2 matrix giivng start and end indices of X for design matrix -// -// * @return NxJ matrix of GLM means -//---------------------------------------------------------------- - -/** - * Get GLM means - * - * @copyright Ethan Alt, Sean Pinkney 2021 \n - * - * @param beta Vector of regression coefficients - * @param X Matrix concatenated design matrix across all margins - * @param J Integer array specifying the number of margins for each distribution - * @param Kj Integer array of the number of covariates for each margin - * @return Matrix mean outcome matrix from marginal GLMs - */ -matrix get_means(vector beta, matrix X, array[] int J, array[] int Kj) { - int N = rows(X); - matrix[N, num_elements(J)] mu; - int pos = 1; - int outcome = 1; - - if (J[1] > 0) { - mu[ : , outcome] = block(X, 1, 1, N, Kj[1]) * segment(beta, 1, Kj[1]); - pos += Kj[1]; - outcome += 1; - } - - if (J[2] > 0) { - mu[ : , outcome] = inv_logit(block(X, 1, pos, N, Kj[2]) * segment(beta, pos, Kj[2])); - pos += Kj[2]; - outcome += 1; - } - - if (J[3] > 0) { - mu[ : , outcome] = exp(block(X, 1, pos, N, Kj[3]) * segment(beta, pos, Kj[3])); - } - - return mu; + return multi_normal_cholesky_truncated_lpdf(u | tmvn_mu, L, lb, ub, lb_ind, ub_ind); } -/** - * Get Bound Indicators - * - * @copyright Ethan Alt 2021 \n - * - * @param Yn Real 2d-array of normal outcomes - * @param Yb Integer 2d-array of Bernoulli outcomes - * @param Yp Integer 2d-array of Poisson outcomes - * @param N Integer number of observations - * @param Jn number of normal outcomes - * @param Jb number of bernoulli outcomes - * @param Jp number of poisson outcomes - * @param J total number of outcomes - * @return N x J x 2 array indicating whether there is a lower/upper bound. - */ -array[] matrix get_bounds_indicator(array[,] real Yn, array[,] int Yb, array[,] int Yp, - int N, int Jn, int Jb, int Jp, int J) { - matrix[J, N] lb_ind = rep_matrix(0, J, N); // initialize to no lower bound - matrix[J, N] ub_ind = rep_matrix(0, J, N); // initialize to no upper bound - array[2] matrix[J, N] out; // output is a 2-dim array of JxN matrices - int j = Jn; // counter over J (normal outcomes are unbounded) - int jb = 0; // counter over Jb - int jp = 0; // counter over Jp - // loop over binomial variables (if any) - // upper bound if y == 0 - // lower bound if y == 1 - while (jb < Jb) { - jb += 1; - j += 1; - for (i in 1 : N) { - if (Yb[i, jb] == 0) { - ub_ind[j, i] = 1; - } else { - lb_ind[j, i] = 1; - } - } - } - // loop over poisson variables (if any) - // always upper bound - // lower bound for y > 0 - while (jp < Jp) { - jp += 1; - j += 1; - for (i in 1 : N) { - ub_ind[j, i] = 1; - if (Yp[i, jp] > 0) { - lb_ind[j, i] = 1; - } - } - } - out[1] = lb_ind; - out[2] = ub_ind; - return out; -} /** @} */ \ No newline at end of file From 38ea1888e36ebd3a6c72048e0dcf5fe7c0fa3797 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 9 Mar 2022 17:00:52 +0800 Subject: [PATCH 02/11] Consistent naming --- examples/copula/R/mixed_copula_example.R | 30 +++---- .../copula/stan/mixed_copula_example.stan | 6 +- ...unctions => gaussian_copula.stanfunctions} | 78 ++++++++----------- 3 files changed, 50 insertions(+), 64 deletions(-) rename functions/copula/{mixed_copula.stanfunctions => gaussian_copula.stanfunctions} (68%) diff --git a/examples/copula/R/mixed_copula_example.R b/examples/copula/R/mixed_copula_example.R index 4d7986fe..6fe247a9 100644 --- a/examples/copula/R/mixed_copula_example.R +++ b/examples/copula/R/mixed_copula_example.R @@ -1,6 +1,6 @@ library(cmdstanr) # setwd("~\helpful_stan_functions") -fp <- file.path("./examples/copula/stan/mixed_copula_example.stan") +fp <- file.path("./examples/copula/stan/gaussian_copula_example.stan") mod <- cmdstan_model(fp, include_paths = c("./functions/copula", "./functions/distribution")) @@ -37,27 +37,27 @@ family.list = list( ) #' Joint analysis of multivariate GLMs of mixed types -#' +#' #' Uses rstan to estimate a multivariate GLM with correlated #' responses via Gaussian copula -#' +#' #' @param formula.list a list of formulas #' @param family.list a list of families giving the density of each formula. Canonical link functions always used. #' @param data data.frame giving all variables in `formula.list` #' @param ... arguments to pass onto `rstan::sampling` -#' +#' #' @return rstan object fitting the posterior get_data = function(formula.list, family.list, data) { - + N = nrow(data) J = length(formula.list) - + ## reorder variables if necessary n.indx = sapply(family.list, function(x) x$family) == 'gaussian' b.indx = sapply(family.list, function(x) x$family) == 'binomial' p.indx = sapply(family.list, function(x) x$family) == 'poisson' new.order = c(which(n.indx), which(b.indx), which(p.indx)) - + ## gtotal number of endpoints and matrix of each outcome type Jn = sum(n.indx) Jb = sum(b.indx) @@ -66,26 +66,26 @@ get_data = function(formula.list, family.list, data) { Yn = Y[, n.indx, drop = F] Yb = Y[, b.indx, drop = F] Yp = Y[, p.indx, drop = F] - + ## store character vector giving dist name and number so results make sense distnames = c( paste0('gaussian[', 1:Jn, ']'), paste0('binomial[', 1:Jb, ']'), paste0('poisson[', 1:Jp, ']') ) - - + + if ( !all( new.order == 1:J ) ) { message("Rearranging order of dependent variables to normal, binomial, poisson") formula.list = formula.list[new.order] family.list = family.list[new.order] } - + ## get design matrices and number of covariates per regression Xlist = lapply(formula.list, model.matrix) Xbig = do.call(cbind, Xlist) K = ncol(Xbig) Kj = sapply(Xlist, ncol) - + ## obtain matrix giving start and end values of x matrix per endpoint Xindx = matrix(nrow = J, ncol = 2) xstart = 1 @@ -94,7 +94,7 @@ get_data = function(formula.list, family.list, data) { Xindx[j, ] = c(xstart, xend) xstart = xend + 1 } - + ## construct stan data standat = list( 'N' = N, @@ -110,7 +110,7 @@ get_data = function(formula.list, family.list, data) { 'Xindx' = Xindx, 'Kj' = Kj ) - + ## fit stan model # fit = rstan::sampling(stanmod, data = standat, ...) # fit = stanmod$sample(data = standat, @@ -124,7 +124,7 @@ get_data = function(formula.list, family.list, data) { # beta.new = c(beta.new, paste0(distnames[j], '[', seq(0,Kj[j] - 1), ']') ) # } # names(fit)[beta.old.indx] = beta.new - # + # ## return stan object return(standat) } diff --git a/examples/copula/stan/mixed_copula_example.stan b/examples/copula/stan/mixed_copula_example.stan index 7c09136c..79ec8f03 100644 --- a/examples/copula/stan/mixed_copula_example.stan +++ b/examples/copula/stan/mixed_copula_example.stan @@ -1,5 +1,5 @@ functions { - #include mixed_copula.stanfunctions + #include gaussian_copula.stanfunctions } data { int N; // number of observations @@ -52,8 +52,8 @@ model { poisson_marginal(Yp[n], mu_p[n], uraw_p[n]) ); - // Increment log-likelihoo - marginals ~ mixed_copula_cholesky(mu_zero, L); + // Increment log-likelihoods + marginals ~ gaussian_copula_cholesky(mu_zero, L); Yn[n] ~ normal(mu_n[n], sigma); } } diff --git a/functions/copula/mixed_copula.stanfunctions b/functions/copula/gaussian_copula.stanfunctions similarity index 68% rename from functions/copula/mixed_copula.stanfunctions rename to functions/copula/gaussian_copula.stanfunctions index a73f44ee..b699222e 100644 --- a/functions/copula/mixed_copula.stanfunctions +++ b/functions/copula/gaussian_copula.stanfunctions @@ -1,10 +1,10 @@ -/** @addtogroup mixed_copula Mixed Discrete-Continuous Gaussian Copula Functions +/** @addtogroup gaussian_copula Mixed Discrete-Continuous Gaussian Copula Functions * * The mixed-discrete normal copula. The method is from * Smith, Michael & Khaled, Mohamad. (2011). Estimation of Copula Models With Discrete Margins via Bayesian Data Augmentation. * Journal of the American Statistical Association. 107. 10.2139/ssrn.1937983. * - * @include \copula\mixed_copula.stanfunctions + * @include \copula\gaussian_copula.stanfunctions * * \ingroup copula * @{ */ @@ -19,26 +19,20 @@ * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param y Real[] Array - * @param mu_glm Row vector - * @param sigma Real[] Array - * @param J Integer - * @return Vector in [0, 1] of continuous outcome + * @param y row_vector of normal outcomes + * @param mu_glm row_vector of regression means + * @param sigma vector of outcome SD's + * @return Matrix containing continuous outcome in [0, 1], + * bounds values and indicators */ - - matrix normal_marginal(row_vector y, row_vector mu_glm, vector sigma) { int J = num_elements(y); matrix[J, 5] rtn; - real lb = negative_infinity(); - real ub = positive_infinity(); - real lb_ind = 0; - real ub_ind = 0; - - for (j in 1 : J) { - real u = Phi_approx((y[j] - mu_glm[j]) / sigma[j]); - rtn[j] = [u, lb, ub, lb_ind, ub_ind]; - } + rtn[ : , 1] = Phi_approx(transpose(y - mu_glm) ./ sigma); + rtn[ : , 2] = rep_vector(negative_infinity(), J); + rtn[ : , 3] = rep_vector(positive_infinity(), J); + rtn[ : , 4] = rep_vector(0, J); + rtn[ : , 5] = rep_vector(0, J); return rtn; } @@ -55,18 +49,17 @@ matrix normal_marginal(row_vector y, row_vector mu_glm, vector sigma) { * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param y int[] Array of integers - * @param mu_glm Row vector - * @param J Integer - * @return Vector in [0, 1] of continuous outcome + * @param y int[] Array of binary outcomes + * @param mu_glm row_vector of regression means + * @param u_raw row_vector of nuisance latent variables + * @return Matrix containing continuous outcome in [0, 1], + * bounds values and indicators */ matrix bernoulli_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { int J = num_elements(y); matrix[J, 5] rtn; row_vector[J] mu_glm_logit = inv_logit(mu_glm); -// u, lb, ub, lb_ind, ub_ind - for (j in 1 : J) { real u = u_raw[j]; real lb; @@ -96,18 +89,17 @@ matrix bernoulli_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { * * The lb and ub: * Always upper bound at inv_Phi( y, mu ) - * If y != 0, lower bound at inv_Phi(y-1, mu) + * If n != 0, lower bound at inv_Phi(n-1, mu) * * @copyright Ethan Alt, Sean Pinkney 2021 \n * * @param n int[] Array of numerator integers * @param d int[] Array of denominator integers - * @param mu_glm Row vector - * @param J Integer - * @return Vector in [0, 1] of continuous outcome + * @param mu_glm row_vector of regression means + * @param u_raw row_vector of nuisance latent variables + * @return Matrix containing continuous outcome in [0, 1], + * bounds values and indicators */ - -// u, lb, ub, lb_ind, ub_ind matrix binomial_marginal(array[] int n, array[] int d, row_vector mu_glm, row_vector u_raw) { int J = num_elements(n); matrix[J, 5] rtn; @@ -141,10 +133,11 @@ matrix binomial_marginal(array[] int n, array[] int d, row_vector mu_glm, row_ve * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param y int[] Array of integers - * @param mu_glm Row vector - * @param J Integer - * @return Vector in [0, 1] of continuous outcome + * @param y int[] Array of integer counts + * @param mu_glm row_vector of regression means + * @param u_raw row_vector of nuisance latent variables + * @return Matrix containing continuous outcome in [0, 1], + * bounds values and indicators */ matrix poisson_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { int J = num_elements(y); @@ -173,21 +166,14 @@ matrix poisson_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param Yn real[,] 2d-Array of Normal depended variates - * @param Yb int[,] 2d-Array of Bernoulli dependent variates - * @param Yp int[,] 2d-Array of Poisson dependent variates - * @param uraw_be Vector of standard uniform variates for Bernoulli margins - * @param uraw_p Vector of standard uniform variates for Poisson margins - * @param mu_glm Matrix of GLM means - * @param sigma Vector of standard deviations for Gaussian variables - * @param L Cholesky Factor Correlation - * @param bounds_ind Matrix giving whether there is an lower([1]) or upper([2]) bound + * + * @param marginals Matrix of marginal calculations, constructed by concatenating the returns + * from the _marginal functions * @param tmvn_mu Vector of means for the Multivariate (Truncated) Gaussian, almost always 0 - * @param J Integer array specifying the number of margins for each distribution - * @param N Integer number of observations + * @param L Cholesky Factor Correlation * @return Real log-probability */ -real mixed_copula_cholesky_lpdf(matrix marginals, vector tmvn_mu, matrix L) { +real gaussian_copula_cholesky_lpdf(matrix marginals, vector tmvn_mu, matrix L) { int J = num_elements(tmvn_mu); vector[J] u = marginals[ : , 1]; vector[J] lb = marginals[ : , 2]; From 1ec852d11a57eacf3f0870d7aba4fed28e065f04 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 9 Mar 2022 17:04:17 +0800 Subject: [PATCH 03/11] Example renaming --- .../R/{mixed_copula_example.R => gaussian_copula_example.R} | 0 .../{mixed_copula_example.stan => gaussian_copula_example.stan} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename examples/copula/R/{mixed_copula_example.R => gaussian_copula_example.R} (100%) rename examples/copula/stan/{mixed_copula_example.stan => gaussian_copula_example.stan} (100%) diff --git a/examples/copula/R/mixed_copula_example.R b/examples/copula/R/gaussian_copula_example.R similarity index 100% rename from examples/copula/R/mixed_copula_example.R rename to examples/copula/R/gaussian_copula_example.R diff --git a/examples/copula/stan/mixed_copula_example.stan b/examples/copula/stan/gaussian_copula_example.stan similarity index 100% rename from examples/copula/stan/mixed_copula_example.stan rename to examples/copula/stan/gaussian_copula_example.stan From 5884950fe0c5257ae37db88c09a1dd5019dc9f03 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 10 Mar 2022 02:07:15 +0800 Subject: [PATCH 04/11] Refactor to use alternate likelihood --- ...e.R => centered_gaussian_copula_example.R} | 10 +- ... => centered_gaussian_copula_example.stan} | 19 +-- ...=> centered_gaussian_copula.stanfunctions} | 147 +++++++++--------- functions/copula/normal_copula.stanfunctions | 34 ++-- 4 files changed, 103 insertions(+), 107 deletions(-) rename examples/copula/R/{gaussian_copula_example.R => centered_gaussian_copula_example.R} (91%) rename examples/copula/stan/{gaussian_copula_example.stan => centered_gaussian_copula_example.stan} (77%) rename functions/copula/{gaussian_copula.stanfunctions => centered_gaussian_copula.stanfunctions} (53%) diff --git a/examples/copula/R/gaussian_copula_example.R b/examples/copula/R/centered_gaussian_copula_example.R similarity index 91% rename from examples/copula/R/gaussian_copula_example.R rename to examples/copula/R/centered_gaussian_copula_example.R index 6fe247a9..633e8745 100644 --- a/examples/copula/R/gaussian_copula_example.R +++ b/examples/copula/R/centered_gaussian_copula_example.R @@ -1,8 +1,10 @@ library(cmdstanr) # setwd("~\helpful_stan_functions") -fp <- file.path("./examples/copula/stan/gaussian_copula_example.stan") -mod <- cmdstan_model(fp, include_paths = c("./functions/copula", - "./functions/distribution")) +fp <- file.path("./examples/copula/stan/centered_gaussian_copula_example.stan") +mod <- cmdstan_model(fp, include_paths = c("./functions", + "./functions/copula", + "./functions/distribution", + "../functions")) ## obtain non-extreme random correlation matrix Gamma = matrix(1, 3, 3) @@ -140,7 +142,7 @@ stan_data[["special"]] <- 1 # stan_data[["special"]] <- 0 mod_out <- mod$sample(data = stan_data, seed = 908453, - # init = 0 , + init = 0 , chains = 2, parallel_chains = 2, iter_warmup = 400, diff --git a/examples/copula/stan/gaussian_copula_example.stan b/examples/copula/stan/centered_gaussian_copula_example.stan similarity index 77% rename from examples/copula/stan/gaussian_copula_example.stan rename to examples/copula/stan/centered_gaussian_copula_example.stan index 79ec8f03..16fdccc3 100644 --- a/examples/copula/stan/gaussian_copula_example.stan +++ b/examples/copula/stan/centered_gaussian_copula_example.stan @@ -1,5 +1,5 @@ functions { - #include gaussian_copula.stanfunctions + #include centered_gaussian_copula.stanfunctions } data { int N; // number of observations @@ -42,20 +42,9 @@ model { beta_n ~ normal(0.0, 10.0); beta_b ~ normal(0.0, 10.0); beta_p ~ normal(0.0, 10.0); - - for (n in 1 : N) { - // For each individual, construct a matrix containing all needed - // marginal calculations (e.g., latent variables, bounds, and indicators) - matrix[J_all, 5] marginals = append_row( - append_row(normal_marginal(Yn[n], mu_n[n], sigma), - bernoulli_marginal(Yb[n], mu_b[n], uraw_b[n])), - poisson_marginal(Yp[n], mu_p[n], uraw_p[n]) - ); - - // Increment log-likelihoods - marginals ~ gaussian_copula_cholesky(mu_zero, L); - Yn[n] ~ normal(mu_n[n], sigma); - } + { normal_marginal(Yn, mu_n, sigma), + bernoulli_marginal(Yb, mu_b, uraw_b), + poisson_marginal(Yp, mu_p, uraw_p) } ~ centered_gaussian_copula_cholesky(L); } generated quantities { corr_matrix[J_all] Gamma = multiply_lower_tri_self_transpose(L); diff --git a/functions/copula/gaussian_copula.stanfunctions b/functions/copula/centered_gaussian_copula.stanfunctions similarity index 53% rename from functions/copula/gaussian_copula.stanfunctions rename to functions/copula/centered_gaussian_copula.stanfunctions index b699222e..d9e90274 100644 --- a/functions/copula/gaussian_copula.stanfunctions +++ b/functions/copula/centered_gaussian_copula.stanfunctions @@ -9,7 +9,7 @@ * \ingroup copula * @{ */ -#include distribution/multi_normal_cholesky_truncated.stanfunctions +#include copula/normal_copula.stanfunctions /** * Normal marginal @@ -25,16 +25,15 @@ * @return Matrix containing continuous outcome in [0, 1], * bounds values and indicators */ -matrix normal_marginal(row_vector y, row_vector mu_glm, vector sigma) { - int J = num_elements(y); - matrix[J, 5] rtn; - rtn[ : , 1] = Phi_approx(transpose(y - mu_glm) ./ sigma); - rtn[ : , 2] = rep_vector(negative_infinity(), J); - rtn[ : , 3] = rep_vector(positive_infinity(), J); - rtn[ : , 4] = rep_vector(0, J); - rtn[ : , 5] = rep_vector(0, J); +array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) { + int N = rows(mu_glm); + int J = cols(mu_glm); + matrix[N, J] adj = rep_matrix(0, N, J); + for (j in 1:J) { + adj[1, j] = normal_lpdf(y[ : , j] | mu_glm[ : , j], sigma[j]); + } - return rtn; + return { Phi_approx((y - mu_glm) ./ rep_matrix(transpose(sigma), N)), adj }; } /** @@ -55,30 +54,18 @@ matrix normal_marginal(row_vector y, row_vector mu_glm, vector sigma) { * @return Matrix containing continuous outcome in [0, 1], * bounds values and indicators */ -matrix bernoulli_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { - int J = num_elements(y); - matrix[J, 5] rtn; - row_vector[J] mu_glm_logit = inv_logit(mu_glm); - - for (j in 1 : J) { - real u = u_raw[j]; - real lb; - real ub; - real lb_ind = 0; - real ub_ind = 0; - if (y[j] == 0) { - lb = negative_infinity(); - ub = inv_Phi(bernoulli_cdf(0 | mu_glm_logit[j])); - ub_ind = 1; - } else { - lb = inv_Phi(bernoulli_cdf(0 | mu_glm_logit[j])); - ub = positive_infinity(); - lb_ind = 1; - } - rtn[j] = [u, lb, ub, lb_ind, ub_ind]; - } - - return rtn; +array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { + int N = rows(mu_glm); + int J = cols(mu_glm); + matrix[N, J] mat_y = to_matrix(y); + matrix[N, J] mat_inv_y = fabs(mat_y - 1); + + matrix[N, J] mu_glm_logit = inv_logit(mu_glm); + matrix[N, J] Lbound = mat_y .* (1 - mu_glm_logit); + matrix[N, J] Ubound = 1 - mat_inv_y .* mu_glm_logit; + matrix[N, J] UmL = Ubound - Lbound; + + return { Lbound + UmL .* u_raw, log(UmL) }; } /** @@ -100,22 +87,24 @@ matrix bernoulli_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { * @return Matrix containing continuous outcome in [0, 1], * bounds values and indicators */ -matrix binomial_marginal(array[] int n, array[] int d, row_vector mu_glm, row_vector u_raw) { - int J = num_elements(n); - matrix[J, 5] rtn; - row_vector[J] mu_glm_logit = inv_logit(mu_glm); - - for (j in 1 : J) { - real u = u_raw[j]; - real lb = negative_infinity(); - real ub = inv_Phi(binomial_cdf(n[j] | d[j], mu_glm_logit[j])); - real lb_ind = 0; - real ub_ind = 1; - if (n[j] > 0) { - lb_ind = 1; - lb = inv_Phi(binomial_cdf(n[j] - 1 | d[j], mu_glm_logit[j])); +array[] matrix binomial_marginal(array[,] int num, array[,] int den, + matrix mu_glm, matrix u_raw) { + int N = rows(mu_glm); + int J = cols(mu_glm); + matrix[N, J] mu_glm_logit = inv_logit(mu_glm); + array[2] matrix[N, J] rtn; + + for (n in 1:N) { + for (j in 1:J) { + real Ubound = binomial_cdf(num[n, j] | den[n, j], mu_glm_logit[n, j]); + real Lbound = 0; + if (num[n, j] > 0) { + Lbound = binomial_cdf(num[n, j] - 1 | den[n, j], mu_glm_logit[n, j]); + } + real UmL = Ubound - Lbound; + rtn[1][n, j] = Lbound + UmL * u_raw[n, j]; + rtn[2][n, j] = log(UmL); } - rtn[j] = [u, lb, ub, lb_ind, ub_ind]; } return rtn; @@ -139,22 +128,23 @@ matrix binomial_marginal(array[] int n, array[] int d, row_vector mu_glm, row_ve * @return Matrix containing continuous outcome in [0, 1], * bounds values and indicators */ -matrix poisson_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { - int J = num_elements(y); - matrix[J, 5] rtn; - row_vector[J] mu_glm_exp = exp(mu_glm); - - for (j in 1 : J) { - real u = u_raw[j]; - real lb = negative_infinity(); - real ub = inv_Phi(poisson_cdf(y[j] | mu_glm_exp[j])); - real lb_ind = 0; - real ub_ind = 1; - if (y[j] > 0) { - lb_ind = 1; - lb = inv_Phi(poisson_cdf(y[j] - 1 | mu_glm_exp[j])); +array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { + int N = rows(mu_glm); + int J = cols(mu_glm); + matrix[N, J] mu_glm_exp = exp(mu_glm); + array[2] matrix[N, J] rtn; + + for (n in 1:N) { + for (j in 1:J) { + real Ubound = poisson_cdf(y[n, j] | mu_glm_exp[n, j]); + real Lbound = 0; + if (y[n, j] > 0) { + Lbound = poisson_cdf(y[n, j] - 1 | mu_glm_exp[n, j]); + } + real UmL = Ubound - Lbound; + rtn[1][n, j] = Lbound + UmL * u_raw[n, j]; + rtn[2][n, j] = log(UmL); } - rtn[j] = [u, lb, ub, lb_ind, ub_ind]; } return rtn; @@ -169,19 +159,26 @@ matrix poisson_marginal(array[] int y, row_vector mu_glm, row_vector u_raw) { * * @param marginals Matrix of marginal calculations, constructed by concatenating the returns * from the _marginal functions - * @param tmvn_mu Vector of means for the Multivariate (Truncated) Gaussian, almost always 0 * @param L Cholesky Factor Correlation * @return Real log-probability */ -real gaussian_copula_cholesky_lpdf(matrix marginals, vector tmvn_mu, matrix L) { - int J = num_elements(tmvn_mu); - vector[J] u = marginals[ : , 1]; - vector[J] lb = marginals[ : , 2]; - vector[J] ub = marginals[ : , 3]; - vector[J] lb_ind = marginals[ : , 4]; - vector[J] ub_ind = marginals[ : , 5]; - - return multi_normal_cholesky_truncated_lpdf(u | tmvn_mu, L, lb, ub, lb_ind, ub_ind); +real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) { + int N = rows(marginals[1][1]); + int J = rows(L); + int N_marginals = size(marginals); + matrix[N, J] U; + real adj = 0; + int pos = 1; + for (m in 1:N_marginals) { + int curr_cols = cols(marginals[m][1]); + for (c in pos:(pos + curr_cols - 1)) { + U[ : , c] = marginals[m][1][ : , c - (pos - 1)]; + } + adj += sum(marginals[m][2]); + pos += curr_cols; + } + + return multi_normal_cholesky_copula_lpdf(U | L) + adj; } -/** @} */ \ No newline at end of file +/** @} */ diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions index 3cbc5faa..7aced250 100644 --- a/functions/copula/normal_copula.stanfunctions +++ b/functions/copula/normal_copula.stanfunctions @@ -21,7 +21,7 @@ * * @copyright Andre Pfeuffer, Sean Pinkney 2017, 2021 \n * https://groups.google.com/g/stan-users/c/hnUtkMYlLhQ/m/XdX3u1vDAAAJ \n - * Accessed and modified Feb. 5, 2021 + * Accessed and modified Feb. 5, 2021 * * @param u Real number on (0,1], not checked but function will return NaN * @param v Real number on (0,1], not checked but function will return NaN @@ -30,7 +30,7 @@ */ real normal_copula_lpdf(real u, real v, real rho) { real rho_sq = square(rho); - + return (0.5 * rho * (-2. * u * v + square(u) * rho + square(v) * rho)) / (-1. + rho_sq) - 0.5 * log1m(rho_sq); } @@ -52,12 +52,12 @@ real normal_copula_lpdf(real u, real v, real rho) { real normal_copula_vector_lpdf(vector u, vector v, real rho) { int N = num_elements(u); real rho_sq = square(rho); - + real a1 = 0.5 * rho; real a2 = rho_sq - 1; real a3 = 0.5 * log1m(rho_sq); real x = -2 * u' * v + rho * (dot_self(u) + dot_self(v)); - + return a1 * x / a2 - N * a3; } @@ -66,22 +66,22 @@ real normal_copula_vector_lpdf(vector u, vector v, real rho) { * * \f{aligned}{ * c(\mathbf{u}) &= \frac{\partial^d C}{\partial \mathbf{\Phi_1}\cdots \partial \mathbf{\Phi_d}} \\ - * &= \frac{1}{\sqrt{\det \Sigma}} \exp \Bigg(-\frac{1}{2} + * &= \frac{1}{\sqrt{\det \Sigma}} \exp \Bigg(-\frac{1}{2} * \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots \\ \mathbf{\Phi}^{-1}(u_d) \end{pmatrix}^T * (\Sigma^{-1} - I) * \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots \\ \mathbf{\Phi}^{-1}(u_d) \end{pmatrix} * \Bigg) \\ * &= \frac{1}{\sqrt{\det \Sigma}} \exp \bigg(-\frac{1}{2}(A^T\Sigma^{-1}A - A^TA) \bigg)\f} - * where + * where * \f[ * A = \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots \\ \mathbf{\Phi}^{-1}(u_d) \end{pmatrix}. * \f] * Note that \f$\det \Sigma = |LL^T| = |L |^2 = \big(\prod_{i=1}^d L_{i,i}\big)^2\f$ so \f$\sqrt{\det \Sigma} = \prod_{i=1}^d L_{i,i}\f$ - * and - * + * and + * *\f{aligned}{ * c(\mathbf{u}) &= \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}(A^T(LL^T)^{-1}A - A^TA)\bigg) \\ - * &= \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}\big((L^{-1}A)^TL^{-1}A - A^TA\big)\bigg) + * &= \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}\big((L^{-1}A)^TL^{-1}A - A^TA\big)\bigg) \f} * where \f$X = L^{-1}A\f$ can be solved with \code{.cpp} X = mdivide_left_tri_low(L, A)\endcode. * @@ -97,7 +97,7 @@ real multi_normal_copula_lpdf(matrix u, matrix L) { int N = cols(u); real inv_sqrt_det_log = sum(log(diagonal(L))); matrix[K, N] x = mdivide_left_tri_low(L, u); - + return -N * inv_sqrt_det_log - 0.5 * sum(columns_dot_self(x) - columns_dot_self(u)); } @@ -122,11 +122,19 @@ real bivariate_normal_copula_cdf(vector u, real rho) { real alpha_u = a * (pv / pu - rho); real alpha_v = a * (pu / pv - rho); real d = 0; - - if (u[1] < 0.5 && u[2] >= 0.5 || u[1] >= 0.5 && u[2] < 0.5) + + if (u[1] < 0.5 && u[2] >= 0.5 || u[1] >= 0.5 && u[2] < 0.5) d = 0.5; - + return avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d; } +real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) { + int N = rows(U); + int J = cols(U); + matrix[N,J] Q = to_matrix(inv_Phi(U), N, J ); + matrix[J,J] Gammainv = chol2inv(L); + return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U)); +} + /** @} */ \ No newline at end of file From c61c4cf85f7b719113071c4cb20e1df68a4fd881 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 10 Mar 2022 11:55:38 +0800 Subject: [PATCH 05/11] Test alternate likelihood --- .../centered_gaussian_copula_example.stan | 17 +++++-- .../centered_gaussian_copula.stanfunctions | 48 ++++++++++++------- functions/copula/normal_copula.stanfunctions | 10 ++-- 3 files changed, 50 insertions(+), 25 deletions(-) diff --git a/examples/copula/stan/centered_gaussian_copula_example.stan b/examples/copula/stan/centered_gaussian_copula_example.stan index 16fdccc3..52721db7 100644 --- a/examples/copula/stan/centered_gaussian_copula_example.stan +++ b/examples/copula/stan/centered_gaussian_copula_example.stan @@ -42,9 +42,20 @@ model { beta_n ~ normal(0.0, 10.0); beta_b ~ normal(0.0, 10.0); beta_p ~ normal(0.0, 10.0); - { normal_marginal(Yn, mu_n, sigma), - bernoulli_marginal(Yb, mu_b, uraw_b), - poisson_marginal(Yp, mu_p, uraw_p) } ~ centered_gaussian_copula_cholesky(L); + L ~ lkj_corr_cholesky(1.0); + + array[2] matrix[N,J[1]] normal_marg = normal_marginal(Yn, mu_n, sigma); + array[2] matrix[N,J[2]] bernoulli_marg = bernoulli_marginal(Yb, mu_b, uraw_b); + array[2] matrix[N,J[3]] poisson_marg = poisson_marginal(Yp, mu_p, uraw_p); + + // Append all Uniform RV's colwise + matrix[N,J_all] U = append_col(append_col(normal_marg[1], bernoulli_marg[1]), poisson_marg[1]); + + // Increment LL for copula + U ~ multi_normal_cholesky_copula(L); + + // Add jacobian adjustment + target += sum(normal_marg[2] + bernoulli_marg[2] + poisson_marg[2]); } generated quantities { corr_matrix[J_all] Gamma = multiply_lower_tri_self_transpose(L); diff --git a/functions/copula/centered_gaussian_copula.stanfunctions b/functions/copula/centered_gaussian_copula.stanfunctions index d9e90274..39c402d4 100644 --- a/functions/copula/centered_gaussian_copula.stanfunctions +++ b/functions/copula/centered_gaussian_copula.stanfunctions @@ -28,12 +28,16 @@ array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) { int N = rows(mu_glm); int J = cols(mu_glm); - matrix[N, J] adj = rep_matrix(0, N, J); + array[2] matrix[N, J] rtn; + for (j in 1:J) { - adj[1, j] = normal_lpdf(y[ : , j] | mu_glm[ : , j], sigma[j]); + for (n in 1:N) { + rtn[1][n,j] = Phi_approx((y[n,j] - mu_glm[n,j]) / sigma[j]); + rtn[2][n,j] = normal_lpdf(y[n,j] | mu_glm[n,j], sigma[j]); + } } - return { Phi_approx((y - mu_glm) ./ rep_matrix(transpose(sigma), N)), adj }; + return rtn; } /** @@ -57,15 +61,25 @@ array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) { array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { int N = rows(mu_glm); int J = cols(mu_glm); - matrix[N, J] mat_y = to_matrix(y); - matrix[N, J] mat_inv_y = fabs(mat_y - 1); - matrix[N, J] mu_glm_logit = inv_logit(mu_glm); - matrix[N, J] Lbound = mat_y .* (1 - mu_glm_logit); - matrix[N, J] Ubound = 1 - mat_inv_y .* mu_glm_logit; - matrix[N, J] UmL = Ubound - Lbound; + array[2] matrix[N, J] rtn; + + for (j in 1:J) { + for (n in 1:N) { + real Lbound = 0; + real Ubound = 1; + if (y[n,j] == 0) { + Ubound = 1 - mu_glm_logit[n,j]; + } else { + Lbound = 1 - mu_glm_logit[n,j]; + } + real UmL = Ubound - Lbound; + rtn[1][n,j] = Lbound + UmL * u_raw[n,j]; + rtn[2][n,j] = log(UmL); + } + } - return { Lbound + UmL .* u_raw, log(UmL) }; + return rtn; } /** @@ -94,8 +108,8 @@ array[] matrix binomial_marginal(array[,] int num, array[,] int den, matrix[N, J] mu_glm_logit = inv_logit(mu_glm); array[2] matrix[N, J] rtn; - for (n in 1:N) { - for (j in 1:J) { + for (j in 1:J) { + for (n in 1:N) { real Ubound = binomial_cdf(num[n, j] | den[n, j], mu_glm_logit[n, j]); real Lbound = 0; if (num[n, j] > 0) { @@ -134,8 +148,8 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { matrix[N, J] mu_glm_exp = exp(mu_glm); array[2] matrix[N, J] rtn; - for (n in 1:N) { - for (j in 1:J) { + for (j in 1:J) { + for (n in 1:N) { real Ubound = poisson_cdf(y[n, j] | mu_glm_exp[n, j]); real Lbound = 0; if (y[n, j] > 0) { @@ -168,11 +182,11 @@ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) int N_marginals = size(marginals); matrix[N, J] U; real adj = 0; - int pos = 1; + int pos = 0; for (m in 1:N_marginals) { int curr_cols = cols(marginals[m][1]); - for (c in pos:(pos + curr_cols - 1)) { - U[ : , c] = marginals[m][1][ : , c - (pos - 1)]; + for (c in 1:curr_cols) { + U[ : , (c + pos)] = marginals[m][1][ : , c]; } adj += sum(marginals[m][2]); pos += curr_cols; diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions index 7aced250..39dafb5d 100644 --- a/functions/copula/normal_copula.stanfunctions +++ b/functions/copula/normal_copula.stanfunctions @@ -130,11 +130,11 @@ real bivariate_normal_copula_cdf(vector u, real rho) { } real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) { - int N = rows(U); - int J = cols(U); - matrix[N,J] Q = to_matrix(inv_Phi(U), N, J ); - matrix[J,J] Gammainv = chol2inv(L); - return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U)); + int N = rows(U); + int J = cols(U); + matrix[N,J] Q = to_matrix( inv_Phi(to_vector(U)), N, J ); + matrix[J,J] Gammainv = chol2inv(L); + return -N * sum(log(diagonal(L))) - 0.5 * sum( add_diag(Gammainv, -1.0) .* crossprod(U) ); } /** @} */ \ No newline at end of file From 6ed9c72ca3d3d317018a056e2a764eaceaf50117 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 10 Mar 2022 12:06:40 +0800 Subject: [PATCH 06/11] Marginal construction to model block for debugging --- .../centered_gaussian_copula_example.stan | 34 ++++++++++++------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/examples/copula/stan/centered_gaussian_copula_example.stan b/examples/copula/stan/centered_gaussian_copula_example.stan index 52721db7..5b43a3d6 100644 --- a/examples/copula/stan/centered_gaussian_copula_example.stan +++ b/examples/copula/stan/centered_gaussian_copula_example.stan @@ -29,13 +29,28 @@ parameters { matrix[N, J[2]] uraw_b; // latent variables for bernoulli outcomes matrix[N, J[3]] uraw_p; // latent variables for Poisson outcomes } +transformed parameters { + array[2] matrix[N,J[1]] normal_marg; + array[2] matrix[N,J[2]] bernoulli_marg; + array[2] matrix[N,J[3]] poisson_marg; + + { + // initialize variables + vector[J[1]] sigma = sqrt(sigmasq); // stdev + // Calculate the means for each regression + matrix[N, J[1]] mu_n = Xn * rep_matrix(beta_n, J[1]); + matrix[N, J[2]] mu_b = Xb * rep_matrix(beta_b, J[2]); + matrix[N, J[3]] mu_p = Xp * rep_matrix(beta_p, J[3]); + + normal_marg = normal_marginal(Yn, mu_n, sigma); + bernoulli_marg = bernoulli_marginal(Yb, mu_b, uraw_b); + poisson_marg = poisson_marginal(Yp, mu_p, uraw_p); + } + + // Append all Uniform RV's colwise + matrix[N,J_all] U = append_col(append_col(normal_marg[1], bernoulli_marg[1]), poisson_marg[1]); +} model { - // initialize variables - vector[J[1]] sigma = sqrt(sigmasq); // stdev - // Calculate the means for each regression - matrix[N, J[1]] mu_n = Xn * rep_matrix(beta_n, J[1]); - matrix[N, J[2]] mu_b = Xb * rep_matrix(beta_b, J[2]); - matrix[N, J[3]] mu_p = Xp * rep_matrix(beta_p, J[3]); sigmasq ~ inv_gamma(1.0e-4, 1.0e-4); // priors for regression coefficients @@ -44,13 +59,6 @@ model { beta_p ~ normal(0.0, 10.0); L ~ lkj_corr_cholesky(1.0); - array[2] matrix[N,J[1]] normal_marg = normal_marginal(Yn, mu_n, sigma); - array[2] matrix[N,J[2]] bernoulli_marg = bernoulli_marginal(Yb, mu_b, uraw_b); - array[2] matrix[N,J[3]] poisson_marg = poisson_marginal(Yp, mu_p, uraw_p); - - // Append all Uniform RV's colwise - matrix[N,J_all] U = append_col(append_col(normal_marg[1], bernoulli_marg[1]), poisson_marg[1]); - // Increment LL for copula U ~ multi_normal_cholesky_copula(L); From 9b4cd5b96968c584b885c8cde6ce1bf7e95f57d7 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 11 Mar 2022 17:35:14 +0800 Subject: [PATCH 07/11] Update inv_Phi usage --- .../centered_gaussian_copula_example.stan | 31 ++++++------------- .../centered_gaussian_copula.stanfunctions | 18 ++++++----- functions/copula/normal_copula.stanfunctions | 18 +++++++---- 3 files changed, 32 insertions(+), 35 deletions(-) diff --git a/examples/copula/stan/centered_gaussian_copula_example.stan b/examples/copula/stan/centered_gaussian_copula_example.stan index 5b43a3d6..f8492a1c 100644 --- a/examples/copula/stan/centered_gaussian_copula_example.stan +++ b/examples/copula/stan/centered_gaussian_copula_example.stan @@ -30,27 +30,14 @@ parameters { matrix[N, J[3]] uraw_p; // latent variables for Poisson outcomes } transformed parameters { - array[2] matrix[N,J[1]] normal_marg; - array[2] matrix[N,J[2]] bernoulli_marg; - array[2] matrix[N,J[3]] poisson_marg; - - { - // initialize variables - vector[J[1]] sigma = sqrt(sigmasq); // stdev - // Calculate the means for each regression - matrix[N, J[1]] mu_n = Xn * rep_matrix(beta_n, J[1]); - matrix[N, J[2]] mu_b = Xb * rep_matrix(beta_b, J[2]); - matrix[N, J[3]] mu_p = Xp * rep_matrix(beta_p, J[3]); - - normal_marg = normal_marginal(Yn, mu_n, sigma); - bernoulli_marg = bernoulli_marginal(Yb, mu_b, uraw_b); - poisson_marg = poisson_marginal(Yp, mu_p, uraw_p); - } - - // Append all Uniform RV's colwise - matrix[N,J_all] U = append_col(append_col(normal_marg[1], bernoulli_marg[1]), poisson_marg[1]); } model { + // initialize variables + vector[J[1]] sigma = sqrt(sigmasq); // stdev + // Calculate the means for each regression + matrix[N, J[1]] mu_n = Xn * rep_matrix(beta_n, J[1]); + matrix[N, J[2]] mu_b = Xb * rep_matrix(beta_b, J[2]); + matrix[N, J[3]] mu_p = Xp * rep_matrix(beta_p, J[3]); sigmasq ~ inv_gamma(1.0e-4, 1.0e-4); // priors for regression coefficients @@ -60,10 +47,10 @@ model { L ~ lkj_corr_cholesky(1.0); // Increment LL for copula - U ~ multi_normal_cholesky_copula(L); + { normal_marginal(Yn, mu_n, sigma), + bernoulli_marginal(Yb, mu_b, uraw_b), + poisson_marginal(Yp, mu_p, uraw_p) } ~ centered_gaussian_copula_cholesky(L); - // Add jacobian adjustment - target += sum(normal_marg[2] + bernoulli_marg[2] + poisson_marg[2]); } generated quantities { corr_matrix[J_all] Gamma = multiply_lower_tri_self_transpose(L); diff --git a/functions/copula/centered_gaussian_copula.stanfunctions b/functions/copula/centered_gaussian_copula.stanfunctions index 39c402d4..0401506c 100644 --- a/functions/copula/centered_gaussian_copula.stanfunctions +++ b/functions/copula/centered_gaussian_copula.stanfunctions @@ -74,7 +74,7 @@ array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { Lbound = 1 - mu_glm_logit[n,j]; } real UmL = Ubound - Lbound; - rtn[1][n,j] = Lbound + UmL * u_raw[n,j]; + rtn[1][n,j] = inv_Phi(Lbound + UmL * u_raw[n,j]); rtn[2][n,j] = log(UmL); } } @@ -116,7 +116,7 @@ array[] matrix binomial_marginal(array[,] int num, array[,] int den, Lbound = binomial_cdf(num[n, j] - 1 | den[n, j], mu_glm_logit[n, j]); } real UmL = Ubound - Lbound; - rtn[1][n, j] = Lbound + UmL * u_raw[n, j]; + rtn[1][n, j] = inv_Phi(Lbound + UmL * u_raw[n, j]); rtn[2][n, j] = log(UmL); } } @@ -156,7 +156,7 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { Lbound = poisson_cdf(y[n, j] - 1 | mu_glm_exp[n, j]); } real UmL = Ubound - Lbound; - rtn[1][n, j] = Lbound + UmL * u_raw[n, j]; + rtn[1][n, j] = inv_Phi(Lbound + UmL * u_raw[n, j]); rtn[2][n, j] = log(UmL); } } @@ -171,19 +171,22 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { * @copyright Ethan Alt, Sean Pinkney 2021 \n * * - * @param marginals Matrix of marginal calculations, constructed by concatenating the returns - * from the _marginal functions + * @param marginals Nested arrays of matrices from marginal calculations * @param L Cholesky Factor Correlation * @return Real log-probability */ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) { + // Extract dimensions of final outcome matrix int N = rows(marginals[1][1]); int J = rows(L); - int N_marginals = size(marginals); matrix[N, J] U; + + // Iterate through marginal arrays, concatenating the outcome matrices by column + // and aggregating the log-likelihoods (from continuous marginals) and jacobian + // adjustments (from discrete marginals) real adj = 0; int pos = 0; - for (m in 1:N_marginals) { + for (m in 1:size(marginals)) { int curr_cols = cols(marginals[m][1]); for (c in 1:curr_cols) { U[ : , (c + pos)] = marginals[m][1][ : , c]; @@ -192,6 +195,7 @@ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) pos += curr_cols; } + // Return the sum of the log-probability for copula outcomes and jacobian adjustments return multi_normal_cholesky_copula_lpdf(U | L) + adj; } diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions index 39dafb5d..f4f9c1fd 100644 --- a/functions/copula/normal_copula.stanfunctions +++ b/functions/copula/normal_copula.stanfunctions @@ -129,12 +129,18 @@ real bivariate_normal_copula_cdf(vector u, real rho) { return avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d; } +/** + * Multivariate Normal Copula lpdf (Cholesky parameterisation) + * + * @param U Matrix of outcomes from marginal calculations + * @param L Cholesky factor of the correlation/covariance matrix + * @return log density of distribution + */ real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) { - int N = rows(U); - int J = cols(U); - matrix[N,J] Q = to_matrix( inv_Phi(to_vector(U)), N, J ); - matrix[J,J] Gammainv = chol2inv(L); - return -N * sum(log(diagonal(L))) - 0.5 * sum( add_diag(Gammainv, -1.0) .* crossprod(U) ); + int N = rows(U); + int J = cols(U); + matrix[J,J] Gammainv = chol2inv(L); + return -N * sum(log(diagonal(L))) - inv(sum(add_diag(Gammainv, -1.0) .* crossprod(U))); } -/** @} */ \ No newline at end of file +/** @} */ From 5ca49db7bd7ed1475d52fef14b0757260796330e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 11 Mar 2022 22:25:06 +0800 Subject: [PATCH 08/11] Update vectorisation, add doc --- .../centered_gaussian_copula.stanfunctions | 82 ++++++++----------- functions/copula/normal_copula.stanfunctions | 14 ---- 2 files changed, 35 insertions(+), 61 deletions(-) diff --git a/functions/copula/centered_gaussian_copula.stanfunctions b/functions/copula/centered_gaussian_copula.stanfunctions index 0401506c..13250b3e 100644 --- a/functions/copula/centered_gaussian_copula.stanfunctions +++ b/functions/copula/centered_gaussian_copula.stanfunctions @@ -19,22 +19,22 @@ * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param y row_vector of normal outcomes + * @param y matrix of normal outcomes * @param mu_glm row_vector of regression means - * @param sigma vector of outcome SD's - * @return Matrix containing continuous outcome in [0, 1], - * bounds values and indicators + * @param matrix vector of outcome SD's + * @return 2D array of matrices containing the random variables + * and jacobian adjustments */ array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) { int N = rows(mu_glm); int J = cols(mu_glm); array[2] matrix[N, J] rtn; + // Initialise the jacobian adjustments to zero, as vectorised lpdf will be used + rtn[2] = rep_matrix(0, N, J); for (j in 1:J) { - for (n in 1:N) { - rtn[1][n,j] = Phi_approx((y[n,j] - mu_glm[n,j]) / sigma[j]); - rtn[2][n,j] = normal_lpdf(y[n,j] | mu_glm[n,j], sigma[j]); - } + rtn[1][ : , j] = (y[: , j] - mu_glm[ : , j]) / sigma[j]; + rtn[2][1, j] = normal_lpdf(y[ : , j] | mu_glm[ : , j], sigma[j]); } return rtn; @@ -52,34 +52,22 @@ array[] matrix normal_marginal(matrix y, matrix mu_glm, vector sigma) { * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param y int[] Array of binary outcomes - * @param mu_glm row_vector of regression means - * @param u_raw row_vector of nuisance latent variables - * @return Matrix containing continuous outcome in [0, 1], - * bounds values and indicators + * @param y int[,] 2d array of binary outcomes + * @param mu_glm matrix of regression means + * @param u_raw matrix of nuisance latent variables + * @return 2D array of matrices containing the random variables + * and jacobian adjustments */ array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { int N = rows(mu_glm); int J = cols(mu_glm); - matrix[N, J] mu_glm_logit = inv_logit(mu_glm); - array[2] matrix[N, J] rtn; + matrix[N, J] matrix_y = to_matrix(y); + matrix[N, J] mu_glm_logit = 1 - inv_logit(mu_glm); - for (j in 1:J) { - for (n in 1:N) { - real Lbound = 0; - real Ubound = 1; - if (y[n,j] == 0) { - Ubound = 1 - mu_glm_logit[n,j]; - } else { - Lbound = 1 - mu_glm_logit[n,j]; - } - real UmL = Ubound - Lbound; - rtn[1][n,j] = inv_Phi(Lbound + UmL * u_raw[n,j]); - rtn[2][n,j] = log(UmL); - } - } + matrix[N, J] Lbound = matrix_y .* mu_glm_logit; + matrix[N, J] UmL = fabs(matrix_y - mu_glm_logit); - return rtn; + return { inv_Phi(Lbound + UmL .* u_raw), log(UmL) }; } /** @@ -94,12 +82,12 @@ array[] matrix bernoulli_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param n int[] Array of numerator integers - * @param d int[] Array of denominator integers - * @param mu_glm row_vector of regression means - * @param u_raw row_vector of nuisance latent variables - * @return Matrix containing continuous outcome in [0, 1], - * bounds values and indicators + * @param num int[,] 2D array of numerator integers + * @param den int[,] 2D array of denominator integers + * @param mu_glm matrix of regression means + * @param u_raw matrix of nuisance latent variables + * @return 2D array of matrices containing the random variables + * and jacobian adjustments */ array[] matrix binomial_marginal(array[,] int num, array[,] int den, matrix mu_glm, matrix u_raw) { @@ -136,11 +124,11 @@ array[] matrix binomial_marginal(array[,] int num, array[,] int den, * * @copyright Ethan Alt, Sean Pinkney 2021 \n * - * @param y int[] Array of integer counts - * @param mu_glm row_vector of regression means - * @param u_raw row_vector of nuisance latent variables - * @return Matrix containing continuous outcome in [0, 1], - * bounds values and indicators + * @param y int[,] 2D array of integer counts + * @param mu_glm matrix of regression means + * @param u_raw matrix of nuisance latent variables + * @return 2D array of matrices containing the random variables + * and jacobian adjustments */ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { int N = rows(mu_glm); @@ -179,24 +167,24 @@ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) // Extract dimensions of final outcome matrix int N = rows(marginals[1][1]); int J = rows(L); - matrix[N, J] U; + matrix[J, N] U; // Iterate through marginal arrays, concatenating the outcome matrices by column // and aggregating the log-likelihoods (from continuous marginals) and jacobian // adjustments (from discrete marginals) real adj = 0; - int pos = 0; + int pos = 1; for (m in 1:size(marginals)) { int curr_cols = cols(marginals[m][1]); - for (c in 1:curr_cols) { - U[ : , (c + pos)] = marginals[m][1][ : , c]; - } + + U[pos:(pos + curr_cols - 1), : ] = transpose(marginals[m][1]); + adj += sum(marginals[m][2]); pos += curr_cols; } // Return the sum of the log-probability for copula outcomes and jacobian adjustments - return multi_normal_cholesky_copula_lpdf(U | L) + adj; + return multi_normal_copula_lpdf(U | L) + adj; } /** @} */ diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions index f4f9c1fd..016b28d9 100644 --- a/functions/copula/normal_copula.stanfunctions +++ b/functions/copula/normal_copula.stanfunctions @@ -129,18 +129,4 @@ real bivariate_normal_copula_cdf(vector u, real rho) { return avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d; } -/** - * Multivariate Normal Copula lpdf (Cholesky parameterisation) - * - * @param U Matrix of outcomes from marginal calculations - * @param L Cholesky factor of the correlation/covariance matrix - * @return log density of distribution - */ -real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) { - int N = rows(U); - int J = cols(U); - matrix[J,J] Gammainv = chol2inv(L); - return -N * sum(log(diagonal(L))) - inv(sum(add_diag(Gammainv, -1.0) .* crossprod(U))); -} - /** @} */ From b84e148febc0f08f6d782ae3fe237cec045cb9ee Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 11 Mar 2022 22:31:08 +0800 Subject: [PATCH 09/11] Update copula cholesky call --- .../copula/centered_gaussian_copula.stanfunctions | 6 +++--- functions/copula/normal_copula.stanfunctions | 14 ++++++++++++++ 2 files changed, 17 insertions(+), 3 deletions(-) diff --git a/functions/copula/centered_gaussian_copula.stanfunctions b/functions/copula/centered_gaussian_copula.stanfunctions index 13250b3e..ad7669a5 100644 --- a/functions/copula/centered_gaussian_copula.stanfunctions +++ b/functions/copula/centered_gaussian_copula.stanfunctions @@ -167,7 +167,7 @@ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) // Extract dimensions of final outcome matrix int N = rows(marginals[1][1]); int J = rows(L); - matrix[J, N] U; + matrix[N, J] U; // Iterate through marginal arrays, concatenating the outcome matrices by column // and aggregating the log-likelihoods (from continuous marginals) and jacobian @@ -177,14 +177,14 @@ real centered_gaussian_copula_cholesky_lpdf(array[,] matrix marginals, matrix L) for (m in 1:size(marginals)) { int curr_cols = cols(marginals[m][1]); - U[pos:(pos + curr_cols - 1), : ] = transpose(marginals[m][1]); + U[ : , pos:(pos + curr_cols - 1)] = marginals[m][1]; adj += sum(marginals[m][2]); pos += curr_cols; } // Return the sum of the log-probability for copula outcomes and jacobian adjustments - return multi_normal_copula_lpdf(U | L) + adj; + return multi_normal_cholesky_copula_lpdf(U | L) + adj; } /** @} */ diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions index 016b28d9..01f8ae4f 100644 --- a/functions/copula/normal_copula.stanfunctions +++ b/functions/copula/normal_copula.stanfunctions @@ -129,4 +129,18 @@ real bivariate_normal_copula_cdf(vector u, real rho) { return avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d; } +/** + * Multivariate Normal Copula lpdf (Cholesky parameterisation) + * + * @param U Matrix of outcomes from marginal calculations + * @param L Cholesky factor of the correlation/covariance matrix + * @return log density of distribution + */ +real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) { + int N = rows(U); + int J = cols(U); + matrix[J,J] Gammainv = chol2inv(L); + return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U)); +} + /** @} */ From 12a47008aea6a6c2c86f96f4eab0c3456994545e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Sat, 12 Mar 2022 11:26:29 +0800 Subject: [PATCH 10/11] Remove stray include path --- examples/copula/R/centered_gaussian_copula_example.R | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/copula/R/centered_gaussian_copula_example.R b/examples/copula/R/centered_gaussian_copula_example.R index 633e8745..a1aef9b1 100644 --- a/examples/copula/R/centered_gaussian_copula_example.R +++ b/examples/copula/R/centered_gaussian_copula_example.R @@ -3,8 +3,7 @@ library(cmdstanr) fp <- file.path("./examples/copula/stan/centered_gaussian_copula_example.stan") mod <- cmdstan_model(fp, include_paths = c("./functions", "./functions/copula", - "./functions/distribution", - "../functions")) + "./functions/distribution")) ## obtain non-extreme random correlation matrix Gamma = matrix(1, 3, 3) From 119d9a4d5ddc7ca2c301bcc50291073aa2b17951 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Sat, 12 Mar 2022 07:14:55 -0500 Subject: [PATCH 11/11] some doc updates --- .../R/centered_gaussian_copula_example.R | 19 ---------- .../centered_gaussian_copula.stanfunctions | 13 +++---- functions/copula/normal_copula.stanfunctions | 38 ++++++++++++++++++- 3 files changed, 43 insertions(+), 27 deletions(-) diff --git a/examples/copula/R/centered_gaussian_copula_example.R b/examples/copula/R/centered_gaussian_copula_example.R index a1aef9b1..ac3abef7 100644 --- a/examples/copula/R/centered_gaussian_copula_example.R +++ b/examples/copula/R/centered_gaussian_copula_example.R @@ -111,27 +111,9 @@ get_data = function(formula.list, family.list, data) { 'Xindx' = Xindx, 'Kj' = Kj ) - - ## fit stan model - # fit = rstan::sampling(stanmod, data = standat, ...) - # fit = stanmod$sample(data = standat, - # parallel_chains = 4, - # iter_warmup = 1000, - # iter_sampling = 1000) - # ## rename betas for clarity about which equation each comes from, e.g., beta[j][0:(Kj-1)] - # beta.old.indx = which(grepl('beta', names(fit))) - # beta.new = character() - # for( j in 1:J ) { - # beta.new = c(beta.new, paste0(distnames[j], '[', seq(0,Kj[j] - 1), ']') ) - # } - # names(fit)[beta.old.indx] = beta.new - # - ## return stan object return(standat) } - - stan_data <- get_data(formula.list[1:3], family.list[1:3], data) stan_data[["J"]] <- rep(1, 3) @@ -146,4 +128,3 @@ mod_out <- mod$sample(data = stan_data, parallel_chains = 2, iter_warmup = 400, iter_sampling = 400) - diff --git a/functions/copula/centered_gaussian_copula.stanfunctions b/functions/copula/centered_gaussian_copula.stanfunctions index ad7669a5..8151e6f8 100644 --- a/functions/copula/centered_gaussian_copula.stanfunctions +++ b/functions/copula/centered_gaussian_copula.stanfunctions @@ -18,7 +18,7 @@ * copula. * * @copyright Ethan Alt, Sean Pinkney 2021 \n - * + * * @param y matrix of normal outcomes * @param mu_glm row_vector of regression means * @param matrix vector of outcome SD's @@ -118,10 +118,10 @@ array[] matrix binomial_marginal(array[,] int num, array[,] int den, * Poisson marginal for mixed continuous-discrete Gaussian * copula. * - * The lb and ub: - * Always upper bound at inv_Phi( y, mu ) - * If y != 0, lower bound at inv_Phi(y-1, mu) - * + * The lower-bound and upper-bound: \n + * The upper-bound is always at \code{.cpp} inv_Phi( y, mu ) \endcode \n + * If \f$y \ne 0\f$, lower-bound at \code{.cpp} inv_Phi(y-1, mu) \endcode. \n + * At \f$y = 0\f$ the lower-bound is \f$ 0 \f$. * @copyright Ethan Alt, Sean Pinkney 2021 \n * * @param y int[,] 2D array of integer counts @@ -156,8 +156,7 @@ array[] matrix poisson_marginal(array[,] int y, matrix mu_glm, matrix u_raw) { /** * Mixed Copula Log-Probability Function * - * @copyright Ethan Alt, Sean Pinkney 2021 \n - * + * @copyright Andrew Johnson 2022 \n * * @param marginals Nested arrays of matrices from marginal calculations * @param L Cholesky Factor Correlation diff --git a/functions/copula/normal_copula.stanfunctions b/functions/copula/normal_copula.stanfunctions index 01f8ae4f..4a2a4811 100644 --- a/functions/copula/normal_copula.stanfunctions +++ b/functions/copula/normal_copula.stanfunctions @@ -129,6 +129,43 @@ real bivariate_normal_copula_cdf(vector u, real rho) { return avg_uv - owens_t(pu, alpha_u) - owens_t(pv, alpha_v) - d; } +/** + * Multivariate Normal Copula log density (Cholesky parameterisation) + * + * \f{aligned}{ + * c(\mathbf{u}) &= \frac{\partial^d C}{\partial \mathbf{\Phi_1}\cdots \partial \mathbf{\Phi_d}} \\ + * &= \prod_{i=1}^{n} \lvert \Sigma \rvert^{-1/2} \exp \left\{-\frac{1}{2} + * \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots \\ \mathbf{\Phi}^{-1}(u_d) \end{pmatrix}^T + * (\Sigma^{-1} - I) + * \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots \\ \mathbf{\Phi}^{-1}(u_d) \end{pmatrix} + * \right\} \\ + * &= \lvert \Sigma \rvert^{-n/2} \exp\left\{ -\frac{1}{2} \text{tr}\left( \left[ \Sigma^{-1} - I \right] \sum_{i=1}^n \Phi^{-1}(u_i) \Phi^{-1}(u_i)' \right) \right\} \hspace{0.5cm} \\ + * &= \lvert L \rvert^{-n} \exp\left\{ -\frac{1}{2} \text{tr}\left( \left[ L^{-T}L^{-1}- I \right] Q'Q \right) \\ + * &= \lvert L \rvert^{-n} \exp\left\{ -\frac{1}{2} \text{tr}\left( L^{-T}L^{-1} Q'Q - Q'Q \right) \right\} \\ + * &= \lvert L \rvert^{-n} \exp\left\{ -\frac{1}{2} \text{tr}\left( QL^{-T}L^{-1}Q' \right) + \text{tr} \left( Q'Q \right) \right\} \\ + * &= \lvert L \rvert^{-n} \exp\left\{ -\frac{1}{2} \left[ \text{tr}\left( [L^{-1}Q'][L^{-1}Q']' \right) + \text{tr} \left( Q'Q \right) \right]\right\} + * \f} + * where + * \f[ + * Q_{n \times d} = \begin{pmatrix} \mathbf{\Phi}^{-1}(u_1) \\ \vdots \\ \mathbf{\Phi}^{-1}(u_d) \end{pmatrix}. + * \f] + * Note that \f$\det \Sigma = |LL^T| = |L |^2 = \big(\prod_{i=1}^d L_{i,i}\big)^2\f$ so \f$\sqrt{\det \Sigma} = \prod_{i=1}^d L_{i,i}\f$ + * and + * + *\f{aligned}{ + * c(\mathbf{u}) &= \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}(A^T(LL^T)^{-1}A - A^TA)\bigg) \\ + * &= \Bigg(\prod_{i=1}^d L_{i,i}\Bigg)^{-1} \exp \bigg(-\frac{1}{2}\big((L^{-1}A)^TL^{-1}A - A^TA\big)\bigg) + \f} + * where \f$X = L^{-1}A\f$ can be solved with \code{.cpp} X = mdivide_left_tri_low(L, A)\endcode. + * + * @copyright Ethan Alt, Sean Pinkney 2021 + * + * @param U Matrix where the number of rows equal observations and the columns are equal to the number of outcomes + * @param L Cholesky factor of the correlation matrix + * @return log density of the distribution + */ + + /** * Multivariate Normal Copula lpdf (Cholesky parameterisation) * @@ -142,5 +179,4 @@ real multi_normal_cholesky_copula_lpdf(matrix U, matrix L) { matrix[J,J] Gammainv = chol2inv(L); return -N * sum(log(diagonal(L))) - 0.5 * sum(add_diag(Gammainv, -1.0) .* crossprod(U)); } - /** @} */