개시하다

본고는 간단한 선형 회귀 모델의 베스 추정과 관련된다.구체적으로 말하면 깁스 샘플을 사용하여python에서 마르코프 체인 몬테카로법(통칭 MC)을 실시했다.pymc 이 모듈에서 MCMC는pymc에서 실행하는 것으로 유명하지만 이번에는pymc를 사용하지 않고 numby와scipy의 실현만 시도했다.이 기사를 쓸 때 나는 고칭의《Bays계산통계학》(조창서점)를 참고했다.

이론

실시하기 전에 간단하지만 이론을 되돌아봐야 한다.추정할 선형 회귀 모델은

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

표시($y i$비설명 변수, $\boldsymbol{x}달러는 k×1의 설명 변수 벡터, $ui$는 평균 0, 분산 $\sigma^2달러의 정적 분포에 따라 서로 독립된 오차항을 나타낸다.현재 $y평균 $\boldsymbol{x}{i]^{'}\boldsymbol{\beta}달러, 분산 $\sigma^2달러의 정적 분포에 따라 유사함수

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

되다
또한 매개 변수의 사전 분포에 관하여

\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는 역감마 분포를 나타냅니다.따라서 베이스의 정리를 상술한 유사도와 예분포에 응용하면 후방향 분포

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

되다
그러면 후방 분포된 내핵은 순조롭게 도출할 수 있으나 상단식의 적분 계산은 해석적으로 진행할 수 없기 때문에 후방 분포의 기준화 상수를 계산할 수 없다(이렇기 때문에 MCC로 샘플링을 해야 한다).따라서 지부스 샘플링을 통해 성공적으로 생성된 후 분포된 무작위 수를 분석하고 기대 분포의 성질을 분석하지만 지부스 샘플링을 하기 위해서는 각 파라미터의 완전한 조건 분포를 도출해야 한다.이것은 방금 도출된 백업 분포의 내핵에서 계산할 수 있으나 과정이 복잡하기 때문에 생략한 후에 결과만 나타낼 수 있다(자세한 내용은 각자《Bays계산통계학》(조창서점) 등을 참고하시오).매개 변수의 전체 조건 분포

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

그렇지만

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

.이 완전한 조건 분포에서 순서대로 교대로 샘플을 채취할 수 있다.

python 기반 설치

위조 데이터를 생성하고 이 데이터에 따라 바이스 평가를 하여 결과의 분포가 원시 파라미터에 가깝는지 검증한다.매우 간단하지만, 석고 샘플링의pytho 코드는 다음과 같다.

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

## 擬似データの生成
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')



## ギブス・サンプリング

# burn inの回数とサンプリング数を指定
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) 

# ハイパーパラメーターの指定
beta0 = np.array([0,0])
B0 = np.array([[1, 0],[0,1]])
n0 = 2
s0 = 2

# 初期値の設定
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])))

샘플링의 결과 베타의 직사각형은 다음과 같다.

MAP는 애초 설정한 값인 0.8에 근접한 것으로 추정된다.잘 됐다.
(세부 사항은 잠시 후에 추기)