統計モデリング PR

【pyro】ベイズ統計モデルの設計手順メモ

ベイズ統計モデルの設計手順
記事内に商品プロモーションを含む場合があります

こんにちは。

現役エンジニアの”はやぶさ”@Cpp_Learningです。最近は統計モデリングを勉強しています。

今回はベイズ統計モデルの設計手順を整理するために、本記事を書きます。

統計モデルの設計手順

私の場合は以下の手順で統計モデリングをしています。

統計モデルの設計手順
  1. データがどんな数式に従うか考える
  2. データがどんな確率分布に従うか考える
  3. 推論結果を考察
  4. 最初に戻り”より良いモデル”を検討

一言で説明すると『データの背景を考える作業』が統計モデリングだと考えています。

【設計手順❶】データがどんな数式に従うか考える

私が扱う実験データの多くは物理法則に従うので、まず最初に運動方程式を考えます。

ただし 情報不足・知識不足などにより、運動方程式の検討が困難な場合、”あてはまり”が良さそうなモデルを検討します。

sin波

例えば上図のデータなら、以下のモデルだと”あてはまり”が良さそうです。

【モデル】

y = w0 + w1 * x + w2 * x^2 + w3 * x^3

ここでいう”あてはまり”が良いモデルとは、下図の実線(緑)のような”ばらつきの中心”を表現できるモデルという意味です。

pyroでベイズ統計モデリング

【設計手順❷】データがどんな確率分布に従うか考える

データが”あてはまり”の良いモデルを中心に”ばらつく”と考えたとき、データがどんな確率分布に従うか(どう”ばらつく”か)を考えます。

私が扱う実験データの多くは正規分布に従うので、まず最初に以下の式を採用します。

【確率分布】

μ = y, σ = 1.0
Y ~ Normal(μ, σ)

※σ=1.0が最適か否かは、この時点では不明です

【モデル】

y = w0 + w1 * x + w2 * x^2 + w3 * x^3

機械学習では上記のw0~w3に適当な初期値を与え、最適化アルゴリズムにより一意のw0~w3を算出しますが、統計モデリングではw0~w3も何らかの確率分布に従うという考え方ができます。

各パラメータがどんな確率分布に従うかはケースバイケースですが、悩む場合は初手”正規分布”で良いと思います(あとから変更できるので)。

【パラメータ】

μ = 0, σ = 10
W0 ~ Normal(μ, σ)
W1 ~ Normal(μ, σ)
W2 ~ Normal(μ, σ)
W3 ~ Normal(μ, σ)

※μ = 0, σ = 10が最適か否かは、この時点では不明です

もちろんドメイン知識や仮説から「W0 = 1の固定値を採用」などもOKです。

  • 機械学習のように各パラメータや結果を一意に決めるのではなく、不確実性を考慮した設計・推論ができるは「統計モデリング」の強みです
  • 全てのパラメータをコンピュータに算出させる必要はなく、ドメイン知識や仮説を組み込んだモデリングが重要です(自論)

【設計手順➌】推論結果を考察

事後分布は?不確実性は?”あてはまり”の良さは?などの推論結果を考察し、σも正規分布に従うと考えるなら、以下のように設計を変更します。

【確率分布】

μ = y
σ ~ Normal(0, 1.0)
Y ~ Normal(μ, σ)

※時間が許す限り何度も設計変更すれば良いと考えています

【設計手順❹】最初に戻り、より良いモデルを検討

最終的に設計したモデルは以下の通りです。

【パラメータ】

μ = 0, σ = 10
W0 ~ Normal(μ, σ)
W1 ~ Normal(μ, σ)
W2 ~ Normal(μ, σ)
W3 ~ Normal(μ, σ)

【モデル】

y = w0 + w1 * x + w2 * x^2 + w3 * x^3

【確率分布】

μ = y
σ ~ Normal(0, 1.0)
Y ~ Normal(μ, σ)

ここで設計した確率分布が事前分布になります。最終的には実験データの分布を考慮し、ベイズの定理やMCMCなどのアルゴリズムを活用して、事後分布を推論します。

この推論を自動化するために、コンピュータとPPLを使います。

