Lösen einer einfachen Integralgleichung durch Zufallsstichprobe

8

Sei eine nichtnegative Funktion. Ich bin daran interessiert, so zu finden, dass Die Einschränkung : Alles was ich tun kann, ist an Punkten in . Ich kann jedoch die Orte, an denen ich zufällig probiere, nach Belieben auswählen . fz[0,1]

0zf(x)dx=1201f(x)dx
f[0,1]f

Fragen:

  1. Ist es möglich, nach endlich vielen Stichproben eine unvoreingenommene Schätzung von zu erhalten ? Wenn ja, was ist die kleinstmögliche Varianz einer solchen Schätzung nach Stichproben?zk
  2. Wenn nicht, welche Verfahren stehen zur Schätzung von zur Verfügung und welche Konvergenzzeiten sind damit verbunden?z

Wie Douglas Zare in den Kommentaren hervorhob, kann dies sehr schwierig sein, wenn die Funktion nahe Null oder sehr groß ist. Glücklicherweise ist die Funktion, für die ich dies verwenden muss, von oben und unten begrenzt. Nehmen wir also an, dass . Darüber hinaus können wir auch annehmen, dass Lipschitz ist oder sogar differenzierbar, wenn dies hilft.1f(x)2f

Robinson
quelle
1
Wenn Sie keine weiteren Informationen haben, können Sie sich sehr schlecht verhalten. Stellen Sie sich vor , dass ist zwischen und , undDurch geringfügige Änderungen an springt der Median von unter auf über . 0 1 / 3 2 / 3 1 / 3 0 f ( x ) d x 1 / 2 f 1 / 3 2 / 3f01/32/301/3f(x) dx1/2.f1/32/3
Douglas Zare
@robinson Könnten Sie weitere Informationen zu bereitstellen ? Oder sind Sie daran interessiert, das Problem für eine beliebige Dichte zu lösen ? fff
@DouglasZare - Danke für den Kommentar; siehe meine Bearbeitung.
Robinson
@Procrastinator - Ich habe die Frage mit etwas mehr Informationen bearbeitet.
Robinson
3
(+1) Für das Update. Wenn man die linke Seite durch die rechte Seite teilt, kann man sehen, dass sich dies darauf reduziert, den Median einer unbekannten Wahrscheinlichkeitsverteilung zu finden, die auf . [0,1]
Kardinal

Antworten:

6

Wie Kardinal in seinem Kommentar ausgeführt hat, kann Ihre Frage wie folgt angepasst werden.

Durch einfache Algebra kann die Integralgleichung als umgeschrieben werden wobei die Wahrscheinlichkeitsdichtefunktion ist, die als G g ( x ) = f ( x )

0zg(x)dx=12,
g
g(x)=f(x)01f(t)dt.

Sei eine Zufallsvariable mit der Dichte . Per Definition ist , sodass Ihre Integralgleichung was bedeutet, dass Ihr Problem wie folgt angegeben werden kann:XgP{Xz}=0zg(x)dx

P{Xz}=12,

"Sei eine Zufallsvariable mit der Dichte . Finde den Median von "XgX

Um den Median von zu schätzen , verwenden Sie eine beliebige Simulationsmethode, um eine Stichprobe von Werten zu zeichnen und den Stichprobenmedian als Schätzung zu verwenden.XX

Eine Möglichkeit besteht darin, den Metropolis-Hastings-Algorithmus zu verwenden, um eine Stichprobe von Punkten mit der gewünschten Verteilung zu erhalten. Aufgrund des Ausdrucks der Akzeptanzwahrscheinlichkeit im Metropolis-Hastings-Algorithmus müssen wir den Wert der Normalisierungskonstante der Dichte . Wir müssen diese Integration also nicht durchführen.01f(t)dtg

Der folgende Code verwendet eine besonders einfache Form des Metropolis-Hastings-Algorithmus, der als Indepence Sampler bekannt ist und einen Vorschlag verwendet, dessen Verteilung nicht vom aktuellen Wert der Kette abhängt. Ich habe unabhängige einheitliche Vorschläge verwendet. Zum Vergleich gibt das Skript das Monte-Carlo-Minimum und das mit der Standardoptimierung gefundene Ergebnis aus. Die Abtastpunkte werden im Vektor gespeichert chain, aber wir verwerfen die ersten Punkte, die die sogenannte "Einbrenn" -Periode der Simulation bilden.10000

BURN_IN = 10000
DRAWS   = 100000

f = function(x) exp(sin(x))

chain = numeric(BURN_IN + DRAWS)

x = 1/2

for (i in 1:(BURN_IN + DRAWS)) {
    y = runif(1) # proposal
    if (runif(1) < min(1, f(y)/f(x))) x = y
    chain[i] = x
}

x_min = median(chain[BURN_IN : (BURN_IN + DRAWS)])

cat("Metropolis minimum found at", x_min, "\n\n")

# MONTE CARLO ENDS HERE. The integrations bellow are just to check the results.

A = integrate(f, 0, 1)$value

F = function(x) (abs(integrate(f, 0, x)$value - A/2))

