Einführung in Monte-Carlo-Methoden

Einführung

Zwei große Klassen numerischer Probleme, die in der Datenanalyse auftreten, sind Optimierungs- und Integrationsprobleme. Die zu einem Modell gehörenden Schätzer lassen sich nicht immer analytisch berechnen, daher greifen wir häufig zu numerischen Verfahren. Eine Möglichkeit, dieses Problem zu umgehen, ist die Simulation. Bei der Monte-Carlo-Schätzung werden hypothetische Ziehungen aus einer Wahrscheinlichkeitsverteilung simuliert, um wichtige Kenngrößen dieser Verteilung zu berechnen.

Die Grundidee der Monte-Carlo-Methode besteht darin, das Integral als Erwartungswert bezüglich einer geeigneten Wahrscheinlichkeitsverteilung zu schreiben und diesen anschließend mit dem Momentenmethode-Schätzer zu approximieren ($E[g(X)] \\approx \\overline{g(X)} = \\dfrac{1}{n}\\sum g(X_{i})$).

Haben wir eine stetige Funktion $g(\\theta)$ und wollen über das Intervall (a,b) integrieren, können wir das Integral als Erwartungswert einer gleichverteilten Zufallsvariablen $U \\sim U[a,b]$ umschreiben, also:

Mit dem Momentenmethode-Schätzer ergibt sich für die Integrationsnäherung:

Wobei die

simulierte Werte aus einer Gleichverteilung sind.

Beispiel 1: Exponentielle Integrationsnäherung

1) Gegeben sei die Funktion $f(x)= e^{x}$. Das Integral über das Intervall [3,5] ist:

Die Monte-Carlo-Näherung des Integrals ist:

#       Declaring the desired function
f = function(x){return(exp(x))}
#       Declaring the absolute error function
error = function(x,y){return(abs(x-y))}
#       The actual integral answer
ans = exp(5)-exp(3)

