0
$\begingroup$

I'm attempting to implement a Metropolis-Hastings Algorithm to evaluate integrals of the following form:

$$I =\frac{1}{\sqrt\pi}\int_{-\infty}^{\infty} {f(x)\exp(-x^2)} \text{d}x$$

Now we can rewrite this as:

$$I =\int_{-\infty}^{\infty} {f(x)P(x)} \text{d}x$$

Where:

$$ P(x) = \mathcal{N}(\mu=0,\,\sigma^{2}=\frac{1}{\sqrt2}) $$

Now I believe we should generate $X^*$ from some density $q(X^*|X_n)$

Generally we use $q(X^*|X_n) \sim \mathcal{N}(X_n,\,\sigma^{2}) $

We accept $X^* \ as\ X_{n+1}$ with probability:

$ \alpha(X^*, X_n)= min(1, \frac{P(X^*)}{P(X)})$

Else $X_{n+1} = X_n$

So that's how we sample values.

Will the integral then tend towards:

$I \approx \frac{1}{N} \sum f(X_i)$

$\endgroup$

marked as duplicate by Xi'an, Michael Chernick, mdewey, John, kjetil b halvorsen Dec 31 '17 at 21:42

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • $\begingroup$ I'm not an expert, but I don't understand why one would want anything besides a uniform distribution . The M-H algo presumes that there is a prior distribution on the probabilities of going from one state to another. But 'dx' means, if you will, sample at random. $\endgroup$ – aginensky Dec 30 '17 at 16:20
  • $\begingroup$ I think we sample from the normal distribution because it is symmetric. If you look at the tag description for metropolis-hastings it tells us: $$\alpha=\min\left(1,\frac{\mathscr{p}(x^\star)\mathscr{q}(x_t|x^\star)}{\mathscr{p}(x_t)\mathscr{q}(x^*|x_t)}\right)$$ This however simplifies if q is symmetric $\endgroup$ – John Meighan Dec 30 '17 at 16:31
  • 2
    $\begingroup$ @John Meighan is this something as a toy example or to study M-H algorithm? Generally you use MH and more generally Monte Carlo integration when you cannot write down the integral in closed form or you are in a high dimensional space. A standard 1-d quadrature rule should work perfectly fine to evaluate this integral for any specified $f(x)$. $\endgroup$ – Lucas Roberts Dec 30 '17 at 16:34
  • $\begingroup$ @LucasRoberts Yea I just want to study M-H algorithm. $\endgroup$ – John Meighan Dec 30 '17 at 17:28
  • $\begingroup$ In going from the first expression for $I$ to the second expression, I believe the factor $1/\sqrt{\pi}$ gets absorbed into $P(x)$. $\endgroup$ – mef Dec 30 '17 at 20:03
1
$\begingroup$

I believe your conjecture is correct.

However, in this case there is a better proposal distribution. As you note the the comments, the acceptance probability is given by \begin{equation} \alpha = \min\left(1, \frac{p(x^*)\,q(x_t|x^*)}{p(x_t)\,q(x^*|x_t)}\right) . \end{equation} Let $q(x^*|x_t) =\mathcal{N}(x^*|0,1/2)$. Then $p(x^*) = q(x^*|x_t)$ and $p(x_t) = q(x_t|x^*)$ and consequently $\alpha \equiv 1$. In other words you always accept with this proposal.

$\endgroup$
  • 1
    $\begingroup$ This is, of course, more efficient...but also kinda cheating because it takes advantage of the fact that we can already draw directly from the target distribution. If that's the case, MH sampling is in no way necessary. So in real applications, we will never have that $p = q$ (or we shouldn't be using MH) $\endgroup$ – Cliff AB Dec 30 '17 at 20:27
  • $\begingroup$ In real world applications you may have a Gibbs sampler, which can be expressed as a Metropolis-Hastings sampler where, because of the cancellation involved, produces $\alpha \equiv 1$. $\endgroup$ – mef Dec 30 '17 at 23:33

Not the answer you're looking for? Browse other questions tagged or ask your own question.