Integracja Metropolis-Hastings - dlaczego moja strategia nie działa?

Załóżmy, że mam funkcję , którą chcę zintegrować Oczywiście przy założeniu, że osiąga zero w punktach końcowych, brak wybuchów, fajna funkcja. Jednym ze sposobów, w jakie się bawiłem, jest użycie algorytmu Metropolis-Hastings do wygenerowania listy próbek z rozkładu proporcjonalnego do , w którym brakuje stałej normalizacyjnej który , a następnie obliczenia niektórych statystyk na tych : $g(x)$

$$\int_{-\infty}^\infty g(x) dx.$$ $g(x)$$x_1, x_2, \dots, x_n$$g(x)$ $$N = \int_{-\infty}^{\infty} g(x)dx$$ $p(x)$$f(x)$$x$ $$\frac{1}{n} \sum_{i=0}^n f(x_i) \approx \int_{-\infty}^\infty f(x)p(x)dx.$$

Ponieważ , mogę podstawić aby anulować g z całki, co powoduje wyrażenie postaci \ frac {1} {N} \ int _ {- \ infty} ^ {\ infty} \ frac {U (x)} {g (x)} g (x) dx = \ frac {1} {N} \ int _ {- \ infty} ^ \ infty U (x) dx. Więc pod warunkiem, że U (x) integruje się z 1 wzdłuż tego regionu, powinienem uzyskać wynik 1 / N , który mógłbym wziąć na zasadzie wzajemności, aby uzyskać odpowiedź, której chcę. Dlatego mógłbym wziąć zakres mojej próbki (aby najskuteczniej wykorzystać punkty) r = x_ \ max - x_ \ min i pozwolić U (x) = 1 / r dla każdej narysowanej próbki. W ten sposób U (x)$p(x) = g(x)/N$$f(x) = U(x)/g(x)$$g$

$$\frac{1}{N}\int_{-\infty}^{\infty}\frac{U(x)}{g(x)} g(x) dx = \frac{1}{N}\int_{-\infty}^\infty U(x) dx.$$ $U(x)$$1$$1/N$$r = x_\max - x_\min$$U(x) = 1/r$$U(x)$ocenia na zero poza regionem, w którym nie ma moich próbek, ale integruje się z $1$ w tym regionie. Więc jeśli teraz przyjmuję oczekiwaną wartość, powinienem otrzymać: $$E\left [\frac{U(x)}{g(x)}\right ] = \frac{1}{N} \approx \frac{1}{n} \sum_{i=0}^n \frac{U(x)}{g(x)}.$$

Próbowałem to przetestować w R dla przykładowej funkcji $g(x) = e^{-x^2}$ . W tym przypadku nie używam Metropolis-Hastings do generowania próbek, ale używam rzeczywistych prawdopodobieństw rnormdo generowania próbek (tylko do testowania). Nie do końca otrzymuję wyniki, których szukam. Zasadniczo pełne wyrażenie tego, co bym obliczał, to:

$$\frac{1}{n(x_{\max} - x_\min)} \sum_{i=0}^n \frac{1}{ e^{-x_i^2}}.$$ To powinno w mojej teorii oceniać na $1/\sqrt{\pi}$ . Zbliża się, ale na pewno nie zbiega się w oczekiwany sposób, czy robię coś złego?

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

Edycja dla CliffAB

Powodem, dla którego używam zakresu jest po prostu łatwe zdefiniowanie funkcji, która nie jest zerowa w regionie, w którym znajdują się moje punkty, ale która integruje się z w zakresie . Pełna specyfikacja funkcji jest następująca: Nie musiałem używać jako tej jednolitej gęstości. Mógłbym użyć innej gęstości zintegrowanej z , na przykład gęstości prawdopodobieństwa Sprawiłoby to jednak, że sumowanie poszczególnych próbek byłoby banalne, tzn [ - ∞ , ∞ ] U ( x ) = { $1$$[-\infty, \infty]$U(x)1P(x)=

$$U(x) = \begin{cases} \frac{1}{x_\max - x_\min} & x_\max > x > x_\min \\ 0 & \text{otherwise.} \end{cases}$$ $U(x)$$1$ $$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}}.$$

Mógłbym wypróbować tę technikę dla innych dystrybucji, które integrują się z . Chciałbym jednak nadal wiedzieć, dlaczego nie działa w przypadku równomiernej dystrybucji.$1$

Jest to najciekawsze pytanie, które dotyczy zagadnienia przybliżenia stałej normalizującej gęstości na podstawie wyjścia MCMC z tej samej gęstości . (Uwaga boczna jest taka, że ​​poprawnym założeniem jest to, że jest liczbą całkowitą, zejście do zera w nieskończoności nie jest wystarczające.)$g$$g$$g$

Moim zdaniem najbardziej odpowiedni wpis na ten temat w odniesieniu do twojej sugestii to artykuł Gelfanda i Deya (1994, JRSS B ), w którym autorzy opracowali bardzo podobne podejście do znalezienia przy generowaniu z . Jednym z rezultatów tego artykułu jest to, że dla dowolnej gęstości prawdopodobieństwa [jest to równoważne twojej ], tak że następująca tożsamość pokazuje, że próbka z może wytworzyć

$$\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\}$$ $$\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}$$ $p$nieobciążony oceny od estymatora próbkowania znaczenie Oczywiście wydajność (szybkość zbieżności, istnienie wariancji i tc.) estymatora zależy od wyboru [ nawet jeśli jego oczekiwania nie]. W ramach bayesowskich opcją zalecaną przez Gelfanda i Deya jest przyjęcie , poprzedniej gęstości. Prowadzi to do gdzie jest funkcją prawdopodobieństwa, ponieważ$1/N$ $$\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)$$ $\hat\eta$$\alpha$$\alpha=\pi$ $$\dfrac{\alpha(x)}{g(x)} = \dfrac{1}{\ell(x)}$$ $\ell(x)$$g(x)=\pi(x)\ell(x)$. Niestety, uzyskany estymator jest średnią harmoniczną estymator , zwany również najgorszym Monte Carlo estymator nigdy przez Radford Neal, z University of Toronto. Więc nie zawsze działa to dobrze. Lub nawet prawie nigdy. $$\hat{N}=\dfrac{n}{\sum_{i=1}^n1\big/\ell(x_i)}$$

Twój pomysł użycia zakresu próbki i jednolitości w tym zakresie jest związany ze średnią średnią harmoniczną: estymator nie ma wariancji choćby dlatego, że z powodu pojawiający się w liczniku (podejrzewam, że zawsze może tak być w przypadku nieograniczonego wsparcia!), a zatem bardzo powoli zbliża się do stałej normalizującej. Na przykład, jeśli ponownie uruchomisz kod kilka razy, otrzymasz bardzo różne wartości liczbowe po 10 different iteracjach. Oznacza to, że nie możesz nawet ufać wielkości odpowiedzi.$(\min(x_i),\max(x_i))$$\exp\{x^2\}$

Ogólną poprawką tego nieskończonego problemu wariancji jest użycie dla bardziej skoncentrowanej gęstości, przy użyciu na przykład kwartyli z próbki , ponieważ następnie pozostaje niższy w tym przedziale.$\alpha$$(q_{.25}(x_i),q_{.75}(x_i))$$g$

Podczas dostosowywania kodu do tej nowej gęstości przybliżenie jest znacznie bliższe :$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

Omawiamy tę metodę szczegółowo w dwóch artykułach z Darrenem Wraithem i Jean-Michelem Marinem .