Überblick

Teilnahme an der Lesesitzung "Einführung in die statistische Modellierung für die Datenanalyse" (sog. Midorimoto) des Unternehmens, Es war ein sehr leicht verständliches Buch, aber leider ist das Implementierungsbeispiel für Mac-Benutzer WinBUGS. Weil es ... war Ich habe den Bayes'schen Schätzungsansatz des verallgemeinerten linearen Modells in Kapitel 9 mit Python + STAN implementiert.

Was ich getan habe

Befolgen Sie ungefähr die folgenden Schritte.

Generieren Sie Dummy-Daten aus der Wahrscheinlichkeitsverteilung basierend auf bestimmten Parametern
Stellen Sie das Vorhersagemodell ein
Schätzen Sie aus den Dummy-Daten und dem Vorhersagemodell die Parameter (posteriore Verteilung), die die Daten mit MCMC generiert haben, und stimmen Sie mit den Antworten überein.

Insbesondere wird die Körpergröße einer bestimmten Pflanze (wobei diskrete Werte in Schritten von 0,1 von 3,0 bis 7,0 angenommen werden) als erklärende Variable verwendet. Schätzen Sie die Wahrscheinlichkeitsverteilung der Anzahl der Seeds (Ganzzahl 0 oder höher), die der Poisson-Verteilung folgt.

Verwendete Tools und Bibliotheken

MCMC-Sampler: STAN (PyStan)
Matrixberechnung, Array-Operation: NumPy
Grafikzeichnung: matplotlib

Schätzmethode

MCMC Die MCMC-Methode ist eine Abkürzung für die Marcov-Ketten-Monte-Carlo-Methode. Im Japanischen wird es als Markov-Ketten-Monte-Carlo-Methode bezeichnet. Markov-Kette Die Monte-Carlo-Methode ist eine Markov-Kette, die die Monte-Carlo-Methode verwendet. ... Ich bin eine Totologie geworden, deshalb erkläre ich etwas mehr.

Markov-Kette

Die Eigenschaft, dass ein bestimmter Zustand nur vom unmittelbar vorhergehenden Zustand abhängt, heißt ** Markov-Eigenschaft **. ** Markov-Kette ** ist ein Wahrscheinlichkeitsmodell, bei dem Zustände, die nur vom unmittelbar vorhergehenden Zustand abhängen, in einer Kette auftreten. Aus Sicht eines Ingenieurs ist es einfacher zu verstehen, wenn Sie es als einen ** Automatentyp ** betrachten.

Monte-Carlo-Methode

Die Monte-Carlo-Methode ist ein allgemeiner Name für numerische Berechnungen und Simulationen unter Verwendung von Zufallszahlen.

Mit anderen Worten, die Markov-Ketten-Monte-Carlo-Methode ** ist eine Methode zum Erzeugen einer Markov-Kette unter Verwendung von Zufallszahlen. Hier bezieht es sich insbesondere auf einen Algorithmus, der unter Verwendung der Eigenschaften der Markov-Kette eine Wahrscheinlichkeitsverteilung (posteriore Verteilung von Parametern) erzeugt.

Bayesianische Schätzung x MCMC

Durch Kombination des Bayes'schen Schätzrahmens (posteriore Verteilung ∝ Wahrscheinlichkeit x vorherige Verteilung) und MCMC, die die Wahrscheinlichkeitsverteilung proportional zur Wahrscheinlichkeit abtastet, Selbst wenn das Modell nicht analytisch gelöst werden kann, kann die posteriore Verteilung geschätzt werden, solange sie durch eine mathematische Formel ausgedrückt werden kann.

Implementierung

Eine detaillierte Erklärung des verallgemeinerten linearen Modells und der MCMC scheint sehr lang zu sein, daher werden wir mit der Implementierung beginnen. "Einführung in die statistische Modellierung für Datenanalyse-verallgemeinertes lineares Modell, hierarchisches Bayes-Modell, MCMC (Science of Probability and Information)" ist leicht zu verstehen.

Generierung von Trainingsdaten

Körpergröße ist 3,0 ~ 7,0, Durchschnitt μ = 1,5 + 0,1 * für jedes Individuum Es wird aus der Poisson-Verteilung der Körpergröße erzeugt.

