Einführung

Dieser Artikel befasst sich mit der Bayes'schen Schätzung eines einfachen linearen Regressionsmodells. Insbesondere wurde die Markov-Ketten-Monte-Carlo-Methode (allgemein bekannt als MCMC) in Python unter Verwendung von Gibbs-Sampling implementiert. Die Ausführung von MCMC in Python ist berühmt für das Modul pymc, aber dieses Mal habe ich es gewagt, es nur mit numpy und scipy zu implementieren, ohne pymc zu verwenden. Beim Schreiben dieses Artikels schrieb Herr Kosumi ["Bay's Computational Statistics" (Asakura Shoten)](https://www.amazon.co.jp/%E3%83%99%E3%82%A4%E3%82 % BA% E8% A8% 88% E7% AE% 97% E7% B5% B1% E8% A8% 88% E5% AD% A6-% E7% B5% B1% E8% A8% 88% E8% A7% A3% E6% 9E% 90% E3% 82% B9% E3% 82% BF% E3% 83% B3% E3% 83% 80% E3% 83% BC% E3% 83% 89-% E5% 8F% A4 % E6% BE% 84-% E8% 8B% B1% E7% 94% B7 / dp / 4254128568 / ref = sr_1_1? Ie = UTF8 & qid = 1485493617 & sr = 8-1 & keywords =% E3% 83% 99% E3% 82% A4 % E3% 82% BA% E8% A8% 88% E7% AE% 97% E7% B5% B1% E8% A8% 88) wurde stark erwähnt.

Theorie

Lassen Sie uns vor der Implementierung auf die Theorie zurückblicken, obwohl sie einfach ist. Das lineare Regressionsmodell, das Sie schätzen möchten

y_i=\boldsymbol{x}_{i}'\boldsymbol{\beta}+u_i, \ \ u_i\sim N(0,\sigma^2)　\ \  (i=1,\cdots,n)

