Trouble with Multivariate Regression

[Nathan Lemoine]

I'm having trouble running a toy examples for multivariate regression. This stan code is:

stanMod = """

data{

    int<lower = 0> N;

    int<lower = 0> K;

    vector[K] Y[N];

    vector[N] D;

    cov_matrix[K] R;

}

parameters{

    vector[K] a;

    vector[K] b;

    cov_matrix[K] Omega;

}

transformed parameters{

    vector[K] mu[N];

    for(k in 1:K){

        for(i in 1:N){

            mu[i, k] <- a[k] + b[k] * D[i];

        }

    }

}

model{

    for(i in 1:N){

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

    }

    Omega ~ inv_wishart(K+1, R);

}

"""

I get this warning:

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:

Error in function stan::prob::inv_wishart_log(N4stan5agrad3varE): WW is not positive definite (7.3809879583185371e+71:0).

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

I'm not sure what I'm doing wrong, if this is the right way to parameterize it. This model works in JAGS so I thought I'd try porting it over to STAN. I've read other threads that suggest the problem is in the Omega (in particular the inv_wishart coming too close to not being positive definite for the covariance matrix). I've tried numerous things, such as changing Omega from a covariance matrix to a regular matrix (relaxing the requirement of positive definiteness I think). The model runs, but the chains do not converge at all. I've tried removing the multivariate response and changing the model to:

for(i in 1:N){

        for(k in 1:K){

            Y[i, k] ~ normal(mu[i], stdev[k]);

        }

    }

I still get an error message:

Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:

Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is -3.364634555144618:0, but must be > 0!

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

but the model runs and the chains do converge. I'm not crazy about this because a) I still get error messages, b) I'm not modeling the covariance between response variables, which is preferable, and c) the parameter estimates are just too similar (almost identical) for a[1]:a[2] and b[1]:b[2], so I'm not sure I trust it.

I've tried reparameterizing to the Cholesky decomposition as in the manual, but most examples I've found are for multivariate normal of parameters, not a response, and I'm pretty sure I did it wrong.

Am I screwing something up?

[Bob Carpenter]

Hopefully Ben can jump in with some insight here.
Some comments below.

On 1/21/14, 3:42 AM, Nathan Lemoine wrote:
> I'm having trouble running a toy examples for multivariate regression. This stan code is:
>
> stanMod = """
> data{
> int<lower = 0> N;
> int<lower = 0> K;
>
> vector[K] Y[N];
> vector[N] D;
>
> cov_matrix[K] R;
> }
>
> parameters{
> vector[K] a;
> vector[K] b;
>
> cov_matrix[K] Omega;
> }
>
> transformed parameters{
> vector[K] mu[N];
>
> for(k in 1:K){
> for(i in 1:N){
> mu[i, k] <- a[k] + b[k] * D[i];
> }
> }
> }
>
> model{
> for(i in 1:N){
> Y[i] ~ multi_normal(mu[i], Omega);
> }
>
> Omega ~ inv_wishart(K+1, R);
> }
> """

We're no longer very keen on the inverse Wishart priors.
See this paper by Ben Goodrich et al.:

http://www.stat.columbia.edu/~gelman/research/unpublished/Visualization.pdf

But the model should still work.

>
> I get this warning:
>
> /Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:/
> /Error in function stan::prob::inv_wishart_log(N4stan5agrad3varE): WW is not positive definite (7.3809879583185371e+71:0)./
> /If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the
> sampler is fine,/
> /but if this warning occurs often then your model may be either severely ill-conditioned or misspecified./

For what it's worth, that's *not* an error message. We just
don't know how to make that point any clearer to users. We don't
want to get rid of it as it's helpful, but we may have to in order
to stop freaking our users out about it.

The dump of the value is a mess here --- we really need to fix
that. But that 7e71 seems to indicate that you're getting
huge values in the covariance matrix.

The key issue is: does the error persist in most iterations? If
so, there's a problem.

I think we need to lower our tolerance for positive definiteness
somehow. We're using the Eigen built-ins to do this, and they're
very picky about what counts as positive definite.

