Einführung

Infolge der Erfassung verschiedener Beobachtungsdaten besteht möglicherweise der Wunsch, die Wahrscheinlichkeitsverteilung hinter dem Beobachtungswert auf das Modell anzuwenden und zu schätzen. Es ist praktisch für die Berechnung verschiedener Statistiken und Simulationsberechnungen. Es gibt Normalverteilung, Exponentialverteilung, Binomialverteilung, Beta-Verteilung usw. als Modelle der Wahrscheinlichkeitsverteilung, diesmal jedoch, während die MCMC-Methode (Markov-Ketten-Monte-Carlo-Methode) "einfach" aus dem Prinzip der Gammaverteilung implementiert wird. Ich würde gerne versuchen, wie viel es montiert werden kann. Als gutes Tool zum Ausführen von MCMC können WinBUGS und PyMC3 Es gibt //docs.pymc.io/) usw., so dass Sie es nicht selbst machen müssen **. Es gibt jedoch viele Dinge, die Sie sehen können, wenn Sie es selbst machen, und ich denke, es ist ein Verdienst, das Verständnis zu fördern.

Was ist Gammaverteilung?

Was ist Gamma-Verteilung? Es ist eine Verteilung, die durch zwei Parameter gekennzeichnet ist: Formparameter $ k $ und Skalierungsparameter $ q $.

P(x) = \frac{1}{\Gamma(k)q^k}x^{k-1}e^{-x/q},\ x \geq 0 \ (k>0,\ q>0)

Es ist eine Verteilung, die nicht negative reelle Zahlen verwendet, und ein ausgezeichneter Typ, der auch Exponentialverteilung, Chi-Quadrat-Verteilung oder Alan-Verteilung ausdrücken kann. Besonders gefällt mir die Tatsache, dass es eine verzerrte Verteilung mit einem breiten Saum ausdrücken kann. Es gibt viele Beobachtungsdaten auf der Welt, die durch nicht negative Werte wie Anzahl, Länge, Gewicht, Zeitintervall, Bevölkerung, Geldbetrag usw. dargestellt werden, und sie scheinen nützlich zu sein, um die Verteilung solcher Variablen auszudrücken. Tatsächlich scheint es in den Ingenieur- und Wirtschaftswissenschaften wie der Wartezeitverteilung und der Einkommensverteilung verwendet zu werden.

Was ist MCMC?

Die MCMC-Methode (Markov-Chain-Monte-Carlo-Methode) ist eine Methode, die effizient eine Stichprobe einer bestimmten Wahrscheinlichkeitsverteilung generiert und die Art der Wahrscheinlichkeitsverteilung aus der Stichprobe schätzt. Die Markov-Kette ist für die Erzeugung der Stichprobe gemäß der Wahrscheinlichkeitsverteilung verantwortlich, und Monte Carlo ist für die Schätzung aus der Stichprobe verantwortlich. Für Details ist das Video des Kommentars von Dr. Yukito Iba des Instituts für Statistische Mathematik sehr leicht zu verstehen. In diesem Artikel verwenden wir MCMC, um die Gammaverteilung zu schätzen, und definieren daher die folgenden Begriffe:

$ \ Theta $: $ N $ Dimensionsparameter, den Sie schätzen möchten
$ P (x | \ theta) $: Wahrscheinlichkeitsverteilungsmodell berechnet nach Parameter $ \ theta $ (ich weiß es, weil es ein Modell ist)
$ X_1, X_2, ..., X_n $: Beobachteter Stichprobenwert. Unabhängigkeit wird jedoch vorausgesetzt.
$ D_n = (X_1, X_2, ..., X_n) $: Beobachtete Daten (Wertesatz)
$ P (\ theta | D_n) $: Wahrscheinlichkeitsverteilung des Parameters $ \ theta $, geschätzt aus den Daten $ D_n $

