Metropolis-Hastings-Integration - warum funktioniert meine Strategie nicht?

Angenommen, ich habe eine Funktion , die ich integrieren möchte: Natürlich unter der Annahme, dass an den Endpunkten auf Null geht, keine Blowups, nette Funktion. Eine Möglichkeit, mit der ich herumgespielt habe, besteht darin, mit dem Metropolis-Hastings-Algorithmus eine Liste der Stichproben aus der zu proportionalen Verteilung zu , in der die Normalisierungskonstante fehlt. das ich , und dann eine Statistik für diese berechnen : $g(x)$

\int_{- \infty}^{\infty} g (x) d x .

$\int_{-\infty}^\infty g(x) dx.$

g (x)

$g(x)$

x_{1}, x_{2}, \dots, x_{n}

$x_1, x_2, \dots, x_n$

g (x)

$g(x)$

N = \int_{- \infty}^{\infty} g (x) d x

$N = \int_{-\infty}^{\infty} g(x)dx$

p (x)

$p(x)$

f (x)

$f(x)$

x

$x$

\frac{1}{n} \sum_{i = 0}^{n} f (x_{i}) \approx \int_{- \infty}^{\infty} f (x) p (x) d x .

$\frac{1}{n} \sum_{i=0}^n f(x_i) \approx \int_{-\infty}^\infty f(x)p(x)dx.$

Da , kann ich , um vom Integral zu streichen , was zu einem Ausdruck der Form Vorausgesetzt also, dass entlang dieser Region zu integriert wird , sollte ich das Ergebnis , das ich einfach als Kehrwert verwenden könnte, um die gewünschte Antwort zu erhalten. Daher könnte ich den Bereich meiner Stichprobe (um die Punkte am effektivsten zu nutzen) und für jede Stichprobe, die ich gezogen habe, lassen. Auf diese Weise $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) d x = \frac{1}{N} \int_{- \infty}^{\infty} U (x) d x .

$\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)

$U(x)$

1

$1$

1 / N

$1/N$

r = x_{max} - x_{min}

$r = x_\max - x_\min$

U (x) = 1 / r

$U(x) = 1/r$

U (x)

$U(x)$ Wertet außerhalb der Region, in der meine Samples nicht vorhanden sind, auf Null aus, integriert aber in dieser Region auf

1

$1$ . Wenn ich nun den erwarteten Wert nehme, sollte ich Folgendes erhalten:

E [\frac{U (x)}{g (x)}] = \frac{1}{N} \approx \frac{1}{n} \sum_{i = 0}^{n} \frac{U (x)}{g (x)} .

$E\left [\frac{U(x)}{g(x)}\right ] = \frac{1}{N} \approx \frac{1}{n} \sum_{i=0}^n \frac{U(x)}{g(x)}.$

Ich habe versucht, dies in R für die Beispielfunktion testen $g(x) = e^{-x^2}$ . In diesem Fall benutze ich Metropolis-Hastings nicht, um die Samples zu generieren, sondern benutze die tatsächlichen Wahrscheinlichkeiten rnorm, um Samples zu generieren (nur zum Testen). Ich verstehe die gewünschten Ergebnisse nicht ganz. Grundsätzlich lautet der vollständige Ausdruck dessen, was ich berechnen würde:

\frac{1}{n (x_{max} - x_{min})} \sum_{i = 0}^{n} \frac{1}{e^{- x_{i}^{2}}} .

$\frac{1}{n(x_{\max} - x_\min)} \sum_{i=0}^n \frac{1}{ e^{-x_i^2}}.$ Dies sollte meiner Theorie nach

1 / \sqrt{π}

$1/\sqrt{\pi}$ . Es kommt näher, aber es konvergiert sicherlich nicht in der erwarteten Weise. Mache ich etwas falsch?

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

Bearbeiten Sie für CliffAB

Der Grund, warum ich den Bereich verwende, ist, dass ich auf einfache Weise eine Funktion definiere, die in dem Bereich, in dem sich meine Punkte befinden, ungleich Null ist, die jedoch im Bereich zu integriert wird . Die vollständige Spezifikation der Funktion lautet: Ich musste als einheitliche Dichte verwenden. Ich hätte eine andere Dichte verwenden können, die zu integriert ist , zum Beispiel die Wahrscheinlichkeitsdichte Dies hätte jedoch das Summieren der einzelnen Proben trivial gemacht, d.h. $1$ $[-\infty, \infty]$

U (x) = {\begin{cases} \frac{1}{x_{max} - x_{min}} & x_{max} > x > x_{min} \\ 0 & otherwise. \end{cases}

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

U (x)

$U(x)$

1

$1$

P (x) = \frac{1}{\sqrt{π}} e^{- x^{2}} .

$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{π}}{e^{- x_{i}^{2}}} = \frac{1}{n} \sum_{i = 0}^{n} \frac{1}{\sqrt{π}} = \frac{1}{\sqrt{π}} .

$\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}}.$

Ich könnte diese Technik für andere Distributionen ausprobieren, die in integriert sind . Ich möchte aber trotzdem wissen, warum es bei einer einheitlichen Verteilung nicht funktioniert. $1$

simulation monte-carlo metropolis-hastings numerical-integration Mike Flynn
quelle

Schauen Sie sich das nur kurz an, deshalb bin ich mir nicht sicher, warum Sie sich für range (x) entschieden haben. Vorausgesetzt, es ist gültig, ist es äußerst ineffizient! Der Bereich einer Stichprobe dieser Größe entspricht in etwa der instabilsten Statistik, die Sie erhalten können.

Cliff AB

