ベイズ線形回帰

Contents

ベイズ線形回帰#

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

# データを作成
n = 1000

from scipy.stats import multivariate_normal
mean = np.array([3, 5])
Sigma = np.array([
    [1, 0.5],
    [0.5, 2],
])
X = multivariate_normal.rvs(mean=mean, cov=Sigma, size=n, random_state=0)

import statsmodels.api as sm
X = sm.add_constant(X)

# 真のパラメータ
beta = np.array([2, 3, 4])

データが均一分散の場合#

# 均一分散の場合
e = np.random.normal(loc=0, scale=1, size=n)
y = X @ beta + e

# 頻度主義
import statsmodels.api as sm
ols = sm.OLS(y, X).fit(cov_type="HC1")
ols.summary()

OLS Regression Results
Dep. Variable:	y	R-squared:	0.980
Model:	OLS	Adj. R-squared:	0.980
Method:	Least Squares	F-statistic:	2.516e+04
Date:	Fri, 04 Apr 2025	Prob (F-statistic):	0.00
Time:	11:03:13	Log-Likelihood:	-1418.4
No. Observations:	1000	AIC:	2843.
Df Residuals:	997	BIC:	2858.
Df Model:	2
Covariance Type:	HC1

	coef	std err	z	P>\|z\|	[0.025	0.975]
const	1.9911	0.133	15.018	0.000	1.731	2.251
x1	3.0074	0.035	85.559	0.000	2.939	3.076
x2	3.9994	0.025	158.700	0.000	3.950	4.049

Omnibus:	1.626	Durbin-Watson:	1.954
Prob(Omnibus):	0.443	Jarque-Bera (JB):	1.596
Skew:	0.039	Prob(JB):	0.450
Kurtosis:	2.821	Cond. No.	25.9

Notes:
[1] Standard Errors are heteroscedasticity robust (HC1)

import pymc as pm
import arviz as az

model = pm.Model()
with model:
    beta0 = pm.Normal("beta0", mu=0, sigma=1)
    beta1 = pm.Normal("beta1", mu=0, sigma=1)
    beta2 = pm.Normal("beta2", mu=0, sigma=1)
    sigma = pm.HalfNormal("sigma", sigma=1)  # 分散なので非負の分布を使う

    # 平均値 mu
    mu = beta0 + beta1 * X[:, 1] + beta2 * X[:, 2]
    # 観測値をもつ確率変数は_obsとする
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

# モデルをGraphvizで表示
pm.model_to_graphviz(model)

../../_images/8db035c7e654249210826758756c5b9f3fabd85298d5aa28f5526469e1117775.svg

# ベイズ線形回帰モデルをサンプリング
with model:
    idata = pm.sample(
        chains=2,
        tune=1000, # バーンイン期間の、捨てるサンプル数
        draws=2000, # 採用するサンプル数
        random_seed=0,
    )

# 各chainsの結果を表示
az.plot_trace(idata, figsize=[4, 4])
plt.tight_layout()
plt.show()

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [beta0, beta1, beta2, sigma]

                                                                                                                   
  Progress                    Draws   Divergences   Step size   Grad evals   Sampling Speed   Elapsed   Remaining  
 ───────────────────────────────────────────────────────────────────────────────────────────────────────────────── 
  ━━━━━━━━━━━━━━━━━━━━━━━━━   3000    0             0.18        15           662.70 draws/s   0:00:04   0:00:00    
  ━━━━━━━━━━━━━━━━━━━━━━━━━   3000    0             0.18        15           683.10 draws/s   0:00:04   0:00:00

Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 5 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

../../_images/8852a233c812c3f224adc7d85f0ca4b2f4ae230316ab268c6ae897e1c426c39b.png

az.plot_posterior(idata)
plt.show()

../../_images/0dab11c8c6316d7a1a64b8446edb08abd91b9f5aefdc69b9c52a226837ee5eb9.png

データが不均一分散の場合#

# 不均一分散の場合
def normalize(x):
    return (x - x.min()) / (x.max() - x.min())

sigma = 1 + normalize(X[:, 1] + X[:, 2]) * 3
e = np.random.normal(loc=0, scale=sigma, size=n)
y = X @ beta + e

頻度主義 & 不均一分散に頑健な誤差推定#

# 頻度主義
import statsmodels.api as sm
ols = sm.OLS(y, X).fit(cov_type="HC1")
ols.summary()

OLS Regression Results
Dep. Variable:	y	R-squared:	0.892
Model:	OLS	Adj. R-squared:	0.892
Method:	Least Squares	F-statistic:	3270.
Date:	Fri, 04 Apr 2025	Prob (F-statistic):	0.00
Time:	11:03:39	Log-Likelihood:	-2312.1
No. Observations:	1000	AIC:	4630.
Df Residuals:	997	BIC:	4645.
Df Model:	2
Covariance Type:	HC1

	coef	std err	z	P>\|z\|	[0.025	0.975]
