Pyro 0.3.0 : Examples : ベイジアン回帰 – イントロダクション (Part 1) (翻訳)
翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 12/18/2018 (v0.3.0)
* 本ページは、Pyro のドキュメント Examples : Bayesian Regression – Introduction (Part 1) を
翻訳した上で適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
ベイジアン回帰 – イントロダクション (Part 1)
回帰は機械学習における最も一般的で基本的な教師あり学習タスクの一つです。次の形式のデータセット $\mathcal{D}$ が与えられたと仮定します :
\[
\mathcal{D} = \{ (X_i, y_i) \} \qquad \text{for}\qquad i=1,2,…,N
\]
線形回帰のゴールは関数を次の形式のデータに fit させることです :
\[
y = w X + b + \epsilon
\]
ここで $w$ と $b$ は学習可能なパラメータで $\epsilon$ は観測ノイズを表します。特に $w$ は重み行列で $b$ はバイアス・ベクトルです。
最初に PyTorch で線形回帰を実装してパラメータ $w$ と $b$ のための点推定を学習しましょう。それからベイジアン線形回帰を実装するために Pyro を使用することによりどのように不確かさを推定に組み込むかを見ましょう。
セットアップ
必要なモジュールをインポートすることから始めましょう :
import os
from functools import partial
import numpy as np
import pandas as pd
import seaborn as sns
import torch
import torch.nn as nn
import matplotlib.pyplot as plt
import pyro
from pyro.distributions import Normal, Uniform, Delta
from pyro.infer import SVI, Trace_ELBO
from pyro.optim import Adam
from pyro.distributions.util import logsumexp
from pyro.infer import EmpiricalMarginal, SVI, Trace_ELBO, TracePredictive
from pyro.infer.mcmc import MCMC, NUTS
import pyro.optim as optim
import pyro.poutine as poutine
# for CI testing
smoke_test = ('CI' in os.environ)
assert pyro.__version__.startswith('0.3.0')
pyro.enable_validation(True)
pyro.set_rng_seed(1)
pyro.enable_validation(True)
データセット
次のサンプルは [1] から改作しました。Terrain Ruggedness Index (直訳: 地形起伏指標) (データセットで起伏の多い変数) により計測された国の topographic heterogeneity (直訳: 位置特異的な不均質) と一人あたりの GDP の関係性を探究したいです。特に、地形起伏 (= terrain ruggedness) あるいは悪い地形 (= bad geography) はアフリカの外ではより貧しい経済的なパフォーマンスに関係しますが、起伏の激しい地形はアフリカ諸国のための収入上は逆の効果を持っていたことが [1] で著者により記されています。データ [2] を見てこの関係性を調査しましょう。データセットから 3 つの特徴に注目します :
- rugged: Terrain Ruggedness Index を定量化する
- cont_africa: 与えられた国がアフリカにあるか否か
- rgdppc_2000: 2000 年についての一人あたりの Real GDP
応答変数 GDP のために対数を取ります、何故ならばそれは指数関数的に変化する傾向があるからです。
DATA_URL = "https://d2fefpcigoriu7.cloudfront.net/datasets/rugged_data.csv" data = pd.read_csv(DATA_URL, encoding="ISO-8859-1") df = data[["cont_africa", "rugged", "rgdppc_2000"]] df = df[np.isfinite(df.rgdppc_2000)] df["rgdppc_2000"] = np.log(df["rgdppc_2000"])
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True)
african_nations = data[data["cont_africa"] == 1]
non_african_nations = data[data["cont_africa"] == 0]
sns.scatterplot(non_african_nations["rugged"],
np.log(non_african_nations["rgdppc_2000"]),
ax=ax[0])
ax[0].set(xlabel="Terrain Ruggedness Index",
ylabel="log GDP (2000)",
title="Non African Nations")
sns.scatterplot(african_nations["rugged"],
np.log(african_nations["rgdppc_2000"]),
ax=ax[1])
ax[1].set(xlabel="Terrain Ruggedness Index",
ylabel="log GDP (2000)",
title="African Nations")
[Text(0, 0.5, 'log GDP (2000)'), Text(0.5, 0, 'Terrain Ruggedness Index'), Text(0.5, 1.0, 'African Nations')]
線形回帰
国の一人あたりの log GDP をデータセットからの 2 つの特徴の関数として予測したいです – 国がアフリカ内であるか、そしてその Terrain Ruggedness Index です。私達の回帰モデルを定義しましょう。このために PyTorch の nn.Module を使用します。入力 \(X\) は size \(N \times 2\) の行列で出力 \(y\) は size \(2 \times 1\) のベクトルです。関数 nn.Linear(p,1) は形式 \(Xw + b\) の線形変換を定義します、そこでは $w$ は重み行列で $b$ は additive バイアスです。起伏と国がアフリカ内か否かの間の相関性を捕捉することを意図した特別な self.factor 項を含みます。
forward() メソッドで非線形性を追加することでこれを容易にロジスティック回帰にできることに注意してください。
class RegressionModel(nn.Module):
def __init__(self, p):
# p = number of features
super(RegressionModel, self).__init__()
self.linear = nn.Linear(p, 1)
self.factor = nn.Parameter(torch.tensor(1.))
def forward(self, x):
return self.linear(x) + (self.factor * x[:, 0] * x[:, 1]).unsqueeze(1)
p = 2 # number of features
regression_model = RegressionModel(p)
訓練
損失として MSE (mean squared error) をそして optimizer として Adam を使用します。私達は上の regression_model ニューラルネットのパラメータを最適化したいです。0.05 の幾分大きな学習率を使用して 500 反復を実行します :
loss_fn = torch.nn.MSELoss(reduction='sum')
optim = torch.optim.Adam(regression_model.parameters(), lr=0.05)
num_iterations = 1000 if not smoke_test else 2
data = torch.tensor(df.values, dtype=torch.float)
x_data, y_data = data[:, :-1], data[:, -1]
def main():
x_data = data[:, :-1]
y_data = data[:, -1]
for j in range(num_iterations):
# run the model forward on the data
y_pred = regression_model(x_data).squeeze(-1)
# calculate the mse loss
loss = loss_fn(y_pred, y_data)
# initialize gradients to zero
optim.zero_grad()
# backpropagate
loss.backward()
# take a gradient step
optim.step()
if (j + 1) % 50 == 0:
print("[iteration %04d] loss: %.4f" % (j + 1, loss.item()))
# Inspect learned parameters
print("Learned parameters:")
for name, param in regression_model.named_parameters():
print(name, param.data.numpy())
main()
[iteration 0050] loss: 3723.7849 [iteration 0100] loss: 1777.9851 [iteration 0150] loss: 1226.7543 [iteration 0200] loss: 927.1514 [iteration 0250] loss: 706.6473 [iteration 0300] loss: 536.0911 [iteration 0350] loss: 408.0940 [iteration 0400] loss: 316.1810 [iteration 0450] loss: 252.9711 [iteration 0500] loss: 211.2545 [iteration 0550] loss: 184.7950 [iteration 0600] loss: 168.6502 [iteration 0650] loss: 159.1673 [iteration 0700] loss: 153.8036 [iteration 0750] loss: 150.8815 [iteration 0800] loss: 149.3482 [iteration 0850] loss: 148.5732 [iteration 0900] loss: 148.1960 [iteration 0950] loss: 148.0193 [iteration 1000] loss: 147.9397 Learned parameters: factor 0.37248382 linear.weight [[-1.90511 -0.18619268]] linear.bias [9.188872]
ベイジアン・モデリング はモデルの不確かさについて推論するためのシステマティックなフレームワークを提供します。単に点推定を学習する代わりに、観測データと調和した変数に渡る分布を学習していきます。
ベイジアン回帰
私達の線形回帰をベイジアンにするためには、パラメータ $w$ と $b$ の上に事前分布を置く必要があります。これらは (どのようなデータも観測する前に) $w$ と $b$ の合理的な値についての私達の事前信念を表わす分布です。
random_module()
これを行なうために、既存のモデルのパラメータを確率変数に ‘lift’ します。Pyro では random_module() を通してこれを行なうことができます、これは与えられた nn.Module を取りそしてそれを同じモジュールに渡る分布に効果的に変更します ; 私達のケースではこれは regressor に渡る分布です。具体的には、元の回帰モデルの各パラメータは提供された事前分布からサンプリングされます。これは vanilla 回帰モデルをベイジアン設定での利用のために再目的化することを可能にします。例えば :
loc = torch.zeros(1, 1)
scale = torch.ones(1, 1)
# define a unit normal prior
prior = Normal(loc, scale)
# overload the parameters in the regression module with samples from the prior
lifted_module = pyro.random_module("regression_module", nn, prior)
# sample a nn from the prior
sampled_reg_model = lifted_module()
モデル
今ではモデルを指定するために必要な総ての構成要素を持っています。最初に重みとバイアス、そして factor に渡る事前分布を定義します。モデルの異なる潜在変数のために使用している事前分布に注意してください。intercept パラメータ上の事前分布は (これにデータから学習されることを望む程度に) 非常に flat です。データへの overfitting を回避するためにregression coefficients 上で weakly regularizing prior を使用しています。
regression_model を random_module でラップして regressor のインスタンス、lifted_reg_model からサンプリングします。それから x_data 上で regressor を実行します。最後に学習された観測ノイズ sigma を伴う観測データ y_data 上で条件付けるために pyro.sample ステートメントへの obs 引数を使用します。
def model(x_data, y_data):
# weight and bias priors
w_prior = Normal(torch.zeros(1, 2), torch.ones(1, 2)).to_event(1)
b_prior = Normal(torch.tensor([[8.]]), torch.tensor([[1000.]])).to_event(1)
f_prior = Normal(0., 1.)
priors = {'linear.weight': w_prior, 'linear.bias': b_prior, 'factor': f_prior}
scale = pyro.sample("sigma", Uniform(0., 10.))
# lift module parameters to random variables sampled from the priors
lifted_module = pyro.random_module("module", regression_model, priors)
# sample a nn (which also samples w and b)
lifted_reg_model = lifted_module()
with pyro.plate("map", len(x_data)):
# run the nn forward on data
prediction_mean = lifted_reg_model(x_data).squeeze(-1)
# condition on the observed data
pyro.sample("obs",
Normal(prediction_mean, scale),
obs=y_data)
return prediction_mean
ガイド
推論を行なうためにはガイド, i.e. 分布の変分ファミリが必要となります。モデルの総ての分布上に対角共分散を持つ正規分布を自動的に配置するために Pyro の autoguide ライブラリ を使用します。内部的には、これはモデルの各 sample() に対応する学習パラメータを持つ Normal 分布を持つガイド関数を定義します。Autoguide はまた AutoDelta で MAP 推定を学習することや AutoGuideList でガイドを組み立てることもサポートします (より多くの情報については docs 参照)。Part II ではどのようにガイドを主導で書くかを探究します。
from pyro.contrib.autoguide import AutoDiagonalNormal guide = AutoDiagonalNormal(model)
推論
推論を行なうためには確率的変分推論 (SVI) を使用します (SVI のイントロダクションのためには、SVI パート I 参照)。ちょうど非ベイジアン線形回帰におけるように、訓練ループの各反復は勾配ステップを取ります、この場合、MSE 損失の代わりに (SVI に渡す Trace_ELBO オブジェクトを構築することにより) ELBO 目的関数を使用するという違いを伴います。
optim = Adam({"lr": 0.03})
svi = SVI(model, guide, optim, loss=Trace_ELBO(), num_samples=1000)
ここで Adam は torch.optim.Adam 回りの薄いラッパーです (議論のためには こちら を参照)。ELBO 勾配ステップを取るためには SVI の step メソッドを単純に呼び出します。step に渡す data 引数は model() と guide() の両者に渡されることに注意してください。完全な訓練ループは次のようなものです :
def train():
pyro.clear_param_store()
for j in range(num_iterations):
# calculate the loss and take a gradient step
loss = svi.step(x_data, y_data)
if j % 100 == 0:
print("[iteration %04d] loss: %.4f" % (j + 1, loss / len(data)))
train()
[iteration 0001] loss: 18.3907 [iteration 0101] loss: 3.1146 [iteration 0201] loss: 3.1067 [iteration 0301] loss: 2.8602 [iteration 0401] loss: 2.7770 [iteration 0501] loss: 2.6181 [iteration 0601] loss: 2.4298 [iteration 0701] loss: 2.0160 [iteration 0801] loss: 1.7814 [iteration 0901] loss: 1.5811
for name, value in pyro.get_param_store().items():
print(name, pyro.param(name))
auto_loc tensor([-2.2026, 0.2936, -1.8873, -0.1607, 9.1753], requires_grad=True) auto_scale tensor([0.2285, 0.0954, 0.1376, 0.0600, 0.1042], grad_fn=)
見れるように、単なる点推定の代わりに、今では学習されたパラメータのための不確かさ推定 (auto_scale) を持ちます。Autoguide は潜在変数を tensor にパックすることに注意してください、このケースでは、モデルでサンプリングされた変数毎に一つのエントリです。
モデル評価
モデルを評価するために、幾つかの予測サンプルを生成して事後分布を見ます。私達の変分分布は完全にパラメータ化されていますので、lift されたモデルを単に foward 実行できます。Pyro で値を登録するために Delta 分布でモデルをラップします。それから svi.run() で posterior 分布内に実行トレースをストアします。
get_marginal = lambda traces, sites:EmpiricalMarginal(traces, sites)._get_samples_and_weights()[0].detach().cpu().numpy()
def summary(traces, sites):
marginal = get_marginal(traces, sites)
site_stats = {}
for i in range(marginal.shape[1]):
site_name = sites[i]
marginal_site = pd.DataFrame(marginal[:, i]).transpose()
describe = partial(pd.Series.describe, percentiles=[.05, 0.25, 0.5, 0.75, 0.95])
site_stats[site_name] = marginal_site.apply(describe, axis=1) \
[["mean", "std", "5%", "25%", "50%", "75%", "95%"]]
return site_stats
def wrapped_model(x_data, y_data):
pyro.sample("prediction", Delta(model(x_data, y_data)))
posterior = svi.run(x_data, y_data)
# posterior predictive distribution we can get samples from
trace_pred = TracePredictive(wrapped_model,
posterior,
num_samples=1000)
post_pred = trace_pred.run(x_data, None)
post_summary = summary(post_pred, sites= ['prediction', 'obs'])
mu = post_summary["prediction"]
y = post_summary["obs"]
predictions = pd.DataFrame({
"cont_africa": x_data[:, 0],
"rugged": x_data[:, 1],
"mu_mean": mu["mean"],
"mu_perc_5": mu["5%"],
"mu_perc_95": mu["95%"],
"y_mean": y["mean"],
"y_perc_5": y["5%"],
"y_perc_95": y["95%"],
"true_gdp": y_data,
})
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True)
african_nations = predictions[predictions["cont_africa"] == 1]
non_african_nations = predictions[predictions["cont_africa"] == 0]
african_nations = african_nations.sort_values(by=["rugged"])
non_african_nations = non_african_nations.sort_values(by=["rugged"])
fig.suptitle("Regression line 90% CI", fontsize=16)
ax[0].plot(non_african_nations["rugged"],
non_african_nations["mu_mean"])
ax[0].fill_between(non_african_nations["rugged"],
non_african_nations["mu_perc_5"],
non_african_nations["mu_perc_95"],
alpha=0.5)
ax[0].plot(non_african_nations["rugged"],
non_african_nations["true_gdp"],
"o")
ax[0].set(xlabel="Terrain Ruggedness Index",
ylabel="log GDP (2000)",
title="Non African Nations")
idx = np.argsort(african_nations["rugged"])
ax[1].plot(african_nations["rugged"],
african_nations["mu_mean"])
ax[1].fill_between(african_nations["rugged"],
african_nations["mu_perc_5"],
african_nations["mu_perc_95"],
alpha=0.5)
ax[1].plot(african_nations["rugged"],
african_nations["true_gdp"],
"o")
ax[1].set(xlabel="Terrain Ruggedness Index",
ylabel="log GDP (2000)",
title="African Nations")
[Text(0, 0.5, 'log GDP (2000)'), Text(0.5, 0, 'Terrain Ruggedness Index'), Text(0.5, 1.0, 'African Nations')]
上の図は回帰ラインの私達の推定において不確かさを示します。起伏のより低い値のために数多くのより多いデータポイントがあることに注意してください、そのようなものとして、最善の fit のラインのためにより少ない小刻みに動く余地があります。これは平均回りの 90% CI で反映されます。データポイントの殆どは実際には 90% CI の外側に在ることも見ることができて、これは想定されたものです、何故ならば sigma により影響される outcome 変数をプロットしていないからです！次にそのようにしましょう。
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True)
fig.suptitle("Posterior predictive distribution with 90% CI", fontsize=16)
ax[0].plot(non_african_nations["rugged"],
non_african_nations["y_mean"])
ax[0].fill_between(non_african_nations["rugged"],
non_african_nations["y_perc_5"],
non_african_nations["y_perc_95"],
alpha=0.5)
ax[0].plot(non_african_nations["rugged"],
non_african_nations["true_gdp"],
"o")
ax[0].set(xlabel="Terrain Ruggedness Index",
ylabel="log GDP (2000)",
title="Non African Nations")
idx = np.argsort(african_nations["rugged"])
ax[1].plot(african_nations["rugged"],
african_nations["y_mean"])
ax[1].fill_between(african_nations["rugged"],
african_nations["y_perc_5"],
african_nations["y_perc_95"],
alpha=0.5)
ax[1].plot(african_nations["rugged"],
african_nations["true_gdp"],
"o")
ax[1].set(xlabel="Terrain Ruggedness Index",
ylabel="log GDP (2000)",
title="African Nations")
[Text(0, 0.5, 'log GDP (2000)'), Text(0.5, 0, 'Terrain Ruggedness Index'), Text(0.5, 1.0, 'African Nations')]
We observe that the outcome from our model and the 90% CI accounts for the majority of the data points that we observe in practice. モデルが妥当な予測を与えるかを見るためにそのような事後予測チェックを行なうことは通常は良いアイデアです。
# we need to prepend `module$$$` to all parameters of nn.Modules since
# that is how they are stored in the ParamStore
weight = get_marginal(posterior, ['module$$$linear.weight']).squeeze(1).squeeze(1)
factor = get_marginal(posterior, ['module$$$factor'])
gamma_within_africa = weight[:, 1] + factor.squeeze(1)
gamma_outside_africa = weight[:, 1]
fig = plt.figure(figsize=(10, 6))
sns.distplot(gamma_within_africa, kde_kws={"label": "African nations"},)
sns.distplot(gamma_outside_africa, kde_kws={"label": "Non-African nations"})
fig.suptitle("Density of Slope : log(GDP) vs. Terrain Ruggedness", fontsize=16)
Text(0.5, 0.98, 'Density of Slope : log(GDP) vs. Terrain Ruggedness')
次のセクションでは、変分推論のためにどのようにガイドを書くかを見てそして HMC を通した推論による結果を比較します。
Github 上の toy データセットを持つサンプルを見てください。
References
- McElreath, D., Statistical Rethinking, Chapter 7, 2016
- Nunn, N. & Puga, D., Ruggedness: The blessing of bad geography in Africa”, Review of Economics and Statistics 94(1), Feb. 2012
以上