Using MCMC to evaluate the expected value of a high-dimensional function

Question

I am working on a research project that is related to optimization and recently had an idea to use MCMC in this setting. Unfortunately, I am fairly new to MCMC methods so I had several questions. I'll start by describing the problem and then asking my questions.

Our problem boils down to estimating the expected value of a cost function $c(\omega)$ where $\omega = (\omega_1,\omega_2,...\omega_h)$ is an $h$-dimentional random variable with a density $f(\omega)$.

In our case, a closed form version of $c(\omega)$ does not exist. This means that we have to use Monte Carlo methods to approximate the expected value. Unfortunately, it turns out that estimates of $E[c(\omega)]$ that are generated using MC or QMC methods have too much variance to be useful within in a practical setting.

One idea that we had to use an importance sampling distribution to generate sample points that will produce a low variance estimate of $E[c(\omega)]$. In our case, the ideal importance sampling distribution, $g(\omega)$, has to be roughly proportional to $c(\omega)f(\omega)$. Seeing how $g(\omega)$ is known up to constant, I am wondering whether I can use MCMC along with the proposal distribution $c(\omega)f(\omega)$ to eventually generate samples from $g(\omega)$ .

My questions here are:

• Can MCMC be used within this setting? If so, what MCMC method would be appropriate? I am working in MATLAB, so I have a preference to anything that already has a MATLAB implementation.

• Are there any techniques that I can use to speed up the burn-in period for MCMC. And how can I tell that the stationary distribution has been reached? In this case, it actually takes a fair bit of time to calculate $c(\omega)$ for a given $\omega$.

Answer 1

I would always remember, that MCMC is just a numerical integration tool (and a rather inefficient one at that). It is not some magic/mystical thing. It is very useful because it is reasonably easy to apply. It does not require much thinking compared to some other numerical integration techniques. For instance, you do not have to do any derivatives. You only have to generate "random numbers".

However, like any numerical integration method, it is not a universal catch all tool. There are conditions when it is useful, and conditions when it isn't.

It may be wiser to set up another technique. Depending on how big $h$ is, and how fast your computer is, and how much time you are prepared to wait for results. A uniform grid may do the job (although this requires small $h$ or a long amount of waiting). The "job" is to evaluate the integral - the equation does not care what meaning you or I attach to the result (and hence it does not care whether we obtained the result randomly or not).

Additionally, if your estimates of $\omega$ are quite accurate, the $f(\omega)$ will be sharply peaked and closely resemble a delta function, so the integral is effectively substituting $\omega\rightarrow\omega_{max}$.

Another numerical integration technique is using a taylor series under the integral. $f(\omega)\approx f(\omega_{max})+(\omega-\omega_{max})f'(\omega_{max})+\frac{1}{2}(\omega-\omega_{max})^{2}f''(\omega_{max})+\dots$

This is a useful strategy when the moments of $\omega$ are easily obtained.

Edwin Jaynes has a nice quote on this:

whenever there is a randomised way of doing something, there is a non-randomised way which yields better results, but requires more thinking

One "more thinking" way is to use "stratified MCMC" to do the integral. So rather than "randomly" pick a spot on the whole parameter space: divide it up into "strata". These "strata" should be picked so that you get a good range of the high part of the integral. Then randomly sample within each strata. But this will require you to write your own code I would imagine (i.e. more thinking).

Answer 2

There isn't any indication that your variables here are correlated so I dont know why you would use MCMC as opposed to regular Monte Carlo. There are many different sampling methods including the mentioned stratified sampling (Latin hypercube) and QMC. Sparse quadrature methods are very good if the dimension of the problem is not too high (not more than 10) since sparse quadrature grids grow geometrically (curse of dimensionality).

But it sounds like you are on the right track with respect to importance sampling. The key here is to choose a biased distribution that has large probability concentrated near your region of interest and that it has thicker tails than the nominal distribution.

I'd like to add that this is an open research problem so if you can come up with something good it would be of great interest to the community!