($ Y_i $ ist eine ungeklärte Variable, $ \ boldsymbol {x} $ ist ein erklärender Variablenvektor von k × 1, $ u_i $ ist ein Durchschnitt von 0 und unabhängig voneinander mit einer Normalverteilung der Varianz $ \ sigma ^ 2 $ Repräsentiert den Fehlerterm gemäß). Zu diesem Zeitpunkt folgt $ y_i $ der Normalverteilung des Mittelwerts $ \ boldsymbol {x} _ {i} ^ {'} \ boldsymbol {\ beta} $ und der Varianz $ \ sigma ^ 2 $, so dass die Wahrscheinlichkeitsfunktion ist

\begin{align}
f（\boldsymbol{y} \ | \ \boldsymbol{\beta},\sigma^2) &= \prod_{i=1}^{n}\frac{1}{(2\pi\sigma^2)^{1/2}}\mathrm{exp}\left\{-\frac{(y_i-\boldsymbol{x}_{i}'\boldsymbol{\beta})^2}{2\sigma^2}\right\} \\
&\propto \left(\frac{1}{\sigma^2}\right)^{n/2}\mathrm{exp}\left\{-\frac{\sum_{i=1}^{n}(y_i-\boldsymbol{x}_{i}'\boldsymbol{\beta})^2}{2\sigma^2}\right\} 
\end{align}

Wird sein. In Bezug auf die vorherige Verteilung der Parameter,

\boldsymbol{\beta} \sim N(\boldsymbol{\beta}_0,\boldsymbol{B}_0), \ \ \sigma^2 \sim IG\left(\frac{n_0}{2},\frac{s_0}{2}\right)

(IG bedeutet jedoch eine inverse Gammaverteilung). Unter Anwendung des Bayes-Theorems auf die obige Wahrscheinlichkeit und vorherige Verteilung ist daher die hintere Verteilung

\pi(\boldsymbol{\beta},\sigma^2 \ | \ \boldsymbol{y}) \propto \left(\frac{1}{\sigma^2}\right)^{n/2}\mathrm{exp}\left\{-\frac{\sum_{i=1}^{n}(y_i-\boldsymbol{x}_{i}'\boldsymbol{\beta})^2}{2\sigma^2}\right\}\mathrm{exp}\left\{-\frac{1}{2}(\boldsymbol{\beta}-\boldsymbol{\beta}_0)'\boldsymbol{B}_{0}^{-1}(\boldsymbol{\beta}-\boldsymbol{\beta}_0)\right\}\left(\frac{1}{\sigma^2}\right)^{n_0/2+1}\mathrm{exp}\left(-\frac{s_0}{2\sigma^2}\right)

Wird sein.

Übrigens, obwohl der Kern der posterioren Verteilung erfolgreich abgeleitet werden kann, kann die standardisierte Konstante der posterioren Verteilung nicht berechnet werden, da die integrale Berechnung der obigen Gleichung nicht analytisch durchgeführt werden kann (weshalb eine Probenahme unter Verwendung von MCMC erforderlich ist). Wird sein). Daher wird die Gibbs-Abtastung verwendet, um erfolgreich Zufallszahlen zu erzeugen, die der posterioren Verteilung folgen und die Eigenschaften der gewünschten Verteilung analysieren. Um jedoch eine Gibbs-Abtastung durchzuführen, muss eine vollständige bedingte Verteilung jedes Parameters abgeleitet werden. Dies kann aus dem zuvor abgeleiteten posterioren Verteilungskern berechnet werden, aber der Prozess ist kompliziert. Lassen Sie ihn daher weg und zeigen Sie nur das Ergebnis an (für Details jede Person ["Bayes Computational Statistics" (Asakura Shoten)) (https: // www). .amazon.co.jp /% E3% 83% 99% E3% 82% A4% E3% 82% BA% E8% A8% 88% E7% AE% 97% E7% B5% B1% E8% A8% 88% E5% AD% A6-% E7% B5% B1% E8% A8% 88% E8% A7% A3% E6% 9E% 90% E3% 82% B9% E3% 82% BF% E3% 83% B3% E3 % 83% 80% E3% 83% BC% E3% 83% 89-% E5% 8F% A4% E6% BE% 84-% E8% 8B% B1% E7% 94% B7 / dp / 4254128568 / ref = sr_1_1 dh = UTF8 & qid = 1485493617 & sr = 8-1 & Schlüsselwörter =% E3% 83% 99% E3% 82% A4% E3% 82% BA% E8% A8% 88% E7% AE% 97% E7% B5% B1% E8% Bitte beziehen Sie sich auf A8% 88) usw.). Die vollständige bedingte Verteilung jedes Parameters

\begin{align}
\pi(\boldsymbol{\beta} \ | \ \sigma^2,\boldsymbol{y}) &= N(\hat{\beta},\hat{\boldsymbol{B}}) \\
\pi(\sigma^2 \ | \ \boldsymbol{\beta},\boldsymbol{y}) &= IG\left(\frac{n+n_0}{2},\frac{\sum_{i=1}^{n}(y_i-\boldsymbol{x}_{i}'\boldsymbol{\beta})^2)+s_0}{2}\right)
\end{align} \\

Jedoch,

\hat{\boldsymbol{B}}^{-1}=\sum_{i=1}^{n}\frac{\boldsymbol{x}_i\boldsymbol{x}_i'}{\sigma^2}+\hat{\boldsymbol{B}}_{0}^{-1}, \ \ \hat{\boldsymbol{\beta}}=\hat{\boldsymbol{B}}\left(\sum_{i=1}^{n}\frac{\boldsymbol{x}_iy_i}{\sigma^2}+\hat{\boldsymbol{B}}_{0}^{-1}\boldsymbol{\beta}_0\right)

Ist. Ausgehend von dieser perfekten bedingten Verteilung kann die Abtastung abwechselnd nacheinander durchgeführt werden.

Implementierung durch Python

Pseudodaten werden erzeugt, eine Bayes'sche Schätzung wird aus den Daten durchgeführt und es wird überprüft, ob sich die Verteilung der Ergebnisse den ursprünglichen Parametern nähert. Sehr einfach, aber der Gibbs-Sampling-Python-Code ist unten.

import numpy as np
from numpy.random import normal, multivariate_normal
from scipy.stats import invgamma
import matplotlib.pyplot as plt

##Pseudodatengenerierung
alpha = 2.5
beta = 0.8
sigma = 10
data = 100

x = np.c_[np.ones(data), rand(data) * data]

y = 2.5 * x[:,0] + 0.8 * x[:,1] + normal(0, sigma, data)

plt.plot(x[:,1],y,'o')



##Gibbs Sampling

#Geben Sie die Anzahl der Einbrennvorgänge und die Anzahl der Proben an
burnin = 10000
it = 2000

beta = np.zeros([burnin + it, 2])
sigma = np.zeros(burnin + it)
B_hat = np.zeros([2,2])
beta_hat = np.zeros(2)
sum1 = np.zeros([2,2])
sum2 = np.zeros(2) 

#Angeben von Hyperparametern
beta0 = np.array([0,0])
B0 = np.array([[1, 0],[0,1]])
n0 = 2
s0 = 2

#Anfangswerteinstellung
beta[0,0] = 1
beta[0,1] = 1
sigma[0] = 1 


for i in range(data):
    sum1[0,0] += x[i,0] ** 2
    sum1[0,1] += x[i,0] * x[i,1]
    sum1[0,1] += x[i,0] * x[i,1]
    sum1[1,1] += x[i,1] ** 2 
    sum2[0] += x[i,0] * y[i]
    sum2[1] += x[i,1] * y[i]

for i in range(burnin + it - 1):
    B_hat = np.linalg.inv(sum1/sigma[i] + np.linalg.inv(B0))
    beta_hat = B_hat @ (sum2/sigma[i] + np.linalg.inv(B0) @ beta0)
    beta[i+1,:] =  multivariate_normal(beta_hat, B_hat)
    sigma[i+1] = invgamma.rvs((data+n0)/2,(np.dot(y-beta[i,0]*x[:,0]-beta[i,1]*x[:,1], y-beta[i,0]*x[:,0]-beta[i,1]*x[:,1])))

Als Ergebnis der Probenahme ist das Beta-Histogramm wie folgt.

Es ist ersichtlich, dass sich die MAP-Schätzung dem anfänglich eingestellten Wert von 0,8 nähert. Ich habs gemacht.

(Details werden später hinzugefügt)

Ich habe versucht, die Bayes'sche lineare Regression durch Gibbs-Sampling in Python zu implementieren

Einführung

Theorie

Implementierung durch Python