introduction

Cet article porte sur l'estimation bayésienne d'un modèle de régression linéaire simple. Plus précisément, la méthode de Monte Carlo par chaîne de Markov (communément appelée MCMC) a été implémentée en python à l'aide de l'échantillonnage de Gibbs. L'exécution de MCMC en python est célèbre pour le module appelé pymc, mais cette fois j'ai osé essayer de l'implémenter uniquement avec numpy et scipy sans utiliser pymc. En écrivant cet article, M. 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 & mots-clés =% E3% 83% 99% E3% 82% A4 % E3% 82% BA% E8% A8% 88% E7% AE% 97% E7% B5% B1% E8% A8% 88) a été largement mentionné.

théorie

Revenons à la théorie, même si elle est simple, avant la mise en œuvre. Le modèle de régression linéaire que vous souhaitez estimer

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

($ Y_i $ est une variable inexpliquée, $ \ boldsymbol {x} $ est un vecteur de variable explicative de k × 1, $ u_i $ est une moyenne de 0, et est indépendant les uns des autres dans une distribution normale avec une variance de $ \ sigma ^ 2 $. Représente le terme d'erreur selon). À ce stade, $ y_i $ suit la distribution normale de la moyenne $ \ boldsymbol {x} _ {i} ^ {'} \ boldsymbol {\ beta} $ et de la variance $ \ sigma ^ 2 $, de sorte que la fonction de vraisemblance est

\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}

Sera. Concernant la distribution a priori des paramètres,

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

(Cependant, IG signifie une distribution gamma inverse). Par conséquent, en appliquant le théorème de Bayes à la probabilité ci-dessus et à la distribution antérieure, la distribution postérieure devient

\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)

Sera.

À propos, bien que le noyau de la distribution postérieure puisse être dérivé avec succès, la constante normalisée de la distribution postérieure ne peut pas être calculée car le calcul intégral de l'équation ci-dessus ne peut pas être effectué analytiquement (c'est pourquoi l'échantillonnage à l'aide de MCMC est nécessaire). Sera). Par conséquent, nous réussissons à générer des nombres aléatoires qui suivent la distribution postérieure par échantillonnage de Gibbs et à analyser les propriétés de la distribution souhaitée, mais afin de réaliser un échantillonnage de Gibbs, nous devons dériver une distribution conditionnelle complète de chaque paramètre. Cela peut être calculé à partir du noyau de distribution postérieur dérivé plus tôt, mais le processus est compliqué, alors omettez-le et affichez uniquement le résultat (pour plus de détails, chaque personne ["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 ? ie = UTF8 & qid = 1485493617 & sr = 8-1 & mots-clés =% E3% 83% 99% E3% 82% A4% E3% 82% BA% E8% A8% 88% E7% AE% 97% E7% B5% B1% E8% Veuillez vous référer à A8% 88) etc.). La distribution conditionnelle complète de chaque paramètre

\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} \\

cependant,

\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)

Est. A partir de cette distribution conditionnelle parfaite, l'échantillonnage peut être effectué alternativement les uns après les autres.

Implémentation par python

Des pseudo données sont générées, une estimation bayésienne est effectuée à partir des données et il est vérifié si la distribution des résultats se rapproche des paramètres d'origine. Très simple, mais le code pytho d'échantillonnage Gibbs est ci-dessous.

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

##Génération de pseudo données
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')



##Échantillonnage Gibbs

#Spécifiez le nombre de brûlures et le nombre d'échantillons
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) 

#Spécification d'hyper paramètres
beta0 = np.array([0,0])
B0 = np.array([[1, 0],[0,1]])
n0 = 2
s0 = 2

#Réglage de la valeur initiale
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])))

À la suite de l'échantillonnage, l'histogramme bêta est le suivant.

On peut voir que l'estimation MAP s'approche de la valeur initialement fixée de 0,8. Je l'ai fait.

(Les détails seront ajoutés plus tard)

J'ai essayé d'implémenter la régression linéaire bayésienne par échantillonnage de Gibbs en python

introduction

théorie

Implémentation par python