> I'm not sure what I'm doing wrong, if this is the right way to parameterize it. This model works in JAGS so I thought
> I'd try porting it over to STAN. I've read other threads that suggest the problem is in the Omega (in particular the
> inv_wishart coming too close to not being positive definite for the covariance matrix). I've tried numerous things, such
> as changing Omega from a covariance matrix to a regular matrix (relaxing the requirement of positive definiteness I
> think).

No, you need to declare Omega as a covariance matrix. In general,
the constraints on a parameter need to match its support in the model.

This error's coming from the multivariate normal distribution, not
from the parameter itself.

The model runs, but the chains do not converge at all. I've tried removing the multivariate response and
> changing the model to:
>
> for(i in 1:N){
> for(k in 1:K){
> Y[i, k] ~ normal(mu[i], stdev[k]);
> }
> }
>
> I still get an error message:
>

> /Informational Message: The current Metropolis proposal is about to be rejected becuase of the following issue:/
> /Error in function stan::prob::normal_log(N4stan5agrad3varE): Scale parameter is -3.364634555144618:0, but must be > 0!/
> /If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the
> sampler is fine,/
> /but if this warning occurs often then your model may be either severely ill-conditioned or misspecified./

>
> but the model runs and the chains do converge. I'm not crazy about this because a) I still get error messages,

See above.

b) I'm
> not modeling the covariance between response variables, which is preferable,

Right --- this one's the deal breaker.

> and c) the parameter estimates are just too
> similar (almost identical) for a[1]:a[2] and b[1]:b[2], so I'm not sure I trust it.

I'd recommend simulating data from the model and seeing
if you can recover the parameters.

> I've tried reparameterizing to the Cholesky decomposition as in the manual, but most examples I've found are for
> multivariate normal of parameters, not a response, and I'm pretty sure I did it wrong.

You can use the Cholesky decomposition on parameters.
It's just going in as an argument to the multivariate normal.

> Am I screwing something up?

Not that I can see from the model, but I don't know what your data
are like. If the covariance really does involve terms on the scale of
7e71, there may be a problem with the data.

- Bob

[Bob Carpenter]

OK, I was too curious to see whether this model really
works or not, so I did the simulation myself.

I've attached a data simulator in R and a reformatted model that work
just fine. For example, running one chain (sorry, Andrew!) of 200 iterations
gives me:

Inference for Stan model: mr.
1 chains, each with iter=200; warmup=100; thin=1;
post-warmup draws per chain=100, total post-warmup draws=100.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
a[1] 1.0 0.0 0.1 0.8 0.9 1.0 1.1 1.3 25 1.1
a[2] -0.4 0.0 0.1 -0.6 -0.5 -0.4 -0.3 -0.3 26 1.1
a[3] -1.1 0.0 0.1 -1.3 -1.2 -1.1 -1.0 -0.9 26 1.1
a[4] -1.5 0.0 0.1 -1.8 -1.6 -1.5 -1.4 -1.3 33 1.0
b[1] 0.6 0.0 0.1 0.3 0.5 0.6 0.7 0.9 24 1.0
b[2] 0.0 0.0 0.1 -0.2 -0.1 0.0 0.0 0.2 25 1.0
b[3] -0.3 0.0 0.1 -0.5 -0.3 -0.3 -0.2 0.0 25 1.0
b[4] -0.7 0.0 0.1 -1.0 -0.8 -0.7 -0.6 -0.5 22 1.1
Omega[1,1] 1.6 0.0 0.2 1.2 1.4 1.6 1.7 2.0 52 1.0
Omega[1,2] -1.1 0.0 0.1 -1.4 -1.2 -1.1 -1.0 -0.9 47 1.0
Omega[1,3] 1.2 0.0 0.2 0.9 1.1 1.2 1.3 1.5 48 1.0
Omega[1,4] -0.9 0.0 0.2 -1.2 -1.0 -0.9 -0.7 -0.6 66 1.0
Omega[2,1] -1.1 0.0 0.1 -1.4 -1.2 -1.1 -1.0 -0.9 47 1.0
Omega[2,2] 0.8 0.0 0.1 0.6 0.7 0.8 0.8 1.0 45 1.0
Omega[2,3] -0.9 0.0 0.1 -1.1 -1.0 -0.9 -0.8 -0.7 46 1.0
Omega[2,4] 0.6 0.0 0.1 0.4 0.5 0.6 0.6 0.8 58 1.0
Omega[3,1] 1.2 0.0 0.2 0.9 1.1 1.2 1.3 1.5 48 1.0
Omega[3,2] -0.9 0.0 0.1 -1.1 -1.0 -0.9 -0.8 -0.7 46 1.0
Omega[3,3] 1.0 0.0 0.1 0.8 0.9 1.0 1.1 1.3 48 1.0
Omega[3,4] -0.6 0.0 0.1 -0.9 -0.7 -0.6 -0.5 -0.4 58 1.0
Omega[4,1] -0.9 0.0 0.2 -1.2 -1.0 -0.9 -0.7 -0.6 66 1.0
Omega[4,2] 0.6 0.0 0.1 0.4 0.5 0.6 0.6 0.8 58 1.0
Omega[4,3] -0.6 0.0 0.1 -0.9 -0.7 -0.6 -0.5 -0.4 58 1.0
Omega[4,4] 1.5 0.0 0.2 1.2 1.4 1.6 1.7 1.9 81 1.0


