4
$\begingroup$

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.

$\endgroup$

closed as off-topic by gung, Tim, Sycorax, Xi'an, Nick Cox Jul 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.

  • 2
    $\begingroup$ 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). $\endgroup$ – Ben Goodrich Jul 16 '15 at 13:47
6
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ All that said, I do not ever use the Wishart or inverse-Wishart distributions with Stan for reasons given at groups.google.com/d/msg/stan-users/fhghq6Rxbqs/-dG0apK0NV4J . If you must utilize the inverse-Wishart, to get decent results from Stan I would integrate it out of the posterior as described at en.wikipedia.org/wiki/… . $\endgroup$ – Ben Goodrich Jul 16 '15 at 13:52
  • 1
    $\begingroup$ For the benefit of OP, I'd like to add that using a parameterization that separates the variances from correlations is superior in Stan because it allows each to move separately, which is important because Stan cares about the geometry of the posterior. $\endgroup$ – Sycorax Jul 16 '15 at 14:05

Not the answer you're looking for? Browse other questions tagged or ask your own question.