Using Stan with PDE-based inverse problems

General
1

I work in geophysics inverse problems and have written my own NUTS code to find the posterior distribution of model parameters wherein the model is a PDE. I’m wondering if it would be possible to use Stan to do this instead?

I think I need:

  1. a way to specify a custom likelihood function. The likelihood involves the norm of the misfit between observed data and simulated (with the PDE) data.
  2. a way to specify a custom gradient for the likelihood. The gradient is not well suited to autodiff.
2

If you had the output of the PDE, z, and the data, y, could you write the likelihood in Stan code?

The easiest way to do this is probably interface C++ with rstan: Interfacing with External C++ Code • rstan . If your PDE is super simple you might go that route.

–

The robust way to do this at this point is fork cmdstan/stan/stanc3/math and roll your own cmdstan. The cmdstanr interface can wrap around cmdstans from custom repos/forks: Install CmdStan or clean and rebuild an existing installation — install_cmdstan • cmdstanr .

This is pretty involved process though. If the code was super simple C++, maybe it’d take about a month to go through this process if everything went your way? You’d want to write your PDE interface and such in the math library. Do tests there of the autodiff and stuff, and modify those makefiles to include your PDE code.

Then you’d need to modify the stanc3 compiler to expose your functions to Stan to the language. I think there’s a tutorial around on how to do this somewhere.

That’d be the costs at least?

1 Like
3

This is a thread where @maxbiostat and some folks went through the process of doing some custom C++: Problem implementing Zipf (power law) distribution -- external C++ . Might be helpful to look through there. There are a lot of compiler error things going on, but you’ll probably get to work through a lot of weird compiler errors however you do this lol.

1 Like
4

Thanks a bunch for the info! It looks like this is going to be pretty involved.

5

Yes. See https://mc-stan.org/docs/2_18/stan-users-guide/user-defined-probability-functions.html. What you need is to feed your PDE solver’s cost function to your user-define lpdf function, this needs to happen in c++ level.

Not really. You’ll have to come up a C++ function that takes value & gradient and assemble Stan’s var variables. I explored the idea in https://metrumrg.com/wp-content/uploads/2018/09/stan_pde_poster.pdf and actually working on something similar for a specific FEM solver.

1 Like
6

How big problems are you going to solve?