Multivariate Multiple Regression (multi_normal + wishart)

[Matt]

Hi!
I've been trying to get a multivariate regression to work in STAN, similar to the model depicted here in JAGS (http://jackman.stanford.edu/classes/BASS/92plain.bug) but with more dependent variables and using Wishart priors. This is the error message I get:

Error in sampler$call_sampler(args_list[[i]]) :

  Error in function stan::prob::multi_normal_cholesky_log(l): i and j must be same. Found i=235, j=3

error occurred during calling the sampler; sampling not done

Below is the code. It runs ok if I replace the multi_normal in the likelihood with a regular normal or t distribution - then I adjust it to loop over the dependent variables and data points, and replace the Omega by a series of independent standard deviations obtained for each variable. I'm not sure what to try next. Any help to get this going would be much appreciated. I've attached a small chunk of the data if it helps.

Matt

MRegMod = '
  data {
    int<lower=1> nData;
    int<lower=1> nVars;
    real x3[nData];
    matrix[nData, nVars] y;
    matrix[nVars, nVars] Imat;
    vector[nVars] b0p;
    vector[nVars] b3p;
  }
  parameters {
    vector[nVars] b0;
    vector[nVars] b3;
    cov_matrix[nVars] Omega;
    cov_matrix[nVars] Omega0;
    cov_matrix[nVars] Omega3;
    real<lower=0,upper=100> sd[nVars];
  }
  transformed parameters {
    matrix[nData, nVars] mu;
    for ( np in 1:nVars){
      for ( i in 1:nData){
        mu[i, np] <- b0[np] + b3[np] * x3[i];
      }
    }
  }

  model {
//Priors
    b0 ~ multi_normal(b0p , Omega0);
    b3 ~ multi_normal(b3p , Omega3);
    Omega0 ~ inv_wishart(nVars+1, Imat);
    Omega3 ~ inv_wishart(nVars+1, Imat);

//Likelihood
    for(i in 1:nData){
      col(y, i) ~ multi_normal(col(mu,i), Omega);
    }
    Omega ~ inv_wishart(nVars+1, Imat);
  }
'

[Bob Carpenter]

I'm combining a users list and dev list post here (see
below).


Matt (user):

Thanks for the careful error report -- it helps us
immensely.

I believe that one of the size tests is failing
(mu must be same size as y, and the rows and columns
of Omega must be the same).

This is puzzling given your model, the declarations
for which are

> matrix[nData, nVars] y;
> matrix[nData, nVars] mu;
> cov_matrix[nVars] Omega;

seem to me to be compatible with the usage:

> col(y, i) ~ multi_normal(col(mu,i), Omega);

Daniel (dev):

We need to get more information into the
test messages that come back. I assume one of
the check_size_match functions is balking here, but
I can't tell which one from the message. Could
we replace "function" with something more error-specific?


Ben (dev):

Could this be in the way we call multi_normal_cholesky
from multi_normal_log? If multi_normal_log doesn't
throw a sizing error, neither should multi_normal_cholesky.
(That is, we want to test size issues as early as possible
for the purpose of user error reporting.)

- Bob

[Marcus Brubaker]

Forgot to include stan-users on this.

Matt, see below for a likely fix to your issue.

---------- Forwarded message ----------
From: Marcus Brubaker <marcus.a...@gmail.com>
Date: Wed, Jan 30, 2013 at 12:23 PM
Subject: Re: [stan-dev] Re: [stan-users] Multivariate Multiple Regression (multi_normal + wishart)
To: stan...@googlegroups.com

On Wed, Jan 30, 2013 at 11:31 AM, Bob Carpenter <ca...@alias-i.com> wrote:

I'm combining a users list and dev list post here (see
below).


Matt (user):

Thanks for the careful error report -- it helps us
immensely.

I believe that one of the size tests is failing
(mu must be same size as y, and the rows and columns
of Omega must be the same).

This is puzzling given your model, the declarations
for which are

