Non-centered parameterization of the multivariate normal distribution

Modeling
, ,
1

Hi,

I am pretty new to stan and am trying to fit truncated multivariate normal where the random variables generated must be positive. But I am facing some issues with it.

  1. When I am trying to use the case below, I cant use the constraint to produce positive random variables and they need to have different upper bounds, since the variance co-variance matrix contains both positive and negative values . For example, use

say, k = 3


data {

  int<lower=1> N;                                            
  int<lower=1> K;                                  
  
  matrix[N, K] x;                                 
  vector[K] zero_vec;
  matrix[K, K] I_mat; 
  
  vector[N] Y_fit;                                          
  
}

parameters {
  
  //beta
  vector alpha;                                           
  vector<lower=0, upper=1> [K] beta;            
  
   
  real<lower=0> tau; 
  real mu_alpha;                                              
  vector<lower=0> [K] mu_beta;
  

  real<lower=0> sigma_alpha;                                   
  vector<lower=0> [K] sigma_beta;                 
  vector<lower=0> [K] sigma_beta_mu;
  
  cholesky_factor_corr [K] L_Chol;  
  cholesky_factor_corr [K] L_Chol_mu;  
  
}

transformed parameters {
  matrix [K, K] Var_Cov;
  matrix [K, K] Var_Cov_mu;

  vector [K] beta_std[J];
  

  #Var_Cov = diag_pre_multiply(sigma_beta, L_Chol) * diag_pre_multiply(sigma_beta, L_Chol)';

  Var_Cov_mu = diag_pre_multiply(sigma_beta_mu, L_Chol_mu) * diag_pre_multiply(sigma_beta_mu, L_Chol_mu)';
  
  beta_std = mu_beta + diag_pre_multiply(sigma_beta, L_Chol) * beta;  
                    
}

model {
  
  mu_alpha ~ student_t(1, 2, 2);
  sigma_alpha ~ cauchy(0, 2.5);
  
  // Prior
  
  beta ~ multi_normal(zero_vec, I_mat);
  mu_beta ~ multi_normal(zero_vec, Var_Cov_mu/0.01);

  L_Chol ~ lkj_corr_cholesky(2);
  L_Chol_mu ~ lkj_corr_cholesky(2);

  sigma_beta ~ cauchy(0, 2.5);
  sigma_beta_mu ~ cauchy(0, 2.5);
  
  alpha ~ normal(mu_alpha, sigma_alpha);
  sigma ~ normal(0,100);
  
  for(i in 1:N){
    Y_fit[i] ~ normal(alpha[i] + beta_std[i, 1] * x[i, 1] + beta_std[i, 2] * x[i, 2]
                      + beta_std[i, 3] * x[i, 3], sigma); 
    
  }
  
}
  1. And if I use this code then the random variables are all positive, but I am having divergences in iteration. I changed my delta to 0.999 and above with 6 to 7 nines and tree depth 30 - 35, but was of no help in case for divergence.

say, k = 3

data {

  int<lower=1> N;                                            
  int<lower=1> K;                                  
  
  matrix[N, K] x;                                 
  vector[K] zero_vec;
  matrix[K, K] I_mat; 
  
  vector[N] Y_fit;                                          
  
}

parameters {
  
  //beta
  vector alpha;                                           
  vector<lower=0, upper=1> [K] beta;            
  
   
  real<lower=0> tau; 
  real mu_alpha;                                              
  vector<lower=0> [K] mu_beta;
  

  real<lower=0> sigma_alpha;                                   
  vector<lower=0> [K] sigma_beta;                 
  vector<lower=0> [K] sigma_beta_mu;
  
  cholesky_factor_corr [K] L_Chol;  
  cholesky_factor_corr [K] L_Chol_mu;  
  
}

transformed parameters {
  matrix [K, K] Var_Cov;
  matrix [K, K] Var_Cov_mu;

  vector [K] beta_std[J];
  
  beta_std = beta;

  Var_Cov = diag_pre_multiply(sigma_beta, L_Chol) * diag_pre_multiply(sigma_beta, L_Chol)';
  Var_Cov_mu = diag_pre_multiply(sigma_beta_mu, L_Chol_mu) * diag_pre_multiply(sigma_beta_mu, L_Chol_mu)';
                    
}

model {
  
  mu_alpha ~ student_t(1, 2, 2);
  sigma_alpha ~ cauchy(0, 2.5);
  
  // Prior
  
  beta ~ multi_normal(mu_beta, Var_Cov);
  mu_beta ~ multi_normal(zero_vec, Var_Cov_mu/0.01);

  L_Chol ~ lkj_corr_cholesky(2);
  L_Chol_mu ~ lkj_corr_cholesky(2);

  sigma_beta ~ cauchy(0, 2.5);
  sigma_beta_mu ~ cauchy(0, 2.5);
  
  alpha ~ normal(mu_alpha, sigma_alpha);
  sigma ~ normal(0,100);
  
  for(i in 1:N){
    Y_fit[i] ~ normal(alpha[i] + beta_std[i, 1] * x[i, 1] + beta_std[i, 2] * x[i, 2]
                      + beta_std[i, 3] * x[i, 3], sigma); 
    
  }
  
}

can you please help me with this.

thank you,
BG

2

Hi,
this is not really my area of expertise, but since nobody else replied, some quick notes:

  • It would be easier to answer your question, if you stripped your Stan code to the minimal working example showing your problem, which here would probably be just a single truncated multivariete normal with everything except the variables of interest given as data. I am for example confused as to which of the multi_normal calls should be truncated, which part of the code is related to the (attempted) truncation and which is not and how do the two versions of the code differ.

  • You can make your code simpler and more efficient by using multi_normal_cholesky, which takes the Cholesky decomposition of the variance-covariance matrix directly.

Best of luck with the model!