# How do I use Stan to fit a covariance matrix? [closed]

I'm new to Stan (and bayesian methods in general), so this is likely very simple.

I'm trying to model some multivariate normal data. All I want to know is the covariance matrix generating the data, assuming the data is centered. I've been using this as my template: https://github.com/jrnold/pygments_bugs/blob/53eecb07aa805df7cedc365b7faab641e8fe541d/examples/stan/jaws.stan

But I can't get it to properly fit a known covariance matrix.

I'm using Matlab to simulate some data. Here's the little code snippet:

x = mvnrnd([0,0],[1,1.5;1.5,3],100);

data = struct('N',100,...
'dim',2,...
'x',x);

fit1 = stan('file','singleWishart.stan','data',data,'iter',1000,'chains',4);


And here's the stan code:

data {
int<lower=0> N; // Sample size
int<lower=1> dim; // Number of dimensions
row_vector[dim] x[N]; // Value for each sample on each dimension
}

parameters {
cov_matrix[dim] cov1;
}

model {
vector[dim] zeros; //For now, means are zero
matrix[dim,dim] identity; //Identity for scaling matrix

zeros <- rep_vector(0, dim);
identity <- diag_matrix(rep_vector(1.0,dim));

cov1 ~ inv_wishart(dim, identity);

for (n in 1:N)
x[N] ~ multi_normal_prec(zeros, cov1);
}


I've tried both cov1 ~ inv_wishart and cov1 ~ wishart (I'm not clear on which one I should be using - I thought wishart, but the examples I've seen seem to use inv_wishart), but neither have worked, giving me ridiculous estimates instead of the expected 1, 1.5, 1.5, 3.

## closed as off-topic by gung♦, Tim♦, Sycorax, Xi'an, Nick CoxJul 16 '15 at 15:20

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Tim, Sycorax, Xi'an, Nick Cox
If this question can be reworded to fit the rules in the help center, please edit the question.

• You should declare and define zeros and identity in a transformed data block of a .stan program so that they only get allocated once at the outset rather than every time the posterior function is evaluated (which could be more than a million times with the default settings). – Ben Goodrich Jul 16 '15 at 13:47

The primary reason that your code does not yield the expected answer is that you are using the multi_normal_prec likelihood rather than the multi_normal likelihood. The former expects a precision matrix (the inverse of a covariance matrix) as its second argument, while the latter expects a covariance matrix.
For what it is worth, you should be able to recover the matrix that generated the data if you use multi_normal_prec for the likelihood, use cov1 ~ wishart(dim, identity) for the prior, and then inspect the distribution of the inverse of cov1. But it would be more straightforward, to change the likelihood to multi_normal and leave cov1 ~ inv_wishart(dim,identity)` as is.