>      matrix[nData, nVars] y;
>      matrix[nData, nVars] mu;
>      cov_matrix[nVars] Omega;

seem to me to be compatible with the usage:

>        col(y, i) ~ multi_normal(col(mu,i), Omega);

I believe this should be: 

row(y, i) ~ multi_normal(row(mu,i), Omega);

In the above, col(y,i) and col(mu,i) are vector[nData] but Omega is cov_matrix[nVars].

Cheers,

Marcus

[Daniel Lee]

Our next release will have more comprehensible error messages for matrix size mismatches.

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.
 
 

[Ben Goodrich]

Hi Matt,

On Wednesday, January 30, 2013 6:54:30 AM UTC-5, Matt wrote:

I'm not sure what to try next. Any help to get this going would be much appreciated. I've attached a small chunk of the data if it helps.

What Marcus said should avoid the error message you got, but in general, calling multi_normal() (particularly repeatedly) is a slow process and to a lesser extent, so is inv_wishart(). The reason why this is sad but true is that these functions take a Cholesky factor and an inverse of their inputs every time the log-posterior is evaluated.

Thus, you would be better off

1) Calculating the Cholesky factor of Omega
2) Inverting it (i.e. doing a triangular solve with the identity matrix) and transposing

These two steps only have to be done once in the model block. Then for each variable ...

3) Premultiplying each column of y - mu by the result of 2)
4) Subtracting half the sum of squares of 3) from the log-posterior

There are some details in the Stan manual about this, particularly in the chapter entitled Optimizing Stan Code.

Ben

 

[Matt]

Thanks Marcus, I tried as you suggested -
 

   for(i in 1:nData){

      row(y, i) ~ multi_normal(row(mu,i), Omega);
    }

and also replacing "nData" with nVars, but I get this error message:

TRANSLATING MODEL 'albMRegMod' FROM Stan CODE TO C++ CODE NOW.
Error in stanc(file = file, model_code = model_code, model_name = model_name,  :
  failed to parse Stan model 'albMRegMod' with error message:
EXPECTATION FAILURE LOCATION: file=input; line=37, column=7

      row(y, i) ~ multi_normal(row(mu,i), Omega);

      ^-- here


DIAGNOSTIC(S) FROM PARSER:
no matches for function name="multi_normal_log"
    arg 0 type=row vector
    arg 1 type=row vector
    arg 2 type=matrix
unknown distribution=multi_normal
Parser expecting: "}"

Which makes me think I'm doing something more fundamentally wrong.


And, thanks for the tips Ben. When I get this working I'll definitely have a look at optimising.

Matt

[Bob Carpenter]

Sorry -- those have to be vectors the way multi_normal
is currently defined.

You can look in the index of the manual to see the
signatures for which each function is defined.
You'll see that multi_normal is only defined for
vectors, not row vectors. So you can either
transpose (easy) or reorganize your data structures (more
efficient).

So

row(y, i) ~ multi_normal(row(mu,i), Omega);

should transpose the rows into (column) vectors:

row(y,i)' ~ multi_normal(row(mu,i)', Omega);

The more efficient approach would be to redefine
y to be an array of vectors, as in:

vector[???] y[nData];
vector[???] mu[nData];

y[i] ~ multi_normal(mu[i], Omega);

The advantage of this latter approach is that once we've
vectorized the multi_normal, you'll be able to write it
all without the loop as:

y ~ multi_normal(mu,Omega);

- Bob

[Matt]

Thanks Bob, that works.

[Mauricio Garnier-Villarreal]

Hi Matt

could you share your stan code? I am running a similar model and it would be helpfull to see your working code

thanks

[vera]

Hello Matt, I am also interested in multiple outcomes...Could you please share your stan code...

[Matthew Albrecht]

Sorry guys, I didn't catch these, I only check the stan list sporadically.

Let me know if you're still interested and I'll see if I can dig up the code from where I left off on these models. I've also implemented something similar to http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3031784/ if you're interested.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[vera]

Hello Matt, Yes I am still interested to see your code...thanks