data { int K; int D; int N; int NN; int M; //samples for marginal effects est int y[N, D]; vector[K] x[N]; matrix[N,K] x_data; matrix[NN,K] x_predict; matrix[M,K] x_m_1; //marginal effect AREA matrix[M,K] x_m_2; //marginal effect CROP matrix[M,K] x_m_3; //marginal effect SEP } parameters { vector[D] beta0; matrix[D,K] beta; cholesky_factor_corr[D] L_Omega; real u[N,D]; // nuisance that absorbs inequality constraints } model { L_Omega ~ lkj_corr_cholesky(2); beta0 ~ student_t(1, 0, 3); to_vector(beta) ~ student_t(1, 0, 1.25); { // likelihood for (n in 1:N){ vector[D] mu; vector[D] z; real prev; mu = beta0 + (beta * x[n]); prev = 0; for (d in 1:D) { // Phi and inv_Phi may overflow and / or be numerically inaccurate real bound; // threshold at which utility = 0 bound = Phi( -(mu[d] + prev) / L_Omega[d, d] ); if (y[n, d] == 1) { real t; t = bound + (1 - bound) * u[n, d]; z[d] = inv_Phi(t); // implies utility is positive target += log1m(bound); // Jacobian adjustment } else { real t; t = bound * u[n, d]; z[d] = inv_Phi(t); // implies utility is negative target += log(bound); // Jacobian adjustment } if (d < D) prev = L_Omega[d+1, 1:d] * head(z, d); // Jacobian adjustments imply z is truncated standard normal // thus utility --- mu + L_Omega * z --- is truncated multivariate normal } } } } generated quantities { matrix[D,N] z_new; matrix[D,NN] z_predict; matrix[D,M] z_m_1; matrix[D,M] z_m_2; matrix[D,M] z_m_3; corr_matrix[D] Omega; Omega = multiply_lower_tri_self_transpose(L_Omega); { matrix[K,D] t_beta; matrix[D,N] beta_x; matrix[D,NN] beta_x_predict; matrix[D,M] beta_x_m1; matrix[D,M] beta_x_m2; matrix[D,M] beta_x_m3; t_beta = beta'; for (d in 1:D){ beta_x[d,] = beta0[d] + (x_data * t_beta[,d])'; beta_x_predict[d,] = beta0[d] + (x_predict * t_beta[,d])'; beta_x_m1[d,] = beta0[d] + (x_m_1 * t_beta[,d])'; beta_x_m2[d,] = beta0[d] + (x_m_2 * t_beta[,d])'; beta_x_m3[d,] = beta0[d] + (x_m_3 * t_beta[,d])'; } for (n in 1:N) z_new[,n] = multi_normal_cholesky_rng(beta_x[,n], L_Omega); for (nn in 1:NN){ z_predict[,nn] = multi_normal_cholesky_rng(beta_x_predict[,nn], L_Omega); } for (m in 1:M) { z_m_1[,m] = multi_normal_cholesky_rng(beta_x_m1[,m], L_Omega); z_m_2[,m] = multi_normal_cholesky_rng(beta_x_m2[,m], L_Omega); z_m_3[,m] = multi_normal_cholesky_rng(beta_x_m3[,m], L_Omega); } } }