def generate(n):
	for i in range(n):
		x = round(random.random() * 4 + 3, 1) # 3.0 ~ 7.Zufallszahl bis 0
		mu = math.exp(1.5 + 0.1*x)
		print (x, np.random.poisson(mu))

"Körpergröße" und "Anzahl der Samen".

STAN mcmc.stan

data {
  int<lower=1> N;
  vector[N] x;
  int<lower=0> y[N];
}
parameters {
  real beta1;
  real beta2;
}
model {
  for (i in 1:N) {
    y[i] ~ poisson(exp(beta1 + beta2 * x[i])); //Poisson-Verteilung x logarithmische Verknüpfungsfunktion
  }
  beta1 ~ normal(0, 1000); //Durchschnitt 0,Normalverteilung der Varianz 1000 ≒ keine Information vor Verteilung
  beta2 ~ normal(0, 1000); //Durchschnitt 0,Normalverteilung der Varianz 1000 ≒ keine Information vor Verteilung
}

Wie man liest

data Dies sind die Daten, die an das Stan-Programm übergeben werden sollen. Deklarieren Sie im Format ** {Datentyp} {Variablenname}; **. Übergeben Sie die Daten von Python, indem Sie den hier geschriebenen Variablennamen angeben.

parameters Dies ist die Variable, die in dem in stan beschriebenen Modell verwendet wird. Diesmal werden der Abschnitt β1 und der Koeffizient β2 der logarithmischen Verknüpfungsfunktion der Poisson-Verteilung innerhalb des STAN eingestellt. Es wird aus einer Nichtinformations-Vorverteilung (einer Normalverteilung mit einer sehr großen ungefähren Varianz) erzeugt.

model Es ist ein Vorhersagemodell. Die Operatoren ** '~' ** bedeuten, dass der linke Term der Wahrscheinlichkeitsverteilung des rechten Terms folgt. Hier folgt die Anzahl der Samen ** y ** einer Poisson-Verteilung mit ** exp (β1 + β2x) ** als Verknüpfungsfunktion (dh logarithmische Verknüpfungsfunktion).

Python

import numpy as np
import pystan
import matplotlib.pyplot as plt

data = np.loadtxt('./data.txt', delimiter=' ')

#Datengenerierung zur Weitergabe an Stan Interface
x = data[:,0] #Numpy-Notation, schneiden Sie die erste Datenspalte aus
y = data[:,1].astype(np.int) #Da es sich um die Zielvariable der Poisson-Verteilung handelt, wird sie in einen ganzzahligen Wert konvertiert.
N = data.shape[0] #Die Anzahl der Daten

stan_data = {'N': N, 'x': x, 'y': y}

fit = pystan.stan(file='./mcmc.stan',\
	data=stan_data, iter=5000, chains=3, thin=10)
    # iter =Nummer jeder Stichprobe
    # chain =Wiederholen Sie n durch iter angegebene Stichprobensätze
    # thin =Probenverdünnungsnummer

Ausführungsergebnis

Wenn es gut geht, wird ein Protokoll wie dieses herauskommen. Da STAN selbst in C ++ implementiert ist, wird die Kompilierung ausgeführt.

INFO:pystan:COMPILING THE C++ CODE FOR MODEL ~ NOW.
Chain 1, Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 0, Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 2, Iteration:    1 / 10000 [  0%]  (Warmup)
Chain 1, Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 2, Iteration: 1000 / 10000 [ 10%]  (Warmup)
Chain 0, Iteration: 1000 / 10000 [ 10%]  (Warmup)
...
Chain 0, Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 2, Iteration: 8000 / 10000 [ 80%]  (Sampling)
Chain 1, Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 0, Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 2, Iteration: 9000 / 10000 [ 90%]  (Sampling)
Chain 1, Iteration: 10000 / 10000 [100%]  (Sampling)# 
...
#  Elapsed Time: 9.51488 seconds (Warm-up)
#                10.7133 seconds (Sampling)
#                20.2282 seconds (Total)
#

Graph

Statistiken

summary