@CliffAB Es gibt nichts Besonderes an mir, den Bereich zu verwenden, abgesehen davon, dass ich eine gleichmäßige Verteilung auf das Intervall definiere, in dem meine Punkte liegen. Siehe Bearbeitungen.

Mike Flynn

Ich werde später genauer darauf eingehen. Es ist jedoch zu berücksichtigen, dass x eine Menge einheitlicher RVs ist und dann als range . Aber wenn x eine Menge von nicht entarteten normalen Wohnmobilen ist, dann als , .

n \to \infty

$n \rightarrow \infty$

(x) \to 1

$(x) \rightarrow 1$

n \to \infty

$n \rightarrow \infty$

range (x) \to \infty

$\text{range}(x) \rightarrow \infty$

Cliff AB

@CliffAB Sie könnten recht gehabt haben, ich denke, der Grund war, dass die Grenzen des Integrals nicht festgelegt wurden, und so wird die Varianz des Schätzers nie konvergieren ...

Mike Flynn

Dies ist eine äußerst interessante Frage, die sich auf das Problem bezieht, eine Normalisierungskonstante einer Dichte auf der Grundlage eines MCMC-Ausgangssignals derselben Dichte anzunähern . (Eine Nebenbemerkung ist, dass die korrekte Annahme ist, dass integrierbar ist und es nicht ausreicht, im Unendlichen auf Null zu gehen.) $g$ $g$ $g$

Meiner Meinung nach ist der relevanteste Eintrag zu diesem Thema in Bezug auf Ihren Vorschlag ist , ein Papier von Gelfand und Dey (1994, JRSS B ), wo die Autoren einen sehr ähnlichen Ansatz entwickeln finden beim Erzeugen aus . Ein Ergebnis dieser Arbeit ist, dass für jede Wahrscheinlichkeitsdichte [dies ist äquivalent zu Ihrem ] die folgende Identität zeigt, dass eine Stichprobe aus eine

\int_{X} g (x) d x

$\int_\mathcal{X} g(x) \,\text{d}x$

p (x) \propto g (x)

$p(x)\propto g(x)$

α (x)

$\alpha(x)$

U (x)

$U(x)$

{x; α (x) > 0} \subset {x; g (x) > 0}

$\{x;\alpha(x)>0\}\subset\{x;g(x)>0\}$

\int_{X} \frac{α (x)}{g (x)} p (x) d x = \int_{X} \frac{α (x)}{N} d x = \frac{1}{N}

$\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

$p$ unvoreingenommene Bewertung von durch den Wichtigkeitsabtastungsschätzer Offensichtlich hängen die Leistungen (Konvergenzgeschwindigkeit, Existenz einer Varianz usw.) des Schätzers von der Wahl von [ obwohl seine Erwartung nicht]. In einem Bayes'schen Gerüst besteht eine von Gelfand und Dey befürwortete Wahl darin , die vorherige Dichte zu nehmen . Dies führt zu wobei die Wahrscheinlichkeitsfunktion ist, da

1 / N

$1/N$

\hat{η} = \frac{1}{n} \sum_{i = 1}^{n} \frac{α (x_{i})}{g (x_{i})} x_{i} \overset{iid}{\sim} p (x)

$\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{η}

$\hat\eta$

α

$\alpha$

α = π

$\alpha=\pi$

\frac{α (x)}{g (x)} = \frac{1}{ℓ (x)}

$\dfrac{\alpha(x)}{g(x)} = \dfrac{1}{\ell(x)}$

ℓ (x)

$\ell(x)$

g (x) = π (x) ℓ (x)

$g(x)=\pi(x)\ell(x)$ . Leider ist der resultierende Schätzer der harmonische Mittelwertschätzer , der auch der schlechteste Monte-Carlo-Schätzer aller Zeiten genannt wird Radford Neal von der University of Toronto. So klappt es nicht immer gut. Oder auch kaum jemals.

\hat{N} = \frac{n}{\sum_{i = 1}^{n} 1 / ℓ (x_{i})}

$\hat{N}=\dfrac{n}{\sum_{i=1}^n1\big/\ell(x_i)}$

Ihre Vorstellung, den Bereich Ihrer Stichprobe und die Uniform über diesen Bereich zu verwenden, hängt mit dem Problem des harmonischen Mittels zusammen: Dieser Schätzer hat keine Varianz, wenn auch nur wegen der erscheint im Zähler (ich vermute, dass dies für eine unbegrenzte Unterstützung immer der Fall sein könnte!) und konvergiert daher sehr langsam zur Normalisierungskonstante. Wenn Sie beispielsweise Ihren Code mehrmals wiederholen, erhalten Sie nach 10 after Iterationen sehr unterschiedliche Zahlenwerte. Dies bedeutet, dass Sie der Größe der Antwort nicht einmal vertrauen können. $(\min(x_i),\max(x_i))$ $\exp\{x^2\}$

Eine generische Lösung für dieses Problem der unendlichen Varianz besteht darin, für eine konzentriertere Dichte zu verwenden, beispielsweise unter Verwendung der Quartile Ihrer Stichprobe , weil dann bleibt über dieses Intervall begrenzt. $\alpha$ $(q_{.25}(x_i),q_{.75}(x_i))$ $g$

Wenn Sie Ihren Code an diese neue Dichte anpassen, ist die Annäherung viel näher an : $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

Wir diskutieren diese Methode ausführlich in zwei Artikeln mit Darren Wraith und mit Jean-Michel Marin .

Xi'an
quelle

Metropolis-Hastings-Integration - warum funktioniert meine Strategie nicht?

Antworten: