Covariance Matrix Sensitivity

[Michael Betancourt]

Where is the covariance matrix type defined? I'm having a hell of a time where it shows up in the model tests as it seems to spit out zero-valued diagonals even for seemingly reasonable (O(1)) values of the unconstrained parameters, momenta, and gradients. For example, models/bugs_examples/vol3/camel/camel2 dies after a single iteration with either the transformed parameter rho becomes a nan (because Sigma[1,1] or Sigma[2,2] is vanishing) or Sigma itself throws a "not positive-definite" exception.

Additionally, if anyone can tell me how to access constrained parameters via the model object I'd be very much appreciative.

Thanks!

-Mike

[Ben Goodrich]

On Sunday, April 7, 2013 5:06:17 PM UTC-4, Michael Betancourt wrote:

Where is the covariance matrix type defined?  I'm having a hell of a time where it shows up in the model tests as it seems to spit out zero-valued diagonals even for seemingly reasonable (O(1)) values of the unconstrained parameters, momenta, and gradients.  For example, models/bugs_examples/vol3/camel/camel2 dies after a single iteration with either the transformed parameter rho becomes a nan (because Sigma[1,1] or Sigma[2,2] is vanishing) or Sigma itself throws a "not positive-definite" exception.

Conceptually, it is defined as LL' where each diagonal element of L is exp(theta). In the camel2 example, can you get it to work if you change the transformed parameters block into a generated quantities block (since rho isn't used in the model)?

Ben

[Michael Betancourt]

For a 2D covariance as in the camel model, to what does the third parameter refer? I'm mostly interested in understanding the numerical sensitivity.

I could push rho to the generated quantities block, but the covariance matrix would still end up being ill-posed. Something odd is going on in the transformation and I'm not sure why -- initial points are well-behaved, the initial log probability is reasonable, and the step-size is appropriate for the gradients.

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Ben Goodrich]

On Sunday, April 7, 2013 7:24:46 PM UTC-4, Michael Betancourt wrote:

For a 2D covariance as in the camel model, to what does the third parameter refer?  I'm mostly interested in understanding the numerical sensitivity.

I think log(L[2,2]).
 

I could push rho to the generated quantities block, but the covariance matrix would still end up being ill-posed.  Something odd is going on in the transformation and I'm not sure why -- initial points are well-behaved, the initial log probability is reasonable, and the step-size is appropriate for the gradients.

Yeah, I can get exceptions when I do that model in your hmc_refactor branch. I don't recall there being a problem with this model in master.

Ben
 

[Michael Betancourt]

Yeah, the covariance seems to be --really-- sensitive to the unconstrained parameters. For the relatively boring parameters

-2.10605 3.1919 -2.23932

the gradients are HUGE,

-243723 75937.1 -242729

Since that position isn't far from the random initialize range of [-2, 2] that might be the source of the problem. But I'll keep looking into it, especially to try to identify the difference between master and hmc_refactor.

[Bob Carpenter]

Did you change the MCMC algorithms? I don't understand why
your samplers would be behaving differently than the ones
currently in master.

The transform is described near the end of the manual in
the transforms chapter.

For the code, you can look at stan/prob/transform to see
where the transform is happening (through io/reader and io/writer).

In the generated code for a model, look at
transform_inits() for where a writer is used
for the transforms for variables read from
a context (base class implemented by file for the
command line and through R for RStan).

You then see the constraining transform getting called
through a reader in the log_prob functions.

The write_array() method takes unconstrained params and
inverse transforms to constrained parameters and then
writes the real ones to an array. This follows the
same format as write_csv(), writing the parameters out
in row-major order in the order they were declared.
The output is constrained. write_csv_header() writes
out the "names" of the parameters, including the
indices of an array --- this is the header dumped out
in the output from the command.

- Bob

[Michael Betancourt]

> Did you change the MCMC algorithms? I don't understand why
> your samplers would be behaving differently than the ones
> currently in master.

The only real change is going from different "effective" step sizes for each parameter to a single step size and a proper mass matrix. This required significant changes in the adaptation, but the new code works for most models.

The problem with the covariance matrix models isn't with the MCMC - the evolution diverges but the resulting proposal would just be rejected. The diverging state, however, causes the covariance transform to throw an exception before the transition can be rejected. I'm not sure why the previous samplers aren't vulnerable to the same problem.

> The transform is described near the end of the manual in
> the transforms chapter.
>
> For the code, you can look at stan/prob/transform to see
> where the transform is happening (through io/reader and io/writer).
>
> In the generated code for a model, look at
> transform_inits() for where a writer is used
> for the transforms for variables read from
> a context (base class implemented by file for the
> command line and through R for RStan).
>
> You then see the constraining transform getting called
> through a reader in the log_prob functions.
>
> The write_array() method takes unconstrained params and
> inverse transforms to constrained parameters and then
> writes the real ones to an array. This follows the
> same format as write_csv(), writing the parameters out
> in row-major order in the order they were declared.
> The output is constrained. write_csv_header() writes
> out the "names" of the parameters, including the
> indices of an array --- this is the header dumped out
> in the output from the command.

Ah, okay.

[Bob Carpenter]

The previous samplers caught exceptions and treated
them as an instant rejection. They also did the same
thing with -infinity and NaN returns.

We have our functions configured with exceptions
so that then we can record an error message to report
when it's caught.

So far, the messages aren't very good, because they're
not located with respect to the user model -- they just
say which function was called and what the problem was
with the parameters it was called with.

To get better messages, we'd need to pass more information
to the functions from the model, such as the line number
of the call in the Stan model (not the C++ code).
Or we could dream up some alternative. (For instance, we could
put a try/catch around every function call :-))

- Bob

[Michael Betancourt]

Ah, found it! The code in master explicitly catches the grad_log_prob exception, a step which I had initially put aside (and then forgot about).