simulated values:

> a
[1] 0.901863802116 -0.354475350437 -1.187679480554 -1.523551737507

> b
[1] 0.57476634644000 -0.00932367605421 -0.29102771809994 -0.68947032081146

> Omega
[,1] [,2] [,3] [,4]
[1,] 1.852199167310 -1.259972460476 1.383790043310 -0.686849949405
[2,] -1.259972460476 0.880673049025 -0.978498101009 0.442150515322
[3,] 1.383790043310 -0.978498101009 1.098653245960 -0.469949555671
[4,] -0.686849949405 0.442150515322 -0.469949555671 1.354169601243


How many dimensions (K) were you using?

I am having problem with RStan's output --- the overlay
update comes before the warning messages, so I can't follow progress
because it scrolls off the screen [Mac OS X Mavericks, RStan 2.1.0.
And then R goes all wonky. (Let's just say I'm not the biggest fan
of R.)

- Bob

[Nathan Lemoine]

Thanks for the help. I thought it was an error message suggesting my results are completely invalid. I went back and re-ran the model and, though I still got the warning, I followed through and the parameter estimates do match those of JAGS, so it appears to be working.

I've attached the R script and datafile so you can see what kind of data I'm using to get this warning. On some level, it does make sense that I've getting large covariances because the two response variables are highly correlated.

I'll keep playing with it to see if I run into more issues or if the model is working fine despite the warning.

[Bob Carpenter]

We expect to get these warnings, especially with highly
constrained models. We can't seem to formulate text that
doesn't lead to these kinds of exchanges on the mailing list,
so I think it's probably best if we just get rid of them.

As warning message itself says, if they do not persist
they are not a problem. In particular they do not affect
the detailed balance of the sampler --- they are just
indicators that there's a Metropolis rejection due to
constraint violations, which is itself due to discretization
of the Hamiltonian dynamics and limited floating-point
precision.

I think the large covariances are coming from the initializations.
I saw them in the simulation, too. Even with a perfect
correlation, at convergencye, you'd only expect the posterior covariance
to be somewhere near the larger variance of the two variables.

- Bob

> --
> 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.

[Ken Arnold]

So if I understand correctly, those warnings are of no concern if they happen rarely, but helpful to localize problems with the model if they happen frequently?

If so, how about keeping count of the number of errors, and only printing the details if that count exceeds a threshold? Or, perhaps, bin the messages by type and report a summary of the types of errors encountered?

-Ken

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

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

--
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+unsubscribe@googlegroups.com.

[John Hall]

Bob,
I've been having this warning since transitioning to the latest version of Stan. Maybe not necessarily removing it, but having different levels of warnings displayed might be helpful (like Matlab's fmincon has a display variable that you can adjust to display increasing levels of details about an optimization).