Was ich hier tun möchte, ist, viele Daten $ D_n $ zu sammeln, und wenn die Wahrscheinlichkeitsverteilung $ P (\ theta | D_n) $ des Parameters $ \ theta $ gut berechnet werden kann, kann der Parameter $ \ theta $ des angenommenen Modells ( Dies bedeutet, dass Sie den erwarteten Wert berechnen können (basierend auf der Ex-post-Wahrscheinlichkeit). Die Wahrscheinlichkeitsverteilung $ P (\ theta | D_n) $ wird durch die folgende Gleichung aus dem Bayes'schen Theorem (oder aus der Normalisierungsbedingung der Wahrscheinlichkeitsverteilung) ausgedrückt.

P(\theta|D_n) = P(D_n | \theta)P(\theta) / \int_{\theta} d\theta P(D_n | \theta)P(\theta)

Hier sollte die Integration für alle $ \ theta $ mit einer Wahrscheinlichkeit ungleich Null für die Erzeugung der Daten $ D_n $ durchgeführt werden. Eine der einfachen Möglichkeiten, diese Verteilung zu schätzen, ist die Gitterberechnung. Mit anderen Worten, wenn $ \ theta $ gleichmäßig auf einem geeigneten Gitter verteilt ist und durch den diskreten Wert $ \ theta_i (i = 1,2, ..., I) $ dargestellt wird,

P(\theta_i|D_n) \simeq P(D_n | \theta_i)P(\theta_i) / \sum_{j \in I}  P(D_n | \theta_j)P(\theta_j)

Das Integral des Nenners kann wie in durch die Summe angenähert werden. Es sieht so aus, wenn ich es in ein Diagramm schreibe.

Um jedoch die Schätzgenauigkeit von $ \ theta $ zu verbessern, muss die Anzahl der diskreten Werte $ I $ groß genug sein. Wenn Sie einfach versuchen, einen $ N $ -Dimensionsparameter mit $ M $ -Netzen darzustellen, benötigen Sie diskrete Werte für $ I = M ^ N $, die mit zunehmender Dimension explodieren.

Also ändert ** MCMC die Denkweise und legt den diskreten Wert $ \ theta_i $ nicht fest, sondern bewegt ihn stetig **. Verschiebe $ k (k = 1,2, ...) $, $ k $ th $ \ theta_i (i = 1, .., I) $ nach $ \ theta_i (k) $, $ \ theta_i ( Sei $ P_k (\ theta) $ die Wahrscheinlichkeitsverteilung, der k) $ folgt. Lassen Sie $ \ pi (x \ rightarrow x ') $ von $ \ theta_i (k) $ nach $ \ theta_i (k + 1) $ wechseln. Dann kann es wie folgt dargestellt werden. P_k_k+1_map1.jpg

Ich habe diese Art des Verschiebens von $ \ pi (x \ rightarrow x ') $ gut entworfen, und nachdem ich es viele Male verschoben habe,

P_k(\theta) \rightarrow P(\theta | D_n) \  ({\rm as \ } k \rightarrow \infty)

Wenn dies zutrifft, können Sie Ihr Ziel erreichen. Als nächstes betrachten wir den Metropolis-Algorithmus als eine der Möglichkeiten, ihn zu verschieben.

Metropolis-Algorithmus

Der Metropolis-Algorithmus ist ein Stichprobenalgorithmus aus einer Wahrscheinlichkeitsverteilung. In der vorherigen Abbildung erhöhen wir die Anzahl der Zeichen.

$ P_ * (\ theta) : Verteilung gemäß dem Modell, das Sie schätzen möchten ( P (\ theta | D_n) $ in diesem Artikel)
$ P_ {o_k} (\ theta) $: Verteilung nach anderen als $ P_ * (\ theta) $ in $ P_k (\ theta) $

Mit anderen Worten, nehmen wir an, $ \ theta_i $ wird durch die Synthese von Verteilungen generiert, die dem Modell folgen und nicht folgen. Nehmen wir zu diesem Zeitpunkt an, dass die Bewegung durch Zuordnen von $ \ pi (x \ rightarrow x ') $ eingeschränkt ist, wie in der folgenden Abbildung gezeigt.

P_k_k+1_map2.jpg

Mit anderen Worten, die Einschränkung besteht darin, dass $ \ theta_i $, das auf ** $ P_ * (\ theta) $ folgte, immer zu $ P_ * (\ theta) $ ** ((1) in der Abbildung) verschoben wird. Nehmen wir außerdem an, dass $ \ theta_i $ nach ** $ P_ {o_k} (\ theta) $ immer nach $ P_ * (\ theta) $ ** (② in der Abbildung) verschoben werden kann. Wenn dann ② mäßig existiert, kann aus der Begrenzung von can gesehen werden, dass mit zunehmendem $ k $ $ P_k (\ theta) \ rightarrow P_ * (\ theta) $ wird. Genau genommen wird (1) als Beständigkeit bezeichnet, (2) als Kontraktivität, und es scheint, dass die Konvergenz der Wahrscheinlichkeitsverteilung aus diesen beiden Bedingungen als Ergodismus bezeichnet wird. (Einzelheiten finden Sie unter Kommentarvideo von Dr. Yukito Iba).

Eine der Zuordnungen $ \ pi (x \ rightarrow x ') $, die diese Eigenschaften erfüllen, ist der Metropolis-Algorithmus, der durch den folgenden Fluss dargestellt wird.

Generieren Sie $ x '$ aus der symmetrischen Wahrscheinlichkeitsverteilung $ Q (x, x') (= Q (x ', x)) $ mit dem aktuellen Punkt $ x $.
Generieren Sie eine einheitliche Zufallszahl $ 0 \ leq r <1 $.
Wenn $ r <{P_ * (x ')} / {P_ * (x)} $, gehen Sie zu $ x \ rightarrow x' . Wenn nicht, bewegen Sie sich nicht ( x \ rightarrow x $).
Wiederholen Sie die Schritte 1 bis 3.

Solange beispielsweise die Wahrscheinlichkeitsverteilung von 1 symmetrisch ist, kann eine Normalverteilung verwendet werden.

Q(x,x') = \frac{1}{\sqrt{2\pi \sigma^2}}exp\left(-\frac{(x-x')^2}{2\sigma^2}\right)

Übrigens wird der auf den nicht symmetrischen Fall ($ Q (x, x ') \ neq Q (x', x) $) erweiterte Algorithmus als Metropolis-Hastings-Algorithmus bezeichnet. Betrachtet man Schritt 3,

$ P_ * (x ') \ geq P_ * (x) $ bewegt sich immer von x nach x'. Mit anderen Worten, wir werden zu dem mit der höheren Wahrscheinlichkeit übergehen. (Entspricht ② in der obigen Abbildung)
$ P_ * (x ') <P_ * (x) $ bewegt sich nach x'oder x, abhängig vom Verhältnis der Wahrscheinlichkeit $ P_ * (x') $ zu $ P_ * (x) $. (Entspricht ① in der obigen Abbildung)

Versuchen Sie mit Python zu rechnen

Lassen Sie uns nun die Gammaverteilung bei der Implementierung von MCMC in Python schätzen.

1. Vorbereitung

Stellen Sie die Bibliothek ein.

import numpy as np
import matplotlib.pyplot as plt
from math import gamma

Definiert eine Gammafunktion, eine Histogrammgenerierung und eine Funktion, die die Form der Wahrscheinlichkeitsverteilung anzeigt.

def gamma_pdf(x, k, q):
    p = 1./(gamma(k) * pow(q, k)) * pow(x, k-1) * np.exp( -x/q )
    return p

def genPdfHist(f, xmin, xmax, N):
    x = np.linspace(xmin,xmax, N)
    p = [f(v) for v in x]
    return x,p

def showPdf(f, xmin, xmax, N, yscale=1.0):
    x, p = genPdfHist(f, xmin, xmax, N)
    y = [e*yscale for e in p]
    plt.plot(x,y)
    plt.show()

2. Bereiten Sie Beobachtungsdaten vor

Als nächstes erstellen wir Beobachtungsdaten unter Verwendung der Gammaverteilung. Definieren Sie eine abzutastende Funktion zum Erstellen von Beobachtungsdaten.

def PlainSampler(pdf, xmin, xmax, N):
    n = 0
    xs = []
    while n < N:
        x = xmin + np.random.rand()*(xmax-xmin)
        p = pdf(x)
        if np.random.rand() < p:
            # accept
            xs.append(x)
            n = n + 1
    return xs

Verwenden Sie die obige Funktion, um Beobachtungsdaten zu generieren. (Bedeutung experimentelle Beobachtungsdaten) Generieren wir zunächst das richtige Modell mit den Parametern $ (k ^ *, q ^ *) = (2,4) $.

k_true = 2.0
q_true = 4.0
gp1 = lambda x : gamma_pdf(x, k_true, q_true)

smp1 = PlainSampler(gp1, 0, 50, 1000)
h1 = plt.hist(smp1, bins=50)
showPdf(gp1, 0, 50, 100, h1[0].sum()/(h1[1][1]-h1[1][0]))
plt.show()

Die Zeile der Funktion showPdf wird mit dem Anpassungsfaktor multipliziert, um die Y-Achse des Histogramms an der Wahrscheinlichkeitsverteilung auszurichten. Die Daten smp1 sind der beobachtete Wert und das Histogramm sieht folgendermaßen aus. Da es mit Zufallszahlen erzeugt wird, ändert sich natürlich die Form des Histogramms jedes Mal, wenn es ausgeführt wird. pdf_sample_k2.0_q4.0.png

3. Bereiten Sie eine Funktion für MCMC vor

Wir werden die Funktionen definieren, die zur Berechnung der MCMC verwendet werden. Erstens die posteriore WahrscheinlichkeitsverteilungP(D_n |\theta)Ist eine Funktion, die berechnet. Gleichzeitige Wahrscheinlichkeit durch weitere Auswahl von 50 Stichprobenwerten aus der ursprünglichen Stichprobe, um zu verhindern, dass der berechnete Wahrscheinlichkeitswert zu klein wird, und um die Datenunabhängigkeit sicherzustellenP(D_n |\theta)=P(X_1 | \theta)P(X_2 | \theta) \cdots P(X_n | \theta)Wird berechnet.

def pdf_gamma_sample(x, sample):
    k = x[0]
    q = x[1]
    #Verhindern Sie, dass P zu klein wird
    sample1 = np.random.choice(sample, 50)
    p = [gamma_pdf(s, k, q) for s in sample1]
    return np.prod(p)

Eine Funktion, die eine vorherige Verteilung (Anfangsverteilung) des Parameters $ \ theta = (k, q) $ ergibt. Da es keine vorherigen Informationen gibt, generieren wir sie aus einer einheitlichen Zufallszahl. Ich werde jedoch nur den Maximalwert angeben.

def genInit( cdim, xdim, xmax):
    x = np.random.rand(cdim, xdim)*xmax
    return x

Als nächstes kommt der Kernteil. Eine Funktion, die MCMC mit einem Mapping $ \ pi (x \ rightarrow x ') $ durch den Metropolis-Algorithmus ausführt.

def MHmap(pdf, x, xdim, sig):
    #Gibbs-Abtastung (nur eine Variable ändern)
    xnew = x + np.eye(1, xdim, np.random.randint(xdim))[0,:] * np.random.normal(loc=0, scale=sig)
    #Begrenzen Sie den Wertebereich von x
    xnew = np.maximum(xnew, 0.01)
    #Ausgewählt nach Metropolis-Kriterien
    r = np.random.rand()
    if r < pdf(xnew)/pdf(x):
        x = xnew
    return x

def MCMC_gamma_01(sample, x0, sig, N):
    #
    pdf1 = lambda x: pdf_gamma_sample(x, sample)
    #
    x = x0
    xs = []
    for i in range(N):
        #Zuordnungspunkte gemäß Metropolis-Standards
        x = [MHmap(pdf1, e, 2, sig) for e in x]
        #Zufälliges Resampling
        s = np.random.randint(0,len(x),(len(x)))
        x = [x[j] for j in s]
        #
        xs.append(x)
    return np.array(xs)

Dies ist der Punkt beim Codieren.

Gibbs-Sampling verwenden (möglicherweise nicht erforderlich)
Die Übergangsfunktion $ Q (x, x ') $ verwendet eine Normalverteilung mit der Standardabweichung $ \ sigma $
Begrenzen Sie den Bereich der Variablen x (springen Sie nicht in die negative Domäne) Nach dem Mapping mit + π erneut abtasten (theoretisch sollte es nicht notwendig sein, aber ich möchte die Möglichkeit verringern, dass zu weit fliegende Punkte nicht leicht zurückkommen).

** Ich bin nicht sicher, ob diese Implementierung optimal ist **, aber ich kann mir vorstellen, dass die tatsächlichen Bibliotheken (PyMC3, WinBUGS usw.) auf verschiedene Arten entwickelt worden wären. Hier wird die Lücke zwischen Theorie und Umsetzung sichtbar.

4. Versuchen Sie es mit MCMC

Lassen Sie uns nun MCMC mit der in 3 erstellten Funktion verwenden, um den Parameter $ \ theta = (k, q) $ des Gammaverteilungsmodells aus den in 2 erstellten Beobachtungsdaten smp1 zu schätzen. Stellen Sie die Größe der Punktgruppe auf 500, die Dimension des Parameters auf 2 Dimensionen und den Maximalwert auf 10 ein. Außerdem sei die Standardabweichung der Übergangsfunktion $ Q (x, x ') $ 3. Führen Sie die MCMC-Berechnung 100 Mal durch.

p0 = genInit( 500, 2, 10)
sig = 3
xs = MCMC_gamma_01(smp1, p0, sig, 100)

Nach einigen Minuten wird das Ergebnis angezeigt. Erstellen wir also ein Diagramm.

fig = plt.figure()
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)
ax1.hist(xs[-1,:,0], bins=50)
ax2.hist(xs[-1,:,1], bins=50)
ax1.set_xlabel("k")
ax2.set_xlabel("q")
k_est = xs[-1,:,0].mean()
q_est = xs[-1,:,1].mean()
print("k*:{}, q*:{}".format(k_est, q_est))
plt.show()

Das Ergebnis ist unten dargestellt. Es ist kein sehr sauberes Histogramm, aber wenn es mit dem Durchschnittswert berechnet wird, ist es $ (\ bar {k}, \ bar {q}) = (2.21, 4.46) $. MCMC_sample_k2.0_q4.0.png

Ein Vergleich der ursprünglich generierten Verteilung und der Verteilung mit $ (\ bar {k}, \ bar {q}) = (2.21, 4.46) $ ist in der folgenden Abbildung dargestellt. Es ist ein wenig anders, aber es scheint, dass die Long-Tail-Verteilung bis zu einem gewissen Grad angepasst werden kann. MCMC_match_k2.0_q4.0.png

5. Andere Ergebnisse

Wenn Sie versuchen, den Parameter $ \ theta = (k, q) $ auf verschiedene Weise zu ändern, scheint er zu passen. Der Eindruck ist, dass es konvergiert, nachdem es ungefähr 100 Mal wiederholt wurde.

Erwägung

Aus dem Obigen wurden die folgenden Trends aus den Ergebnissen der Schätzung der Gammaverteilung durch einfaches Implementieren von MCMC gesehen.

Es ist (bis zu einem gewissen Grad) möglich, die Parameter der ursprünglichen Gammaverteilung aus dem Probenmittelwert unter Verwendung der aus der Gammaverteilung erzeugten Beobachtungsdaten abzuschätzen.
Es gab jedoch eine Tendenz, von der Anfangsverteilung als der vorherigen Verteilung, die MCMC gegeben wurde, und der Übergangsfunktion $ Q (x, x ') $ (Dispersion) ** abzuhängen.
Es scheint nützlich zu sein, da es für eine Funktion mit starker Nichtlinearität verwendet werden kann, die als Gammaverteilung bezeichnet wird, und im Prinzip auch dann verwendet werden kann, wenn die Parameterdimension hoch ist.

Außerdem ...

Der Teil, der das Wahrscheinlichkeitsverhältnis basierend auf dem Modell aus den Daten $ D_n $ berechnet, wurde sehr einfach geschrieben, aber es scheint, dass wir verschiedene Kennzahlen entwickeln können, z. B. Kennzahlen, wenn die Wahrscheinlichkeit sehr klein wird und sich beschleunigt.
Die Schätzung der Wahrscheinlichkeitsverteilung hinter dem als Gruppe von Punkten gestreuten Probensatz scheint den Partikel- und Welleneigenschaften von Licht etwas ähnlich zu sein, und ich glaube, dass sie zusammenhängen.
Was passiert, wenn sich die Wahrscheinlichkeitsverteilung im Laufe der Zeit ändert? Ist auch interessant.

Referenzlink

Ich habe auf die folgende Seite verwiesen.

[PYTHON] Versuchen Sie, die Parameter der Gammaverteilung zu schätzen, während Sie einfach MCMC implementieren