あとは仮説検証・検証結果・追加情報・有識者からのアドバイスなど、あらゆる情報を考慮して、上記設計のアップデートを繰り返します

以上がベイズ統計モデルの設計手順です。

実践!Pyroでベイズ統計モデリング

以降からはPyroによるベイズ統計モデリングを実践します。以下の記事を先に読んでおくと本記事の内容(ソースコード含む)をスッと理解できると思います。

NumPyroでベイズ回帰モデルをつくる
【ベイズ統計モデリング入門】NumPyroで回帰モデルをつくるこんにちは。 現役エンジニアの”はやぶさ”@Cpp_Learningです。最近は統計モデリングを勉強しています。 今回は確率...
Pyroと変分推論でベイズ回帰モデルをつくる
Pyroと変分推論でベイズ回帰モデルをつくるこんにちは。 現役エンジニアの”はやぶさ”@Cpp_Learningです。最近は統計モデリングを勉強しています。 今回は確率...

以降からコード書いていきます。

Import

まずはimportから

※描画(データを可視化するときの)スタイルは自由に変更してください

サンプルデータ

適当なサンプルを用意します。

sin波

このデータが今回の対象です。

説明変数と目的変数 -torch.tensorに変換

PyroのバックエンドがPyTorchなので、扱うデータをtorch.tensorに変換します。

ベイズ統計モデルの設計

今回はSinにノイズを混ぜたデータ(y = sin(x) + ε)に対し、回帰モデルを設計します。

真のモデルは y = sin(x) ですが、「1周期分のデータからsinを特定するのは困難」という前提で、今回は y = w0 + w1 * x + w2 * x^2 + w3 * x^3 だと”あてはまり”が良いと考えて設計しました。

MCMCの実行

今回はMCMCで事後分布を推論します。

変分推論でも良いのですが、下記の参考記事にある各手法の良い点/悪い点を考慮し、十分なサンプル数や演算時間を確保できるなら、初手MCMCで良いと思いました。

事後分布の可視化

以下のコードで各パラメータの基本統計量を出力できます。

mean std median 5.0% 95.0% n_eff r_hat
w0 -0.26 0.04 -0.26 -0.33 -0.19 514.54 1.00
w1 1.98 0.06 1.98 1.88 2.07 403.04 1.00
w2 0.90 0.02 -0.90 -0.94 -0.86 409.40 1.00
w3 0.09 0.00 0.09 0.09 0.10 432.71 1.00
sigma 0.11 0.01 0.11 0.10 0.13 778.02 1.00

MCMCで算出されたパラメータは以下のコードで取得できます。

以下のコードで事後分布を可視化します。

事後分布

W3 や σ は”ばらつき”が小さいので、再設計するときは固定値にしても良いかもしれません。

推論

以下のコードで X=0~6.5(50サンプル) のときの Y を予測します。

以下のコードで予測結果を可視化します。

pyroでベイズ統計モデリング

今回は上記の予測に満足したことにしますが…

  • x > 6.5 のデータを追加
  • 説明変数を追加
  • 有識者からのアドバイスを考慮してモデルを見直す
  • サンプル数を追加後に確率分布を見直す

など

上記を考慮し、繰り返し統計モデリングにトライすることが多いです。

以上でベイズ統計モデリング実践も終了です。

スポンサーリンク

まとめ

頭の中を整理しながら、ベイズ統計モデルの設計手順をドキュメント化しました。

私の統計モデリングに関する経験値や知見が増えたタイミングで、より良い設計手順を再検討し、本記事の内容も更新できればと思います。

なので…

くるる
くるる
僕が考えた最強の統計モデリング

みたいな知見がある人は情報共有して頂けると嬉しいです。

はやぶさ
はやぶさ
統計モデリングに関する情報共有をしながら、一緒に成長できると嬉しいです。よろしくお願いします。

以下 統計モデリング関連のオススメ書籍

PICK UP BOOKS

  • 数理モデル入門
    数理モデル
  • Jetoson Nano 超入門
    Jetoson Nano
  • 図解速習DEEP LEARNING
    DEEP LEARNING
  • Pythonによる因果分析
    Python