const	2.1773	0.336	6.481	0.000	1.519	2.836
x1	2.9092	0.085	34.419	0.000	2.744	3.075
x2	4.0083	0.060	66.504	0.000	3.890	4.126

Omnibus:	14.211	Durbin-Watson:	2.074
Prob(Omnibus):	0.001	Jarque-Bera (JB):	16.433
Skew:	0.214	Prob(JB):	0.000270
Kurtosis:	3.460	Cond. No.	25.9

Notes:
[1] Standard Errors are heteroscedasticity robust (HC1)

↑ 切片の推定にバイアスが入っている

均一分散を想定したベイズ線形回帰#

import pymc as pm
import arviz as az

model = pm.Model()
with model:
    beta0 = pm.Normal("beta0", mu=0, sigma=1)
    beta1 = pm.Normal("beta1", mu=0, sigma=1)
    beta2 = pm.Normal("beta2", mu=0, sigma=1)
    sigma = pm.HalfNormal("sigma", sigma=1)  # 分散なので非負の分布を使う

    # 平均値 mu
    mu = beta0 + beta1 * X[:, 1] + beta2 * X[:, 2]
    # 観測値をもつ確率変数は_obsとする
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

# モデルをGraphvizで表示
pm.model_to_graphviz(model)

../../_images/8db035c7e654249210826758756c5b9f3fabd85298d5aa28f5526469e1117775.svg

# ベイズ線形回帰モデルをサンプリング
with model:
    idata = pm.sample(
        chains=2,
        tune=1000, # バーンイン期間の、捨てるサンプル数
        draws=2000, # 採用するサンプル数
        random_seed=0,
    )

# 各chainsの結果を表示
az.plot_trace(idata, figsize=[4, 4])
plt.tight_layout()
plt.show()

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [beta0, beta1, beta2, sigma]

                                                                                                                   
  Progress                    Draws   Divergences   Step size   Grad evals   Sampling Speed   Elapsed   Remaining  
 ───────────────────────────────────────────────────────────────────────────────────────────────────────────────── 
  ━━━━━━━━━━━━━━━━━━━━━━━━━   3000    0             0.17        15           710.75 draws/s   0:00:04   0:00:00    
  ━━━━━━━━━━━━━━━━━━━━━━━━━   3000    0             0.23        31           709.00 draws/s   0:00:04   0:00:00

Sampling 2 chains for 1_000 tune and 2_000 draw iterations (2_000 + 4_000 draws total) took 4 seconds.

We recommend running at least 4 chains for robust computation of convergence diagnostics

../../_images/f3fd536ed5a8e33a3e671360fa1a1ff85cf4ec456c61fe3bb0febac448e1cd9e.png

az.plot_posterior(idata)
plt.show()

../../_images/817689515c4a3165c43266ef24dc0843997cc940151137236cf04ee1a1df9d39.png

不均一分散を想定したベイズ線形回帰（WIP）#

分散をxの関数にしたかった。以下コードで推定できるが個々の $σ_{i}$ が別々に推定される形になって結果が見づらい。もっといい表し方はないものか。

import pymc as pm
import arviz as az

model = pm.Model()
with model:
    beta0 = pm.Normal("beta0", mu=0, sigma=1)
    beta1 = pm.Normal("beta1", mu=0, sigma=1)
    beta2 = pm.Normal("beta2", mu=0, sigma=1)

    # 誤差分散にも線形モデルを入れる
    w0 = pm.Normal("w0", mu=0, sigma=1)
    w1 = pm.Normal("w1", mu=0, sigma=1)
    w2 = pm.Normal("w2", mu=0, sigma=1)
    lam = pm.math.exp(w0 + w1 * X[:, 1] + w2 * X[:, 2])
    sigma = pm.Exponential("sigma", lam=lam)  # 分散なので非負の分布を使う

    # 平均値 mu
    mu = beta0 + beta1 * X[:, 1] + beta2 * X[:, 2]
    # 観測値をもつ確率変数は_obsとする
    y_obs = pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

# モデルをGraphvizで表示
pm.model_to_graphviz(model)

# ベイズ線形回帰モデルをサンプリング
with model:
    idata = pm.sample(
        chains=2,
        tune=1000, # バーンイン期間の、捨てるサンプル数
        draws=2000, # 採用するサンプル数
        random_seed=0,
    )

# 各chainsの結果を表示
az.plot_trace(idata, figsize=[4, 4])
plt.tight_layout()
plt.show()

az.plot_posterior(idata)
plt.show()