Diagonal covariance matrix

Modeling
1

I want to write in my stan code a diagonal covariance matrix D; when I did the following, stan gave me an error of
Rejecting initial value:
Error evaluating the log probability at the initial value.
multi_normal_lpdf: LDLT_Factor of covariance parameter is not positive definite. last conditional variance is 0.

parameters {
    vector<lower=0>[K] sigma_d[N];
} 
transformed parameters {
  matrix[K, K] D[N]; 
  for (n in 1:N){
    D[n] = diag_matrix(sigma_d[n]);
  }
}
model {
   for(n in 1:N){
   Y[n] ~ multi_normal(mu[n], D[n]);
    }
}

when I used Cholesky decomposition, the error above went away, but D is not diagonal with zeros on off-diagnonal, and I also received this error: [1] “Error in sampler$call_sampler(args_list[[i]]) : " " c++ exception (unknown reason)”

My code with Cholesky decomposition is below:

parameters {
  cholesky_factor_corr[K] D_Omega[N]; 
  vector<lower=0>[K] D_tau; 
} 
transformed parameters {
  matrix[K, K] D[N]; 
  for (n in 1:N){
    D[n] = diag_pre_multiply(D_tau, D_Omega[n]);
  }
}
model {
   for(n in 1:N){
   Y[n] ~ multi_normal_cholesky(mu[n], D[n]);
    }
}

So how can I use Cholesky decomposition to define a diagonal covariance matrix in Stan?

2

Don’t use a diagonal covariance matrix with multi_normal_*. If the covariance matrix is diagonal, then the random variables are conditionally independent of each other and you can use univariate normal distributions.

1 Like
3

In fact, mu[n] is a vector of size V, and the variance is Theta[n] * D[n] * Theta[n]’, which is V by V. In this case, I couldn’t use univariate normal distributions, right?

4

More detailed and accurate code:

parameters {
    vector[V] mu[N];
    vector<lower=0>[K] sigma_d[N];
    matrix[V, K] Theta[N];

} 
transformed parameters {
  matrix[K, K] D[N]; 
  for (n in 1:N){
    D[n] = diag_matrix(sigma_d[n]);
  }
}
model {
   for(n in 1:N){
   Y[n] ~ multi_normal(mu[n], Theta[n]*D[n]*Theta[n]');
    }
}
5

In this case, you can’t simply use the univariate normal, but without more structure, NUTS is not going to be able to sample from the posterior distribution efficiently enough to make inferences. You can post-multiply Theta[n] by any diagonal matrix and pre-multiply D[n] by the inverse and get the same covariance matrix for Y[n] and hence the same log posterior-kernel. You could just get rid of sigma_d[n] and let it be absorbed into the columns of Theta[n].

6

I am very new to stan programming; do you mean the following? Also, how can I put constraint on Theta so that Theta * Theta’ = I (Identity matrix)?

parameters {
    vector[V] mu[N];
    vector<lower=0>[K] sigma_d[N];
    matrix[V, K] Theta[N];

} 
transformed parameters {
  matrix[V, K] Theta_d[N];
  matrix[K, V] D_Theta[N];
  for (n in 1:N){
    Theta_d[n] = diag_post_multiply(Theta[n], rep_vector(1.0, K));
    D_Theta[n] = diag_pre_multiply(sigma_d[n], Theta[n]');
}
model {
  // priors
  for (n in 1:N){
    mu[n] ~ normal(0, 5); 
    sigma_d[n] ~ cauchy(0,1);

    for(k in 1:K){          
      Theta[n,k] ~ normal(0, 3);
    }
  }

   for(n in 1:N){
     Y[n] ~ multi_normal(mu[n], Theta_d[n] * D_Theta[n]);
    }
}