Intégration Metropolis-Hastings - pourquoi ma stratégie ne fonctionne-t-elle pas?

Supposons que j'ai une fonction $g(x)$ que je souhaite intégrer

$$\int_{-\infty}^\infty g(x) dx.$$ Bien sûr, en supposant que $g(x)$ passe à zéro aux points d'extrémité, pas d'explosions, belle fonction. Une façon avec laquelle j'ai joué est d'utiliser l'algorithme Metropolis-Hastings pour générer une liste d'échantillons $x_1, x_2, \dots, x_n$ partir de la distribution proportionnelle à $g(x)$, Qui manque la constante de normalisation $$N = \int_{-\infty}^{\infty} g(x)dx$$ que nous appellerons $p(x)$ , puis en calculant une certaine statistique $f(x)$ sur ces $x$ « s: $$\frac{1}{n} \sum_{i=0}^n f(x_i) \approx \int_{-\infty}^\infty f(x)p(x)dx.$$

Puisque $p(x) = g(x)/N$ , je peux substituer dans $f(x) = U(x)/g(x)$ pour annuler $g$ de l'intégrale, résultant en une expression de la forme

$$\frac{1}{N}\int_{-\infty}^{\infty}\frac{U(x)}{g(x)} g(x) dx = \frac{1}{N}\int_{-\infty}^\infty U(x) dx.$$ Donc, à condition que$U(x)$s'intègre à$1$long de cette région, je devrais obtenir le résultat$1/N$, que je pourrais simplement prendre l'inverse pour obtenir la réponse que je veux. Par conséquent, je pouvais prendre la plage de mon échantillon (pour utiliser le plus efficacement possible les points)$r = x_\max - x_\min$et laisser$U(x) = 1/r$pour chaque échantillon que j'ai tiré. De cette façon s'évalue à zéro en dehors de la région où mes échantillons ne sont pas, mais s'intègre à 1 dans cette région. Donc, si je prends maintenant la valeur attendue, je devrais obtenir: E [ U ( x $U(x)$$1$ $$E\left [\frac{U(x)}{g(x)}\right ] = \frac{1}{N} \approx \frac{1}{n} \sum_{i=0}^n \frac{U(x)}{g(x)}.$$

J'ai essayé de tester cela dans R pour l'exemple de fonction . Dans ce cas, je n'utilise pas Metropolis-Hastings pour générer les échantillons mais j'utilise les probabilités réelles avec pour générer des échantillons (juste pour tester). Je n'obtiens pas tout à fait les résultats que je recherche. Fondamentalement, l'expression complète de ce que je calculerais est: $g(x) = e^{-x^2}$rnorm Cela devrait dans ma théorie évaluer à1/

$$\frac{1}{n(x_{\max} - x_\min)} \sum_{i=0}^n \frac{1}{ e^{-x_i^2}}.$$ . Il se rapproche mais il ne converge certainement pas de la manière attendue, est-ce que je fais quelque chose de mal?$1/\sqrt{\pi}$

ys = rnorm(1000000, 0, 1/sqrt(2))
r = max(ys) - min(ys)
sum(sapply(ys, function(x) 1/( r * exp(-x^2))))/length(ys)
## evaluates to 0.6019741. 1/sqrt(pi) = 0.5641896

Modifier pour CliffAB

La raison pour laquelle j'utilise la plage est juste pour définir facilement une fonction qui est non nulle sur la région où se trouvent mes points, mais qui s'intègre à sur la plage [ - ∞ , ∞ ] . La spécification complète de la fonction est: U ( x ) = { $1$$[-\infty, \infty]$ Je n'ai pas eu à utiliserU(x)comme cette densité uniforme. J'aurais pu utiliser une autre densité intégrée à1, par exemple la densité de probabilité P(x)=

$$U(x) = \begin{cases} \frac{1}{x_\max - x_\min} & x_\max > x > x_\min \\ 0 & \text{otherwise.} \end{cases}$$ $U(x)$$1$ Cependant, cela aurait rendu la somme des échantillons individuels triviale, c'est-à-dire $$P(x) = \frac{1}{\sqrt{\pi}} e^{-x^2}.$$ $$\frac{1}{n} \sum_{i=0}^n \frac{P(x)}{g(x)} = \frac{1}{n} \sum_{i=0}^n \frac{e^{-x_i^2}/\sqrt{\pi}}{e^{-x_i^2} } = \frac{1}{n} \sum_{i=0}^n \frac{1}{\sqrt{\pi}} = \frac{1}{\sqrt{\pi}}.$$

Je pourrais essayer cette technique pour d'autres distributions qui s'intègrent à . Cependant, je voudrais quand même savoir pourquoi cela ne fonctionne pas pour une distribution uniforme.$1$

$g$$g$$g$

$$\int_\mathcal{X} g(x) \,\text{d}x$$ $p(x)\propto g(x)$$\alpha(x)$$U(x)$ $$\{x;\alpha(x)>0\}\subset\{x;g(x)>0\}$$ the following identity $$\int_\mathcal{X} \dfrac{\alpha(x)}{g(x)}p(x) \,\text{d}x=\int_\mathcal{X} \dfrac{\alpha(x)}{N} \,\text{d}x=\dfrac{1}{N}$$ shows that a sample from $p$ can produce an unbiased evaluation of $1/N$ by the importance sampling estimator $$\hat\eta=\frac{1}{n}\sum_{i=1}^n \dfrac{\alpha(x_i)}{g(x_i)}\qquad x_i\stackrel{\text{iid}}{\sim}p(x)$$ Obviously, the performances (convergence speed, existence of a variance, &tc.) of the estimator $\hat\eta$ do depend on the choice of $\alpha$ [even though its expectation does not]. In a Bayesian framework, a choice advocated by Gelfand and Dey is to take $\alpha=\pi$, the prior density. This leads to $$\dfrac{\alpha(x)}{g(x)} = \dfrac{1}{\ell(x)}$$ where $\ell(x)$ is the likelihood function, since $g(x)=\pi(x)\ell(x)$. Unfortunately, the resulting estimator $$\hat{N}=\dfrac{n}{\sum_{i=1}^n1\big/\ell(x_i)}$$ is the harmonic mean estimator, also called the worst Monte Carlo estimator ever by Radford Neal, from the University of Toronto. So it does not always work out nicely. Or even hardly ever.

Your idea of using the range of your sample $(\min(x_i),\max(x_i))$ and the uniform over that range is connected with the harmonic mean issue: this estimator does not have a variance if only because because of the $\exp\{x^2\}$ appearing in the numerator (I suspect it could always be the case for an unbounded support!) and it thus converges very slowly to the normalising constant. For instance, if you rerun your code several times, you get very different numerical values after 10⁶ iterations. This means you cannot even trust the magnitude of the answer.

A generic fix to this infinite variance issue is to use for $\alpha$ a more concentrated density, using for instance the quartiles of your sample $(q_{.25}(x_i),q_{.75}(x_i))$, because $g$ then remains lower-bounded over this interval.

When adapting your code to this new density, the approximation is much closer to $1/\sqrt{\pi}$:

ys = rnorm(1e6, 0, 1/sqrt(2))
r = quantile(ys,.75) - quantile(ys,.25)
yc=ys[(ys>quantile(ys,.25))&(ys<quantile(ys,.75))]
sum(sapply(yc, function(x) 1/( r * exp(-x^2))))/length(ys)
## evaluates to 0.5649015. 1/sqrt(pi) = 0.5641896

We discuss this method in details in two papers with Darren Wraith and with Jean-Michel Marin.