# expensive, numerical posterior

64 views

### Tamas Papp

Sep 4, 2015, 3:46:03 AM9/4/15
to Stan users mailing list
Hi,

specific model/Stan feature, but I hope I could get some help or
suggestions here; as I want to use Stan for something I have not seen it
used for so I don't know whether my plans are feasible or even make
sense.

I have an economic model for which the likelihood/posterior p(theta|y)
is

1. expensive to calculate (takes about 15s on a reasonably good
computer for a single point),

2. is only available through a numerical algorithm (the economic model
needs to be solved for the law of motion), so there is no derivative
information and I cannot code it in Stan.

Traditionally in economics RMW (Random Walk Metropolis) is suggested for
these problems, but I was wondering if I could somehow leverage the
algorithms of Stan for it, for more efficient sampling.

I was thinking about the following:

a. obtain an approximation q(theta|y) to p in a representation that Stan
can calculate efficiently,

b. use rejection sampling for the result to obtain samples from p.

Reading the approximation: looking at the literature the algorithm
suggested is a normal approximation around some mode. But the problem
seems to be severely multimodal and I doubt that this would be efficient.

I think I can obtain a state-space form, extending with unobserved state
x and make it a Kalman filter where matrices depend on theta, for which
I would use a linear, Pade or a global approximation then marginalize.

But I am not clear on where the approximation should be "good": ideally
it would be best on the region where most of the posterior is, but then
I only know that once I do the rejection sampling. Would it then make
sense to iterate (find another approximation knowing the posterior,
etc)?

Suggestions, pointers to the literature, etc, would be very
welcome. Sorry if the question is vague.

Best,

Tamas

### Michael Betancourt

Sep 4, 2015, 5:11:28 AM9/4/15
No, no, no, no!

The problem with these kinds of approximations (often called emulators
in the literature) is that interpolations in high dimensions are ill-posed —
there are a huge number of possible interpolations and only a singular
few are actually consistent with the truth. So what happens in practice
is you try to interpolate your few numerical evaluations but most of the
resulting proposals fail.

What you really want to do is open up the black box of your economic
model. Go back to the diff eqs, for example, and figure out the gradients
from there and run HMC. You can already see that we’re starting to
do this with explicit ODE models, but you’ll just have to figure out what
those ODEs are. Yes, those gradients can be expensive but the sampling
speed up you get over something like RWM is huge and makes it all
the worth while.
> --
> 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.

### Tamas Papp

Sep 4, 2015, 5:31:45 AM9/4/15
The problem is that most economic models are forward-looking and solving
for the coefficients of the ODEs is the nontrivial step, imposed by
equilibrium conditions, usually done by iterative numerical methods. I
don't think one can code that in Stan, but I will look into it.

On the other hand, I think can provide a linear approximation to the
model, and then an approximation of the parameters of that as a function
of the model parameters, which I believe is good or at least OK.

Thanks,

Tamas

### Michael Betancourt

Sep 4, 2015, 6:15:00 AM9/4/15

> The problem is that most economic models are forward-looking and solving
> for the coefficients of the ODEs is the nontrivial step, imposed by
> equilibrium conditions, usually done by iterative numerical methods. I
> don't think one can code that in Stan, but I will look into it.

There are two ways of using these models. One is to code up the solver
in Stan and let Stan automatically differentiate through the solver. This
works but tends to be slow. Another is to differentiate the original system
and then figure out how to solve for those derivatives. We do this
automatically for ODEs in Stan, but can be done analytically for all kinds
of systems such as implicit constraints that require root solvers.

> On the other hand, I think can provide a linear approximation to the
> model, and then an approximation of the parameters of that as a function
> of the model parameters, which I believe is good or at least OK.

Then go ahead and try it. But keep in mind that for anything but really simple
models (which wouldn’t then require such complex solvers) these approximations
fail hard pretty quickly.

These are not easy problems — but if your likelihood is a black box
then there’s only much you can do to explore it efficiently. You _have_
to open that box, which is usually a good idea for reproducible science
anyways.

### Bob Carpenter

Sep 4, 2015, 12:50:55 PM9/4/15

> On Sep 4, 2015, at 6:15 AM, Michael Betancourt <betan...@gmail.com> wrote:
>
>
>> The problem is that most economic models are forward-looking and solving
>> for the coefficients of the ODEs is the nontrivial step, imposed by
>> equilibrium conditions, usually done by iterative numerical methods. I
>> don't think one can code that in Stan, but I will look into it.
>
> There are two ways of using these models. One is to code up the solver
> in Stan and let Stan automatically differentiate through the solver. This
> works but tends to be slow. Another is to differentiate the original system
> and then figure out how to solve for those derivatives. We do this
> automatically for ODEs in Stan, but can be done analytically for all kinds
> of systems such as implicit constraints that require root solvers.

Yup --- this is on our to-do list and will be part of
our next grant proposal.

>> On the other hand, I think can provide a linear approximation to the
>> model, and then an approximation of the parameters of that as a function
>> of the model parameters, which I believe is good or at least OK.
>
> Then go ahead and try it. But keep in mind that for anything but really simple
> models (which wouldn’t then require such complex solvers) these approximations
> fail hard pretty quickly.
>
> These are not easy problems — but if your likelihood is a black box
> then there’s only much you can do to explore it efficiently. You _have_
> to open that box, which is usually a good idea for reproducible science
> anyways.

I think the source is often available for these "black boxes". Just
not an algorithm that's easy to differentiate through. Or to translate to
Stan if it's person-years of specialized code. At least that
was the case with all the energy models we were talking about with Phil
Price.

- Bob