Metropolis Algorithm for approximating integrals

Question

I've been trying to approximate $\int_a^b f(x)\hspace{1pt} dx$ with the weight function $W(x)=A(e^{x^2}-1)$, $A$ is a constant of normalization (when $x\in[0,1]$ then $A=2.16145...$). I don't have any problems when the integration limits are $x \in [0,1]$, but the general case when $x \in [a,b]$ is cumbersome.

This is the procedure I followed for the $x \in [0,1]$ case:

1. Choose a number $x_0$ in $[0,1]$
2. Now set $x_t$ (trial step) to be $x_t=x_0 +$ Random[{-h,h}]
3. Then set r to be: $r=\frac{W(x_t)}{W(x_0)}$ and compare it to some random number $\Omega_i \in \{0,1\}$
4. If $r\geq\Omega_i$ then $x_t$ is accepted (i.e. $x_0=x_t$), but if $r<\Omega_i$ then $x_t$ is set to be $x_t= x_0$.
5. Repeat discarding the $m$ first steps. (That is, to remove any dependence on the starting point, the algorithm is run for some large number of steps which are then discarded).
6. Now calculate the integral $\int_a^b f(x)\hspace{1pt} dx$ as:

$\int_0^1 f(x)\hspace{1pt} dx=\int_0^1 \frac{f(x)}{W(x)}\hspace{1pt}W(x)dx\approx \frac{1}{(N-m)} \sum\limits_{i=1}^{N-m} \frac{f(x_i)}{W(x_i)}$.

I run the program for little $N$ to fixate h such as the acceptance ratio of $x_t$ is $\approx 50$% and then I re-run for a large $N$ with this $h$.

With $W(x)=A(e^{x^2}-1)$, $f(x)=x^2$, $N=1\times 10^6$, $m=1\times 10^5$ and $h=0.8$ (which gives an acceptance ratio of $50.8$%), I get $\int_0^1 x^2 dx=0.333305$, which is close enough to the actual value $\frac{1}{3}$;

What do I have to do in order for accounting the general integral limits? Do you have any suggestion? Do you suggest any book or paper I can read? Thank you in advance for your time.

P.S. English is not my first language.

Answer 1

The same exact procedure (which qualifies as importance sampling implemented via Metropolis-Hastings) applies to any integration range $(a,b)$.

• Select an importance probability density $w(\cdot)$ over $(a,b)$ [for instance, $w(x)\propto e^{x^2}-1$]
• Simulate a Markov chain $(X_t)_t$ with stationary distribution $w(\cdot)$
• Approximate the integral $\int_a^b f(x)\text{d}x$ by the empirical average$$\frac{1}{T}\sum_{t=1}^T f(x_t)/w(x_t)$$or, if the normalising constant is unknown,$$\sum_{t=1}^T f(x_t)/w(x_t)\Bigg/\sum_{t=1}^T 1/w(x_t)(b-a)$$

The bottom approximation is based on the identity \begin{align*}1&=\int_a^b \frac{1}{b-a}\text{d}x\\&=\int_a^b \frac{1}{b-a}\frac{w(x)}{w(x)}\text{d}x\\&=\int_a^b \frac{A}{(b-a)w(x)}\frac{w(x)}{A}\text{d}x\\&=\mathbb{E}\left[\frac{A}{(b-a)w(X)}\right]\end{align*} (It is a safe version of the harmonic mean estimator.)

As for suggesting textbooks on MCMC integrations, there are many to pick from

• Computational Bayesian Statistics by Antónia Amaral Turkman, Carlos Daniel Paulino, and Peter Müller
• Statistical Rethinking: A Bayesian Course with Examples in R and Stan by Richard McElreath
• Bayesian Data Analysis by Andrew Gelman, John Carlin, Hal Stern, David Dunson, Aki Vehtari, and Don Rubin

(avoiding citing my own books!)

Answer 2

As Xi'an pointed out, MH algorithm is appropriate for this task. One neat part about MH is that the normalizing constant isn't required; algorithm starts with an initial position, which we'll call current. A disturbance is observed, ie you increase by +/- x (you can tune this amount, I've used 0.25 below.) And so a new variable called movement is current + disturbance. Now, we evaluate if the move is reasonable. In MH, we look at the likelihood ratio of observations, using the PDF. For your purposes, this is simply the y-value associated with your function. After comparing these values we'll get a number in the range [0,inf). Now, we observe a random event bounded by [0,1]. If the event is <= to the ratio, we accept the movement, otherwise, we remain at the current position. After some predefined number of iterations, we return the samples.

Here is some python code to illustrate.

def func(x,a,b):
    if (a<x<b):
        return np.exp((x**2)-1)
    else:
        return 0
def areaMCMC(a,b,hops=10_000):
    current = np.random.uniform(a,b)
    samples = []
    for i in range(hops):
        samples.append(current)
        disturbance = np.random.uniform(-0.25,0.25)
        movement = current + disturbance
    current_lik = func(x=current,a=a,b=b)
    movement_lik = func(x=movement,a=a,b=b)
    ratio = movement_lik/current_lik

    event = np.random.uniform(0,1)
    if event &lt;= ratio:
        current = movement
return samples


samples = areaMCMC(a=0,b=3)
import matplotlib,pyplot as plt
plt.hist(samples,density=True)

sum(samples)/len(samples)
>>>
2.7379448212941697

The sum we've just observed is the approximate area of your function in the range (0,3). Note, if you want this range to be inclusive, not exclusive, adjust func() to use <= rather than <. However, at a high enough number of iterations, the effects of the difference should be moot.

Lastly, and as touched on earlier, adjust 0.25 disturbance as necessary. When sampling (0,3) 0.25 is probably appropriate. When sampling (-1000,1000), you probably want bigger disturbances to adequately explore the function.

Edit1: For your case, it might not be true that you can disregard the normalizing constant. When approximating probability distributions, we already know that the distribution integrates to 1 and since we're only using the PDF for likelihood of a current position relative to another, we don't need the normalizing constant; it doesn't affect shape, only scale. But in your case, the area does not sum to 1. Nonetheless, if you multiply the output of this code by the normalizing constant, you'll scale the estimated area correctly.

Edit2: To account for integration limits, as you've asked, you need to return a y-value of 0 for any x-value supplied outside of this range. In other words, if the probability distribution is not defined here, the likelihood returned by the PDF will be 0.

Hope this helps!