Anyway, I had been using the following code (based on a thread from January of last year):

//For testing, I had set S_prior to the identity matrix, V_prior to S_prior times 100, and nu_prior to N.
data {
  int<lower=0> N; //# of variables
  int<lower=0> T; //# of observations
  vector[N] y[T]; //data, formatted as NxT
  vector[N] mu_prior;  //prior on mean
  cov_matrix[N] V_prior; //prior uncertainty on mean
  cov_matrix[N] S_prior; //prior on covariance
  int<lower=0> nu_prior; //prior degrees of freedom for covariance
}
transformed data {
  matrix[N, N] L_V_prior; //Haven't had any luck getting cholesky_factor_cov[N] working
  matrix[N, N] L_S_prior;
  matrix[N, N] inv_L_S_prior;
  int<lower=0> nu_post;
 
  L_V_prior <- cholesky_decompose(V_prior);
  L_S_prior <- cholesky_decompose(S_prior);
  inv_L_S_prior <- inverse(L_S_prior);
  nu_post <- T + nu_prior;
}
parameters {
  vector[N] mu;
  cov_matrix[N] sigma;
}
transformed parameters {
  matrix[N, N] L_sigma;
  L_sigma <- cholesky_decompose(sigma);
}
model {
  mu ~ multi_normal_cholesky(mu_prior, L_V_prior);
  sigma ~ inv_wishart(nu_post, inv_L_S_prior' * inv_L_S_prior);
  for (t in 1:T) {
    y[t] ~ multi_normal_cholesky(mu, L_sigma);
  }
}

You suggest not using the inverse Wishart, but still using the Cholesky. I was able to change the end to

parameters {
  vector[N] mu;
  cov_matrix[N] inv_sigma;
}
transformed parameters {
  matrix[N, N] L_inv_sigma;
  L_inv_sigma <- cholesky_decompose(inv_sigma);
}
model {
  mu ~ multi_normal_cholesky(mu_prior, L_V_prior);
  inv_sigma ~ wishart(nu_post, L_S_prior * L_S_prior');
  for (t in 1:T) {
    y[t] ~ multi_normal_prec(mu, L_inv_sigma * L_inv_sigma');
  }
}

and get approximately what I was getting with the original code (I had also tried with just inv_sigma in the multi_normal_prec statement). However, this version is pretty slow since it's using multi_normal_prec instead of multi_normal_cholesky. I wasn't able to get multi_normal_cholesky working when converting inv_sigma back to sigma and then L_sigma.

Any tips? Or general best practices advice for this type of model?

[Bob Carpenter]

On 1/23/14, 9:40 PM, John Hall wrote:
> Bob,
> I've been having this warning since transitioning to the latest version of Stan. Maybe not necessarily removing it, but
> having different levels of warnings displayed might be helpful (like Matlab's fmincon has a display variable that you
> can adjust to display increasing levels of details about an optimization).

I was taken enough with this idea to create an issue:

https://github.com/stan-dev/stan/issues/512

Congratuatlions, it's the 2^9-th issue :-)

We keep improving the way we handle this, but we're not quite
to where we want to be yet.

I'd recommend scaling a correlation matrix. Give the correlation_matrix
an lkj_correlation prior, see the doc. Then you can use the same
kind of priors you'd use for any scale parameter for the scales.

We're about to roll out 2.2.0 (I think we're bypassing 2.1.1), and
it has diag_pre_multiply and diag_post_multiply
functions which can be used to neatly scale.

So it'd be:

parameters {
corr_matrix[N] Omega;
vector<lower=0> sigma;
}
transformed parameters {
cholesky_factor_cov[N] L;
L <- cholesky_decompose(Omega);
for (i in 1:N)
for (j in 1:i)
L[i,j] <- L[i,j] * sigma[i];
}
model {
Omega ~ lkj_corr(nu); // nu controls concentration around unit
sigma ~ ... whatever you like for scales ...;
...

I'll drop a section on this into the next manual
release.

We're going to introduce a correlation matrix Cholesky
factor type which should be more efficient, but have the
same effect as the above.

For a discussion and a link to Ben's paper, see

http://andrewgelman.com/2012/08/30/visualizing-distributions-of-covariance-matrices/

- Bob

> stan-users+...@googlegroups.com <javascript:>.
> > For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.

[John Hall]

Thanks for the quick reply. I remember something like that in the manual; I'll give that a try when I have some time. The only thing I'm not sure on with this approach is if there will ever be a meaningful correlation between the sigmas when you do it that way (assume for instance you use a inverse gamma). If you do the estimation with the inverse Wishart, it might be that if you plot the posterior for two variances they are correlated. I have actually tried something similar in a Gibbs sampling approach but since each inverse gamma draw was independent it took out the correlations in the variances.

Also, if you're going to recommend working with the correlation matrices, it seems like a next step would be to consider working with copulas. If the multivariate dependence were set up with copulas (even a normal copula), then it presumably would be easier to switch from normal marginals to a t marginal, or assume a constant dependence with stochastic volatility.

On a separate note (not to bog you guys down with too many issues), while the extract function can allow you to only select certain outputs, my understanding is that the stan function simulates every parameter and transformed parameter (and generated quantities) from a .stan file. However, I may not always need to output all the transformed parameters (and generated quantities). In the above example, I may not really need to output the simulated stochastic volatility of each series in each point in time. As the number of variables or length of the series increases, it may cause memory issues on less powerful machines. Similarly, it might be convenient if you could input the output of a previous stan function call and then output whichever other transformed parameters or generated quantities you may need. In that example, it would be like fitting the model with one stan call and outputting only some of the parameters, then calling the stan function again to calculate the posterior predictive given the previous call. This way you wouldn't calculate the generated quantities unless you need them and you don't have to rely on separate stan files (one with generated quantities, one without) when testing.

[Bob Carpenter]

On 1/24/14, 2:05 AM, John Hall wrote:
> Thanks for the quick reply. I remember something like that in the manual; I'll give that a try when I have some time.
> The only thing I'm not sure on with this approach is if there will ever be a meaningful correlation between the sigmas
> when you do it that way (assume for instance you use a inverse gamma). If you do the estimation with the inverse
> Wishart, it might be that if you plot the posterior for two variances they are correlated. I have actually tried
> something similar in a Gibbs sampling approach but since each inverse gamma draw was independent it took out the
> correlations in the variances.

That's really two issues.

1. The LKJ prior still allows correlations. You can
make it uniform on the correlation matrix. Or you
can create a prior with more probability mass around
the identity correlation matrix.

2. Gibbs sampling has a lot of trouble with the multivariate
priors because of the correlations --- it basically takes
forever to converge to the correct posterior. You'll find HMC
does much better with this.

> Also, if you're going to recommend working with the correlation matrices, it seems like a next step would be to consider
> working with copulas. If the multivariate dependence were set up with copulas (even a normal copula), then it presumably
> would be easier to switch from normal marginals to a t marginal, or assume a constant dependence with stochastic volatility.

I've heard of copulas, but don't understand them. Is there
a simple reference for the constraints on parameters and densities
involved? If it's a fixed number of parameters --- that is, it's
not a non-parametric method like Dirichlet processes --- then we
should be able to code it in Stan. Either in the Stan modeling
language or in the underlying C++ --- the latter is more efficient
and stable.

> On a separate note (not to bog you guys down with too many issues), while the extract function can allow you to only
> select certain outputs, my understanding is that the stan function simulates every parameter and transformed parameter
> (and generated quantities) from a .stan file. However, I may not always need to output all the transformed parameters
> (and generated quantities). In the above example, I may not really need to output the simulated stochastic volatility of
> each series in each point in time. As the number of variables or length of the series increases, it may cause memory
> issues on less powerful machines.

You can write the output to files so you won't blow out during sampling,
but then the posterior analysis can be a problem.

That's already on our to-do list. It's not an immediate priority,
but we'll get to it eventually.

It's similar to the BUGS/JAGS monitor specification.

> Similarly, it might be convenient if you could input the output of a previous stan
> function call and then output whichever other transformed parameters or generated quantities you may need.

That's another very good idea. It would force us to modularize the generated
model code a bit to support it, but the underlying code generator
is already modularized that way.

https://github.com/stan-dev/stan/issues/513

We've also been thinking about issues like that in order to do posterior
predictive checks more seamlessly.

- Bob

[John Hall]

Here are some references (the first is a more gentle introduction, the file is below the cached heading):
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.160.8266
Nelson's "An Introduction to Copulas" is a popular reference that you could find at an academic library
http://arxiv.org/abs/1112.4204

Copulas are the multivariate dependence in a joint distribution. Sklar's theorem is the most important result, which says that any joint distribution can be expressed as some univariate distributions times the copula. The benefit of this approach is that you can model the univariate marginal distributions separately from the multivariate dependence. As such, copula estimation often requires the inverse pdf for the marginals, which might be one barrier to implementing them in Stan presently (I didn't see inverse pdfs in the manual).

I used student t marginals and a normal copula as a good example above of why copulas are helpful. Basically, if you were to instead fit a multivariate student t distribution, then you effectively assume the degrees of freedom parameter is the same across every variable. The student t + normal copula would often be fit by first fitting each of the variables with a student t distribution with each having its own degrees of freedom. Then you would create a new multivariate variable applying the appropriate student t inverse pdf to each variable. This puts each variable between 0 and 1. To fit the normal copula to this, I might apply a standard normal cdf to that and estimate the resulting correlation matrix.

[Bob Carpenter]

Thanks for the quick and clear explanation of what
they're trying to do.

You're right that we don't have the inverse cumulative distribution
functions (CDFs). We also don't have any way to first fit just marginals
and use those marginal fits to fit correlations, but I'll take a look
at the Michael Smith paper on Bayesian approaches to copulas (the
arXiv link).

As an aside, we're using the Joe et al. vine approach in our transforms
for modeling correlation matrices; see the chapter of the manual on
transforms and Jacobians near the end of the manual.

- Bob

[Ravi Sastry]

John,

I'm trying to implement precisely the model you described: multivariate data with t marginals and a normal copula. Have you implemented this model? Any advice on the best way to proceed?

I'm using a scale mixture representation, so each observation is augmented with a gamma distributed scaling parameter. I then construct the covariance matrix Sigma[t] for each time period, and sample using

y[t] ~ multi_normal(mu, Sigma[t])

(Everything is actually implemented in terms of Cholesky factors, and I sample using multi_normal_cholesky.) Overall, I'd say the results are pretty good for recovering the scale and location parameters, much less so for the degrees of freedom parameters -- even with a highly informative sample for which MLE works quickly and accurately.

Any thoughts or advice? Thanks.

Ravi

[Bob Carpenter]

On Aug 13, 2014, at 6:46 PM, Ravi Sastry <ravis...@gmail.com> wrote:

> John,
>
> I'm trying to implement precisely the model you described: multivariate data with t marginals and a normal copula. Have you implemented this model? Any advice on the best way to proceed?
>
> I'm using a scale mixture representation, so each observation is augmented with a gamma distributed scaling parameter. I then construct the covariance matrix Sigma[t] for each time period, and sample using
>
> y[t] ~ multi_normal(mu, Sigma[t])
>
> (Everything is actually implemented in terms of Cholesky factors, and I sample using multi_normal_cholesky.) Overall, I'd say the results are pretty good for recovering the scale and location parameters, much less so for the degrees of freedom parameters -- even with a highly informative sample for which MLE works quickly and accurately.

How is the effective sample size for degrees of freedom? Is the simulated
value in the posterior 95% interval?

If both of these tests look OK, then the question is how close you expect
the posterior mean to be to the posterior mode. Stan returns the posterior
mean and median in its summaries of MCMC output, not the posterior mode.

- Bob

[John Hall]

Ravi,
I don't think I ever got this to work. Just kinda gave up on it. Would be interested if you get it to work.

[Ravi Sastry]

I generated a 1000x2 data matrix where the respective degrees of freedom were 6 and 15. Correlation was 0.25. Here are the relevant output lines:

            Mean     MCSE   StdDev        5%       50%       95%  N_Eff  N_Eff/s    R_hat

nu[1]    4.0e+00  2.4e-04  4.0e-03   4.0e+00   4.0e+00   4.0e+00    268  3.2e-01  1.0e+00

nu[2]    4.0e+00  2.4e-04  3.8e-03   4.0e+00   4.0e+00   4.0e+00    260  3.1e-01  1.0e+00

R_hat is 1 in both cases, so it seems to think it's converged. But in both cases, it's decided that nu=4. That happens to be the lower bound on my prior for nu. This suggests a coding error on my part, but I can't find it. (A modified version of the model for univariate t data properly recovers the degrees of freedom using the same scale mixture representation.)

Ravi

[Bob Carpenter]

OK --- that's definitely not a problem with too much posterior
uncertainty.

An R_hat != 1 indicates non-convergence, but R_hat = 1 does not
guarantee convergence. Has the lp__ parameter converged?

Feel free to send your code and we can take a look to see if
there's anything obviously wrong in the way Stan's being used.

- Bob

> --
> 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.

[Ravi Sastry]

R_hat for lp__ is 2.1. So that hasn't converged... I'll try running with more samples, but the run I'm reporting had 10,000 samples.

My code is attached to my original post. I certainly can't rule out coding error... I'd be thrilled if you found a mistake. But I was mostly hoping you could comment on whether this is the best approach to implement this particular model.

Ravi

[John Hall]

Ravi,
Just to give you a little more information. I think when I had first tried this, I didn't have much understanding of Hamiltonian Monte Carlo. I still don't (compared to Gibbs Sampling), but I'm slightly more familiar with it than I used to be. The good thing about Gibbs sampling is that you can pretty much do whatever you want. Like I could throw in a PCA step to a Gibbs sampler and the only problem is that it slows everything down. You (likely) couldn't use PCA within the context of Stan because analytic formulas for PCA are ridiculous (I've tried symbolic math toolboxes for even 5x5 matrices and the formula look like just a sea of variables).

In Stan, you need to be able to provide analytic information for the log likelihood (or at least something proportional to it). The problem for a copula approach is that the copula log likelihood will have some quantile (inverse cdf) functions, depending on what the marginals are, that are not easy to calculate derivatives for. For instance, suppose you fit a t distribution to each univariate distribution. To fit the normal copula, you would then need to apply the cdf function for the relevant t to each univariate distribution and then apply the normal inverse cdf function to get normally distributed variables, which you could then estimate the correlation matrix from. Last time I checked, Stan has cdfs but not inverse cdfs, so at best you could get the uniform marginals (assuming Stan can handle taking derivatives of them, which I haven't tried). You'd then need some other approach to convert the uniform marginals to values, but even if you could write a user-defined function to do it (which you couldn't before) Stan would still need the derivatives of the functions for the HMC to work properly, I think at least.

[Luc Coffeng]

My 50 cents on inverse normal CDF: have you considered approximating the (inverse) erf with elementary functions? It would be ugly, but probably possible to differentiate.

[Herra Huu]

I'm not sure if this fix everything, but at least this part of the Ravi's model is problematic:

    vector<lower=0>[J] sigma ;      // scale of multivariate vector
    vector<lower=0>[J] nu ;  
    ...
    sigma[j] ~ uniform(0,100) ;
    nu[j] ~ uniform(4,30) ;

The problem is that Stan is generating parameters outside of the support of parameters (that is, from interval (0, inf) and not (0,100) and (4,30)). Easy fix would be to declare parameters as (impling the same uniform prior as in the original model, no need to explicitly write them down):

    vector<lower=0,upper=100>[J] sigma ;      // scale of multivariate vector
    vector<lower=4,upper=30>[J] nu ;  

[Ravi Sastry]

John and Luc,

I was (of course) hoping to let Stan generate the full joint log density and gradient on its own. But if my current approach doesn't work out, and since there is no support yet for copula distributions, I will need to explicitly specify the log density as you suggest. Those approximations for erf() may prove essential. 

Herra,

I'll make the change to my declared bounds to see if it fixes things. I hope so! But my simpler model with univariate t data using the scale mixture representation had the same bounds and priors on the relevant parameters, and that converged without issue. Do you know if the Stan documentation warns against specifying a prior whose support doesn't exactly match the declared bounds? (As you see from my code, I assumed that having prior support a subset of the declared bounds would be OK.)

Thank you all for your time and help!

Ravi

[Bob Carpenter]

Yes, the doc says this all over the place.

- Bob

You can only do it the other way around.

For instance, a half-Cauchy prior can be defined by

parameters {
real<lower=0> sigma;
...
}
model
sigma ~ cauchy(0,2.5);
...

You only need an explicit truncation when either one of the bounds or
one of the parameters of the distribution is a parameter --- otherwise, you
can just drop it, as above, where it would be T[0,].

But if you do this:

parameters {
real sigma;
...
}
model {
sigma ~ lognormal(0,1);
...

you'll almost certainly run into problems sampling when sigma goes below 0
either doing initialization or during the Hamiltonian dynamics. That causes
a rejection, and if it's an init, that's fatal.

- Bob

[Ravi Sastry]

After making the appropriate changes to my parameter declarations and priors, my model is working. I've tested it on simulated data up to 20 dimensions.

Code is attached if you're interested. Any suggestions to help express the model more efficiently would be appreciated.

Ravi

[Bob Carpenter]

1. Rather than this:

for(j in 1:J) {
for(t in 1:T)
lambda[t,j] ~ gamma(nu[j]/2,nu[j]/2) ;

you can use:

for (t in 1:T)

lambda[t] ~ gamma(nu / 2, nu / 2);

Even better would be to not divide nu by 2 twice, but to do it once
and assign it to a temporary.

2. Replace

for (t in 1:T) {

for (j in 1:J)
lambdainvsqrt[t][j] <- inv_sqrt(lambda[t][j]) ;
LSigma[t] <- LRho * diag_matrix(sigma) * diag_matrix(lambdainvsqrt[t]) ;
}

with

{
matrix[J,J] LRho_times_diag_matrix_sigma;
LRho_sigma <- diag_post_multiply(LRho, sigma);

for (t in 1:T) {

for (j in 1:J)
lambdainvsqrt[t][j] <- inv_sqrt(lambda[t][j]) ;
LSigma[t] <- diag_post_multiply(LRho_sigma, lambda_inv_sqrt[t]);
}

It may be even faster in a loop due to the sparsity. In just
two dimensions, you can just write

LSigma[t,1,1] <- ...
LSigma[t,2,2] <- ...
LSigma[t,2,1] <- ...
LSigma[t,1,2] <- 0; // lower triangular

with ... replaced by the reuslt of all the above matrix multiplication.
Make sure to save intermediate computations.

3. If you don't need to save lambda_inv_sqrt[t], make it a local in the
same scope.



4. This isn't a speed issue per se, but Andrew's been very down on the constrained
uniform priors (though we still have some in our manual):

vector<lower=0,upper=100>[J] sigma ; // scale of multivariate vector

vector<lower=2,upper=30>[J] nu ; // degrees of freedom of multivariate vector

Rather than this, you can just put a real prior on sigma and nu and leave
them unconstrained (well, they need lower=0 at least).

5. Given that you only have a two-dimensional Cholesky factor, you could
just parameterize directly on the correlation rho.

- Bob

> --
> 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.

> <multiT.stan>

[Ravi Sastry]

Bob,

Thanks for your very helpful comments! It's running much faster with moderate dimensionality (J <= 20). 

My real dataset will have J>=100, so we'll see how well it scales up. 

Ravi