cat("Optimize minimum found at", optimize(F, c(0, 1))$minimum, "\n")

Hier sind einige Ergebnisse:

Metropolis minimum found at 0.6005409 
Optimize minimum found at 0.601365

Dieser Code ist nur als Ausgangspunkt für das gedacht, was Sie wirklich brauchen. Daher mit Vorsicht verwenden.

Zen
quelle
Danke für deine Antwort. Ich kenne R nicht, daher bin ich mir nicht sicher, wie ich analysieren soll, was Sie tun. Könnten Sie Ihre Prozedur in Worten / Formeln angeben? Vielen Dank. Insbesondere frage ich mich, ob Sie die Einschränkung respektieren, dass das einzige, was Sie tun können, die Bewertung von f ist - Sie dürfen beispielsweise nicht integrieren (obwohl Sie sicherlich Monte-Carlo-Näherungen an Integrale basierend auf bilden können zufällige Bewertungen). f
Robinson
Ja, ich bewerte nur , um die Monte-Carlo-Schätzung zu erhalten. f
Zen
Der Code ist nur ein Beispiel. Die R-Syntax ähnelt anderen Sprachen. Eine bestimmte Aussage, die Sie nicht verstehen? Schauen Sie sich die Wikipedia-Seite zum Metropolis-Hastings-Algorithmus an. Natürlich ist die allgemeine Idee wichtiger. Sie können mit jeder verfügbaren Methode aus dem probieren . f/f
Zen
Haben Sie einen Einführungskurs in stochastische Prozesse besucht, der zeitdiskrete Markov-Ketten abdeckt?
Zen
1
Übrigens: Zauderer der Welt, vereinigt euch! Aber nicht heute ...
Zen
3

Die Qualität der integralen Approximation, zumindest im Fall von 1D, ist gegeben durch (Satz 2.10 in Niederreiter (1992) ): wobei ist der Kontinuitätsmodul der Funktion (bezogen auf die und für Lipshitz-Funktionen leicht ) und ist die (extreme) Diskrepanz oder die maximale Differenz zwischen dem Anteil der Treffer durch die Sequenz

|1Nn=1Nf(xn)01f(u)du|ω(f;DN(x1,,xN))
ω(f;t)=sup{|f(u)f(v)|:u,v[0,1],|uv|t,t>0}
DN(x1,,xN)=supu|1Nn1{xn[0,u)}u|=12N+maxn|xn2n12N|
x1,,xNeines halboffenen Intervalls und seines Lebesgue-Maßes . Der erste Ausdruck ist die Definition, und der zweite Ausdruck ist die Eigenschaft der 1D-Sequenzen in (Satz 2.6 im selben Buch).[0,u)u[0,1]

Um den Fehler in der integralen Approximation, zumindest in der rechten Seite Ihrer Gleichung, zu minimieren, müssen Sie offensichtlich . Scheiß auf die zufälligen Auswertungen, sie laufen Gefahr, eine zufällige Lücke bei einem wichtigen Merkmal der Funktion zu haben.xn=(2n1)/2N

Ein großer Nachteil dieses Ansatzes besteht darin, dass Sie sich auf einen Wert von , um diese gleichmäßig verteilte Sequenz zu erzeugen. Wenn Sie mit der Qualität der Annäherung nicht zufrieden sind, können Sie lediglich den Wert von verdoppeln und alle Mittelpunkte der zuvor erstellten Intervalle erreichen.NN

Wenn Sie eine Lösung suchen, mit der Sie die Anzahl der Punkte schrittweise erhöhen können, können Sie dieses Buch weiter lesen und sich über Van-der-Corput-Sequenzen und radikale Umkehrungen informieren. Siehe Sequenzen mit geringer Diskrepanz auf Wikipedia, es enthält alle Details.

Update: nach zu lösen , definieren Sie die Teilsumme Finden Sie so, dass und interpolieren Sie, um Diese Interpolation setzt voraus, dass stetig ist. Wenn zusätzlich zweimal differenzierbar ist, dann diese Näherung durch Integrieren der Erweiterung zweiter Ordnung, um und , und Lösen einer kubischen Gleichung für .z

Sk=1Nn=1kf(2n12N).
k
Sk12SN<Sk+1,
f()f()Sk-1Sk+2z
zN=2k12N+SN/2SkN(Sk+1Sk).
f()f()Sk1Sk+2z
StasK
quelle
1
Ich mag den Kern davon. Ich denke, es wäre hilfreich, die Strategie, die Sie zur Lösung der OP-Frage vorschlagen, deutlicher zu machen. Gegenwärtig lautet die Antwort (für mich) meistens so, als würde sie sich mit der Berechnung der RHS der Gleichung in der Frage befassen.
Kardinal
1
(+1) Nettes Update. kann einfach als Riemann-Summen-Näherung an das Integral angesehen werden, wobei wir den Wert von am Mittelpunkt jedes durch die Partition definierten Intervalls anstelle des linken oder rechten Endpunkts verwenden. :-) fSNf
Kardinal
Ja; Es ist jedoch interessant, dass diese Riemannsche Summe diese Optimalitätsbegründung hat.
StasK