set.seed(6971)
#       number of iterations
n = 10^2
#       simulated uniform data
x= runif(n,3,5)
#       MonteCarlo approximation
MCa= (5-3)*mean(f(x))
#       Approximation error
e = error(ans,MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)
set.seed(6971)
for(k in 3:6){
  n = 10^k
  x = runif(n,3,5)
  mca = (5-3)*mean(f(x))
  rest = rbind(rest,c(n,mca,error(ans,mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = \"Integral Monte Carlo approximation results\",
      col.names =c(\"Number of simulations\",\"Monte Carlo approximation\",\"Error approximation\"))

Generalisierte Monte-Carlo-Näherung

Im allgemeinen Fall lautet die Integrationsnäherung für eine gegebene Verteilung f:

Ein Algorithmus zur Konstruktion von $\\widehat{I}$ lässt sich in folgenden Schritten beschreiben:

1) Erzeuge aus einer f-Verteilung.

2) Berechne:

3) Bestimme den Stichprobenmittelwert: $$\\overline{I} = \\dfrac{1}{n}\\sum_{k=1}^{n}\\dfrac{g(\\theta_k)}{f(\\theta_k)}$$

Im nächsten Abschnitt zeigen wir eine einfache Monte-Carlo-Funktion, um zu veranschaulichen, wie der Algorithmus arbeitet. Dabei sind a und b die Parameter der Gleichverteilung, n die Anzahl der gewünschten Simulationen und f die zu integrierende Funktion.

# The simple Monte Carlo function
MCaf = function(n,a,b,f){
x = runif(n,a,b)
MCa = (b-a)*mean(f(x))
return(MCa)
}

Monte-Carlo-Methoden in der Bayesianischen Datenanalyse

Die Grundidee der Bayesianischen Datenanalyse ist, ein Modell ( z. B. ein Regressions- oder Zeitreihenmodell ) mit einem Ansatz der Bayesianischen Inferenz zu schätzen. Wir nehmen an, dass unsere interessierenden Parameter eine theoretische Verteilung haben. Diese Verteilung (Posterior) wird mithilfe der Verteilung der beobachteten Daten (Likelihood) und der vorherigen bzw. externen Informationen über unsere Parameter (Prior-Verteilung) mittels des Satzes von Bayes aktualisiert.

$$p(\\theta / X) \\text{ } \\alpha\\text{ } p(X/\\theta)p(\\theta)$$ Dabei gilt:

$p(\\theta/X)$ ist die Posterior-Verteilung des Parameters.

$p(X/\\theta)$ ist die Stichprobenverteilung der beobachteten Daten ( Likelihood ).

$p(\\theta)$ ist die Prior-Verteilung des Parameters.

Das Hauptproblem im Bayesianischen Ansatz ist die Schätzung der Posterior-Verteilung. Markov-Chain-Monte-Carlo-Methoden ( mcmc ) erzeugen eine Stichprobe aus der Posterior-Verteilung und approximieren Erwartungswerte, Wahrscheinlichkeiten oder Quantile mithilfe von Monte-Carlo-Methoden.

In den nächsten beiden Abschnitten geben wir zwei Beispiele, um Wahrscheinlichkeiten und Quantile einer theoretischen Verteilung zu approximieren. Das oben skizzierte Vorgehen entspricht den gängigen Methoden im Bayesianischen Ansatz.

Beispiel 2: Wahrscheinlichkeitsnäherung einer Gammaverteilung

Angenommen, wir möchten die Wahrscheinlichkeit berechnen, dass eine Zufallsvariable $\\theta$ zwischen 0 und 5 liegt, also $P(0 < \\theta < 5)$, wobei $\\theta$ gammaverteilt ist mit den Parametern a = 2 und b = 1/3 ($\\theta \\sim Gamma(a = 2,b = 1/3)$). Dann gilt für die Wahrscheinlichkeit:

Wobei $I_{[0,5]}(\\theta) = 1$ ist, falls $\\theta$ im Intervall [0,5] liegt. Die Idee der Monte-Carlo-Näherung ist, die Anzahl der Beobachtungen zu zählen, die in [0,5] fallen, und diese durch die Gesamtzahl der simulierten Daten zu teilen.

set.seed(6972)
#       number of iterations
n = 10^2
#       simulated uniform data
x= rgamma(n,shape = 2,1/3)
#       MonteCarlo approximation
MCa= mean(x <= 5)
#       Approximation error
e = error(pgamma(5,2,1/3),MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)

for(k in 3:6){
  n = 10^k
  x= rgamma(n,shape = 2,1/3)
  mca= mean(x <= 5)
 rest = rbind(rest,c(n,mca,error(pgamma(5,2,1/3),mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = \"Probability Monte Carlo approximation results\",
      col.names =c(\"Number of simulations\",\"Monte Carlo approximation\",\"Error approximation\"))

Beispiel 3: Quantilsnäherung einer Normalverteilung

Angenommen, wir möchten das 0,95-Quantil einer Zufallsvariable $\\theta$ berechnen, die normalverteilt ist mit $\\mu = 20$ und $\\sigma = 3$ ($\\theta \\sim normal(\\mu = 20,\\sigma^{2} = 9)$). Dann lautet das 0,95-Quantil:

Die Grundidee ist, den größten Stichprobenwert zu finden, dessen Wahrscheinlichkeit höchstens 0,95 beträgt. Die Monte-Carlo-Quantilsnäherung wird geschätzt, indem wir das Quantil der simulierten Daten mit quantile() berechnen.

set.seed(6973)
#       number of iterations
n = 10^2
#       simulated uniform data
  x= rnorm(n,20,3)
#       MonteCarlo approximation
MCa= quantile(x,0.95)
#       Approximation error
e = error(qnorm(0.95,20,3),MCa)

rest =  data.frame(n = n,MCapprox = MCa,error = e)

for(k in 3:6){
  n = 10^k
  x= rnorm(n,20,3)
  mca= quantile(x,0.95)
 rest = rbind(rest,c(n,mca,error(qnorm(0.95,20,3),mca) ) )
}

kable(rest,digits = 5,align = 'c',caption = \"Quantile Monte Carlo approximation results\",
      col.names =c(\"Number of simulations\",\"Monte Carlo approximation\",\"Error approximation\"),row.names = FALSE)

Diskussion und Fazit

Monte-Carlo-Näherungen sind ein alternatives Werkzeug zur Integrationsnäherung und in der Bayesianischen Inferenz unverzichtbar, insbesondere bei anspruchsvollen, komplexen Modellen. Wie alle drei Beispiele zeigen, liefern Monte-Carlo-Methoden exzellente Approximationen, erfordern dafür jedoch sehr viele Simulationen, um den Approximationsfehler nahe null zu bringen.

Literatur

1) Introducing Monte Carlo methods with R, Springer 2004, Christian P. Robert and George Casella.

2) Handbook of Markov Chain Monte Carlo, Chapman and Hall, Steve Brooks, Andrew Gelman, Galin L. Jones, and Xiao-Li Meng.

3) Introduction to mathematical Statistics, Pearson, Robert V. Hogg, Joseph W. Mckean, and Allen T. Craig.

4) Statistical Inference An Integrated Approach, Chapman and Hall, Helio S. Migon, Dani Gamerman, Francisco Louzada.