Inference for Stan model: anon_model_381b30a63720cfb3906aa9ce3e051d13.
3 chains, each with iter=10000; warmup=5000; thin=10; 
post-warmup draws per chain=500, total post-warmup draws=1500.

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
beta1   1.41  2.4e-3   0.05   1.31   1.38   1.41   1.45   1.52    481    1.0
beta2   0.12  4.6e-4   0.01    0.1   0.11   0.12   0.13   0.14    478    1.0
lp__  7821.4    0.04   0.97 7818.8 7821.1 7821.7 7822.1 7822.4    496    1.0

Samples were drawn using NUTS(diag_e) at Tue Feb  9 23:31:02 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Wenn die Methode fit.summary () ausgeführt wird, lautet die Ausgabe wie folgt. Ein Blick auf β1 zeigt zunächst, dass der Durchschnitt der Proben 1,41 beträgt und die Wahrscheinlichkeit, dass er im Bereich von 1,31 bis 1,52 liegt, bei 95% liegt. (In Bayes'schen Begriffen scheint es sich um eine Kreditabteilung zu handeln.)

Da es diesmal nur einen Berg in der Verteilung gibt, werde ich den repräsentativen Wert im Durchschnitt betrachten. β1 (1,5) ist 1,41 und β2 (0,1) ist 0,12, was sich nicht signifikant unterscheidet, und die numerischen Werte können geschätzt werden.

Rhat Beim Aufrufen von Stan's Code aus Python habe ich das Sampling dreimal mit dem Parameter * chain * wiederholt. Bei der Schätzung der posterioren Verteilung von Parametern durch MCMC werden die Anfangswerte von Parametern entsprechend bestimmt. Je nach Modell kann die geschätzte Wahrscheinlichkeitsverteilung je nach Anfangswert unterschiedlich sein. Wiederholen Sie die Abtastung mehrmals und prüfen Sie, ob die Wahrscheinlichkeitsverteilung um ** Rhat ** konvergiert, wodurch die Variation zwischen den Mengen quantifiziert wird. Es scheint in Ordnung zu sein, wenn es 1.1 oder weniger ist, aber dieses Mal kann beurteilt werden, dass es kein Problem gibt, da sowohl Beta1 als auch Beta2 kleiner sind.

Nachschlagewerk

Da der Schwerpunkt auf der Einführung der Implementierung liegt, ist die Bedeutung jedes Arguments wie * thin * beim Aufruf von STAN und Warum wird die erste Hälfte der Probenahme für * Warmup * verwendet? Darüber hinaus ist MCMC in erster Linie ein allgemeiner Begriff für Methoden, und was ist der spezifische Algorithmus? Wenn Sie mehr wissen möchten, lesen Sie bitte Folgendes.

[Einführung in die statistische Modellierung für das Datenanalyse-verallgemeinerte lineare Modell, Hierarchisches Bayes-Modell, MCMC (Wissenschaft der Wahrscheinlichkeit und Information)](http://www.amazon.co.jp/%E3%83%87%E3 % 83% BC% E3% 82% BF% E8% A7% A3% E6% 9E% 90% E3% 81% AE% E3% 81% 9F% E3% 82% 81% E3% 81% AE% E7% B5 % B1% E8% A8% 88% E3% 83% A2% E3% 83% 87% E3% 83% AA% E3% 83% B3% E3% 82% B0% E5% 85% A5% E9% 96% 80 % E2% 80% 95% E2% 80% 95% E4% B8% 80% E8% 88% AC% E5% 8C% 96% E7% B7% 9A% E5% BD% A2% E3% 83% A2% E3 % 83% 87% E3% 83% AB% E3% 83% BB% E9% 9A% 8E% E5% B1% A4% E3% 83% 99% E3% 82% A4% E3% 82% BA% E3% 83 % A2% E3% 83% 87% E3% 83% AB% E3% 83% BBMCMC-% E7% A2% BA% E7% 8E% 87% E3% 81% A8% E6% 83% 85% E5% A0% B1% E3% 81% AE% E7% A7% 91% E5% AD% A6-% E4% B9% 85% E4% BF% 9D-% E6% 8B% 93% E5% BC% A5 / dp / 400006973X)

[PYTHON] Die erste Markov-Ketten-Monte-Carlo-Methode von PyStan