PyTorch : Pyro examples : ベイジアン回帰 (翻訳)
翻訳 : (株)クラスキャット セールスインフォメーション
更新日時 : 11/20/2018 (v0.2.1)
作成日時 : 10/28/2018 (v0.2.1)
* 本ページは、Pyro のドキュメント Examples : Bayesian Regression を翻訳した上で適宜、補足説明したものです：
* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。
ベイジアン回帰
回帰は機械学習における最も一般的で基本的な教師あり学習タスクの一つです。次の形式のデータセット $\mathcal{D}$ が与えられたと仮定します :
\[
\mathcal{D} = \{ (X_i, y_i) \} \qquad \text{for}\qquad i=1,2,…,N
\]
線形回帰の目標は関数を次の形式のデータにフィットさせることです :
\[
y = w X + b + \epsilon
\]
ここで $w$ と $b$ は学習パラメータで $\epsilon$ は観測ノイズを表します。特に $w$ は重み行列で $b$ はバイアス・ベクトルです。
最初に PyTorch で線形回帰を実装してパラメータ $w$ と $b$ のための点推定を学習しましょう。それからベイジアン線形回帰を実装するために Pyro を使用することによりどのように不確かさを推定に組み込むかを見ましょう。
セットアップ
いつものように、必要なモジュールをインポートすることから始めましょう :
import os
import numpy as np
import torch
import torch.nn as nn
import pyro
from pyro.distributions import Normal
from pyro.infer import SVI, Trace_ELBO
from pyro.optim import Adam
# for CI testing
smoke_test = ('CI' in os.environ)
pyro.enable_validation(True)
データ
一つの特徴と w=3 と b=1 を持つ toy データセットを次のように生成します :
N = 100 # size of toy data
def build_linear_dataset(N, p=1, noise_std=0.01):
X = np.random.rand(N, p)
# w = 3
w = 3 * np.ones(p)
# b = 1
y = np.matmul(X, w) + np.repeat(1, N) + np.random.normal(0, noise_std, size=N)
y = y.reshape(N, 1)
X, y = torch.tensor(X).type(torch.Tensor), torch.tensor(y).type(torch.Tensor)
data = torch.cat((X, y), 1)
assert data.shape == (N, p + 1)
return data
固定された観測ノイズ $\sigma = 0.1$ を持つデータを生成することに注意してください。
回帰
さて回帰モデルを定義しましょう。このために PyTorch の nn.Module を使用します。入力 $X$ はサイズ $N \times p$ の行列で出力 $y$ はサイズ $p \times 1$ のベクトルです。関数 nn.Linear(p, 1) は形式 $Xw + b$ の線形変換を定義します、ここで $w$ は重み行列で $b$ は追加的なバイアスです。見て取れるように、forward メソッドで非線形を追加することでこれを容易にロジスティック回帰にすることができます。
class RegressionModel(nn.Module):
def __init__(self, p):
# p = number of features
super(RegressionModel, self).__init__()
self.linear = nn.Linear(p, 1)
def forward(self, x):
return self.linear(x)
regression_model = RegressionModel(1)
訓練
損失として MSE をそして optimizer として Adam を使用します。私達は上の regression_model ニューラルネットのパラメータを最適化したいです。0.01 の幾分大きな学習率を使用して 500 反復を実行します :
loss_fn = torch.nn.MSELoss(size_average=False)
optim = torch.optim.Adam(regression_model.parameters(), lr=0.05)
num_iterations = 1000 if not smoke_test else 2
def main():
data = build_linear_dataset(N)
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("%s: %.3f" % (name, param.data.numpy()))
if __name__ == '__main__':
main()
サンプル出力:
[iteration 0400] loss: 0.0105 [iteration 0450] loss: 0.0096 [iteration 0500] loss: 0.0095 [iteration 0550] loss: 0.0095 [iteration 0600] loss: 0.0095 [iteration 0650] loss: 0.0095 [iteration 0700] loss: 0.0095 [iteration 0750] loss: 0.0095 [iteration 0800] loss: 0.0095 [iteration 0850] loss: 0.0095 [iteration 0900] loss: 0.0095 [iteration 0950] loss: 0.0095 [iteration 1000] loss: 0.0095 Learned parameters: linear.weight: 3.004 linear.bias: 0.997
悪くはありません – regressor は $w = 3, b = 1$ の正解に非常に近いパラメータを学習したことを見て取れます。しかしこれらの点推定にどの程度の信頼をおけるのでしょう？
ベイジアン・モデリング は モデルの不確かさについて推論するためのシステマティックなフレームワークを提供します。単に点推定を学習する代わりに、観測されたデータと調和したパラメータ $w$ と $b$ の値に渡る 分布 を学習していきます。
ベイジアン回帰
私達の線形回帰をベイジアンにするためには、$w$ と $b$ 上に事前分布を置く必要があります。これらは (どのようなデータも観測する前に) $w$ と $b$ のための私達の事前信念を表わす分布です。
random_module()
これを行なうために、パラメータ $w$ と $b$ を確率変数に ‘lift’ します。Pyro では random_module() を通してこれを行なうことができます、これは与えられた nn.Module を取りそしてそれを同じモジュールに渡る分布に効果的に変えます; 私達のケースではこれは regressor (独立変数) に渡る分布です。具体的には、元の回帰モデルの各パラメータは提供された事前分布からサンプリングされます。これは vanilla 回帰モデルをベイジアン設定での利用のために再目的化することを可能にします。例えば :
loc = torch.zeros(1, 1)
scale = torch.ones(1, 1)
# 単位正規事前分布を定義する
prior = Normal(loc, scale)
# 回帰モジュールのパラメータを事前分布からのサンプリングでオーバーロードする
lifted_module = pyro.random_module("regression_module", regression_model, prior)
# 事前分布から regressor をサンプリングする
sampled_reg_model = lifted_module()
モデル
今ではモデルを指定するために必要な総ての構成要素を持っています。最初に w と b に渡る事前分布を定義します。パラメータについては先験的には不確かなので、比較的広い事前分布 N(μ=0,σ=10) を使用します。それから regression_model を random_module でラップして regressor のインスタンス、lifted_reg_model からサンプリングします。それから入力 x_data 上で regressor に forward を実行します。最後に観測されたデータ y_data 上で条件付けるために pyro.sample ステートメントに obs 引数を使用します。データを生成するために使用された同じ固定された観測ノイズを使用することに注意してください。
def model(data):
# パラメータに渡る単位正規事前分布を作成する
loc, scale = torch.zeros(1, 1), 10 * torch.ones(1, 1)
bias_loc, bias_scale = torch.zeros(1), 10 * torch.ones(1)
w_prior = Normal(loc, scale).independent(1)
b_prior = Normal(bias_loc, bias_scale).independent(1)
priors = {'linear.weight': w_prior, 'linear.bias': b_prior}
# lift module parameters to random variables sampled from the priors
lifted_module = pyro.random_module("module", regression_model, priors)
# sample a regressor (which also samples w and b)
lifted_reg_model = lifted_module()
with pyro.iarange("map", N):
x_data = data[:, :-1]
y_data = data[:, -1]
# run the regressor forward conditioned on data
prediction_mean = lifted_reg_model(x_data).squeeze(-1)
# condition on the observed data
pyro.sample("obs",
Normal(prediction_mean, 0.1 * torch.ones(data.size(0))),
obs=y_data)
ガイド
推論を行なうためにはガイド, i.e. $w$ と $b$ に渡る分布のパラメータ化された族が必要です。ガイドを書き出すことはモデルの構築への近い類推で進めます、ガイド・パラメータは訓練可能であるという主要な違いがあります。これを行なうために pyro.param() を使用して ParamStore でガイド・パラメータを登録します。
softplus = torch.nn.Softplus()
def guide(data):
# 変分パラメータを定義します。
w_loc = torch.randn(1, 1)
# note that we initialize our scales to be pretty narrow
w_log_sig = torch.tensor(-3.0 * torch.ones(1, 1) + 0.05 * torch.randn(1, 1))
b_loc = torch.randn(1)
b_log_sig = torch.tensor(-3.0 * torch.ones(1) + 0.05 * torch.randn(1))
# param store に学習可能な params を登録します
mw_param = pyro.param("guide_mean_weight", w_loc)
sw_param = softplus(pyro.param("guide_log_scale_weight", w_log_sig))
mb_param = pyro.param("guide_mean_bias", b_loc)
sb_param = softplus(pyro.param("guide_log_scale_bias", b_log_sig))
# w と b に対するガイド分布
w_dist = Normal(mw_param, sw_param).independent(1)
b_dist = Normal(mb_param, sb_param).independent(1)
dists = {'linear.weight': w_dist, 'linear.bias': b_dist}
# overload the parameters in the module with random samples
# from the guide distributions
lifted_module = pyro.random_module("module", regression_model, dists)
# sample a regressor (which also samples w and b)
return lifted_module()
ガイド分布の両者のためにガウシアンを選択していることに注意してください。また正値性を確かなものにするために、softplus() 変換を通して各対数スケールを渡してします (正値性を保証するための代替案は exp() 変換です)。
推論
推論を行なうためには確率的変分推論 (SVI) を使用します (SVI のイントロダクションのためには、SVI パート I 参照)。ちょうど非ベイジアン線形回帰におけるように、訓練ループの各反復は勾配ステップを取ります、この場合、MSE 損失の代わりに (SVI に渡す Trace_ELBO オブジェクトを構築することで) ELBO 目的関数を使用するという違いを伴います。
optim = Adam({"lr": 0.05})
svi = SVI(model, guide, optim, loss=Trace_ELBO())
ここで Adam は torch.optim.Adam 回りの薄いラッパーです (議論のためには こちら を参照)。完全な訓練ループは次のようなものです :
def main():
pyro.clear_param_store()
data = build_linear_dataset(N)
for j in range(num_iterations):
# calculate the loss and take a gradient step
loss = svi.step(data)
if j % 100 == 0:
print("[iteration %04d] loss: %.4f" % (j + 1, loss / float(N)))
if __name__ == '__main__':
main()
ELBO 勾配ステップを取るためには SVI の step メソッドを単純に呼び出します。step に渡す data 引数は model() と guide() の両者に渡されることに注意してください。
結果を確認する
学習した変分パラメータを前の結果と比較してみましょう :
for name in pyro.get_param_store().get_all_param_names():
print("[%s]: %.3f" % (name, pyro.param(name).data.numpy()))
サンプル出力 :
[guide_log_scale_weight]: -3.217 [guide_log_scale_bias]: -3.164 [guide_mean_weight]: 2.966 [guide_mean_bias]: 0.941
見て取れるように、私達のパラメータ推定の平均は前に学習した値に非常に近いです。けれども今は、単なる点推定の代わりに、guide_log_scale_weight と guide_log_scale_bias は不確かな推定を提供します(スケールはここでは対数スケールですので、値がより負になれば、幅はより狭くなることに注意してください)。
最後に、新しいテストデータでその予測精度をチェックすることでモデルを評価しましょう。これは点評価として知られます。事後分布から 20 ニューラルネットをサンプリングして、それらをテストデータ上で実行して、それからそれらの予測に渡り平均して正解と比べた予測値の MSE を計算します。
X = np.linspace(6, 7, num=20)
y = 3 * X + 1
X, y = X.reshape((20, 1)), y.reshape((20, 1))
x_data, y_data = torch.tensor(X).type(torch.Tensor), torch.tensor(y).type(torch.Tensor)
loss = nn.MSELoss()
y_preds = torch.zeros(20, 1)
for i in range(20):
# guide does not require the data
sampled_reg_model = guide(None)
# run the regression model and add prediction to total
y_preds = y_preds + sampled_reg_model(x_data)
# take the average of the predictions
y_preds = y_preds / 20
print("Loss: ", loss(y_preds, y_data).item())
サンプル出力 :
Loss: 0.00025596367777325213
Github 上の完全なコードを見てください。
以上