Averaging generated covariance matrices

General
1

Hello all, I am running a code in which one of the parameters is a covariance matrix, Stan gives the generated samples in the form of an array with dimensions [n,p,q], where n number of generated observations and p=q is the dimension of the covariance matrix. Using the extract function

out_par<-extract(out)
covm<-out_par$Sigma

Now I need to obtain the average covariance matrix, as the posterior estimate of the Sigma, I am not sure how to do it, I could average every single cell of the matrix, but I don’t think that’s efficient. thanks in advance.

2

To what purpose? If as an intermediate step in the computation of another quantity, best to actually compute said quantity in each sample from the posterior rather than trying to average the posterior then computing the quantity.

1 Like
3

The purpose/interest is covariance matrix itself, since I will use it to compute the correlation matrix, which is simple matrix operation of the covariance matrix.

4

Right, so you should compute the correlation matrix in each sample of the posterior, giving you n correlation matrices that you can then summarize however you’d like (ex. looking at the posterior for a given cell).

1 Like
5

Hi Mike, isn’t easier to just compute the correlation matrix afterwards, instead of adding an extra line to the code? The code is simulating a covariance matrix \Sigma=(\sigma_{ij}) (i,j=1,\cdots,p) which the prior is a Wishart, if I have 2000 mcmc steps (excluding the burn-in’s) I will have 2000 matrices, I just want an average of them, then I can compute the correlation matrix C=D^{-1/2}\Sigma D^{-1/2}, where D=diag(\Sigma) and D^{1/2}=(1/\sigma_{ii}).

6

Given parameter \theta and transform f(\theta), it is only under select circumstances that f(E(\theta)) (function applied to the mean posterior) equals E(f(\theta)) (mean of the function as applied to each sample of the posterior, the thing we’re typically interested in for inference). Instead of keeping track of when those circumstances hold (I’m not very math-y), I think it’s a good default to always compute the latter rather than the former.

2 Likes
7

Right, however I still have to average the 2000 generated correlation matrices, I am not sure how to do it? cheers.

8

Once you have the correlation matrices, what you do next again depends on what you want to ultimately achieve. If you are interested in a given element of the matrix, you can look at it’s full posterior or summaries (mean, quantiles, etc) thereof. If you want to use that correlation matrix as an input to another computation (ex. get the eigenvalues for factor analysis, etc), then you perform that computation in each sample from the posterior, yielding one or more new quantities that you now have a full/proper posterior for.

2 Likes
9

Hi, You’re right we should include the correlation matrix in the model, but I am not sure how to tell that to Stan, the model is like

 model {
 mu ~ normal(0,100);
  for(i in 1:N){
   X[i,] ~ multi_normal(mu,Sigma);
  }
    Sigma ~ inv_wishart(k,R);
 }

k and R are given.

The correlation matrix C can be obtained by C=D^{-1/2}\Sigma D^{-1/2}, where D=diag(\Sigma) and D^{1/2}=(1/\sigma_{ii}). But I am not sure how to code it, Stan does not seem to have a function to extract the diagonal of a matrix. I would appreciate any help. thanks.

10

Here is a function I coded years ago that I believe works. Don’t ask me what I was thinking re the naming.

 matrix cov2cor(matrix V){ 
   int p = rows(V);
   vector[p] Is = sqrt(1.0 ./ diagonal(V));
   matrix[p,p] Ism;
   matrix[p,p] r;
   for(i in 1:p) Ism[i,] = Is';
    r =  Ism .* V .* Ism';
   for(i in 1:p) r[i,i]=1;
   return r;
 }
4 Likes
11

Same thing but faster

 matrix cov2cor(matrix V) {
   int p = rows(V);
   vector[p] Is = inv_sqrt(diagonal(V));
   
   return quad_form_diag(V, Is);
 }
5 Likes
12

Thanks you all, it worked fine.

1 Like