これなら分かるBayes統計

記事の内容


  • 例題: 線形回帰問題
    • 問題設定
    • 【余談】最尤推定
  • Bayes統計の考え方
    • 推定・予測の流れ
    • 事前分布の設定
    • 事後分布の実現
    • 事後分布の利用
    • 【余談】モデル選択
    • 【余談】怪しげなBayes統計
    • 【余談】 統計物理との関連
  • 推定をやってみよう!
    • 確率モデルの設定
    • 事前分布の設定
    • 尤度の計算
    • 事後分布の計算
    • 事後分布の観察
    • 予測分布の計算
    • 予測分布の観察
  • ソースコード
    • 共通のコード
    • アニメーションのコード

今回は, 線形回帰を例題に, Bayes統計学を概観していきましょう.

例題: 線形回帰問題

問題設定

次のようなデータを考えます. データ点が20点あります.

このデータを用いて予測を行います. 次に同じような入力が与えられたとき, 出力を予測します. データを見ると, 直線的な関係が読み取れます. ということで, このデータをよく説明する直線を求めましょう. よって, 線形回帰問題として予測を行います.

直線を求めたいので, 式にしておきましょう.

\[ y \simeq w_0 + w_1 x \]

もしこの直線がデータをよく近似するならば, データ点のy座標は次のように近似できます.

\[ y_i \simeq w_0 + w_1 x_i \quad (i=1,\cdots,N) \]

ここで, データ点の数を\(N\), データ点を\((x_i, y_i)\)とします. しかし, データ点は厳密に直線上にあるわけではないので, その差を\(\epsilon_i\)とおきます. よって, \[ y_i = w_0 + w_1 x_i + \epsilon_i, \quad (i=1,\cdots,N) \] となります. さらに次のような仮定をおきます.

\[ \epsilon_i \sim \mathrm{N} (0,\beta^{-1}) \]

それでは, \(w_0, w_1\)を求めましょう. これが求まれば, 直線が特定できます. 直線が特定できれば, その式をもとに予測ができます. ここではBayes統計を用いて推定します.

【余談】最尤推定

ここではBayes統計を用いますが, 他にも方法はあります. 最尤推定です. 最尤推定とは, 尤度関数を最大にするパラメータ値を推定値とする方法です. 尤度関数はパラメータのデータへの適合具合を表します. 今回は最尤推定は用いません.

Bayes統計の考え方

推定・予測の流れ

統計学では, 観測値が真の分布からの実現値であると考えます. この真の分布を推測することが目標です.推定・予測の流れは次の通りです.

  • 確率モデルを設定する(このモデルで真の分布を推測したい).
  • 事前分布を設定する(決め方は様々存在する).
  • 尤度関数を計算する(確率モデルから得られる).
  • 事後分布を計算する(パラメータの推定値が得られる).
  • 予測分布を計算する(予測まで行う場合には計算する).

事前分布の設定

事前分布の決め方には様々あります. 現時点で「これが決定版」というものは無さそうです... . 以下のような決め方があります.

  • 計算しやすいように共役事前分布を用いる.
  • パラメータ値に関する事前知識を利用する.
  • 無情報事前分布や弱情報事前分布を用いる.

事後分布の計算

事後分布は普通, 次のBayesの定理から計算します(Bayesの定理を表に出さない議論もあります).

\[ p(w|D) \propto p(D|w)p(w) \]

\(p(w|D) \)は事後分布, \(p(D|w)\)は尤度関数, \(p(w)\)は事前分布を表します. また, \[ D = \left\{ (x_i, y_i) | i=1,\cdots, N \right\} \] をデータとします.

ただし, 事後分布が複雑で扱いにくい場合もあります. この場合, MCMCや変分推論などの近似推論法を用います.

事後分布の利用

事後分布は1つの決まった確定値というわけではありません. パラメータに関する様々な情報を含んでいるため, 利用法も様々あります. 例えば, 次のような方法があります.

  • 事後分布の平均をパラメータの推定値とする.
  • 事後分布を最大にするパラメータ値を推定値とする.
  • MCMCのサンプル平均を推定値とする.

【余談】モデル選択

データの解析に対して, 確率モデルは複数考えられます. これらのモデルのうちどれが最適かを選ぶ方法があります. 次の3つが代表的なものです.

  • 周辺尤度を最大化する.
  • 情報量規準(AIC, WAIC)を用いる.
  • クロスバリデーション(LOOCV, k-fold-CV)を用いる.

【余談】怪しげなBayes統計

数理的なBayes統計と以下の言説には関係がありません. この記事でも触れません.

  • 頻度主義は間違っており, Bayes主義が正しい.
  • Bayes統計は間違っており, 頻度主義が正しい.
  • 哲学的にBayes統計は間違っている(または正しい).
  • 人間の脳はBayes的である.

【余談】 統計物理との関連

Bayes統計と統計物理には密接な関係があります. 分配関数などの用語が普通に使用されます. また, MCMCの一つ, HMC(Hamiltonianモンテカルロ法)では, Hamilton系も出現し, 解析力学とも関係があります.

推定をやってみよう!

例題に戻ります. 推定・予測を実践します.

確率モデルの設定

確率モデルは, 例題での仮定から自然に出てきます. 先ほどの例題では, \[ y = w_0 + w_1 x + \epsilon, \quad \epsilon \sim \mathrm{N} (0,\beta^{-1}) \] としました. したがって, \[ y \sim \mathrm{N} (w_0 + w_1 x,\beta^{-1}) \] を確率モデルとします.

事前分布の設定

事前分布は, 正規分布にします.

\[ \begin{bmatrix} w_0 \\ w_1\end{bmatrix} \sim \mathrm{N}(0,\alpha^{-1}I) \]

尤度の計算

尤度関数を計算します. 尤度関数は確率モデルにデータを代入したものの積で定義されます. 今回の例題の場合は次のようになります.

\[ p(D|w_0, w_1,\beta) = \prod_{i=1}^{N}p(y_i|x_i, w_0, w_1) \propto \exp\left[ -\frac{\beta}{2} \sum_{i=1}^{N}\left\{ y_i-(w_0+w_1x_i)\right\}^2\right] \]

事後分布の計算

Bayesの定理を用いて事後分布を計算します. 計算過程は省略しますが, 次のようになります.

\begin{align} \begin{bmatrix}w_0\\w_1\end{bmatrix}|D &\sim \mathrm{N}(\mu_{\mathrm{post}}, \Sigma_{\mathrm{post}}) \\ \Sigma_{\mathrm{post}}^{-1} &= \begin{bmatrix} \alpha + \beta N & \beta \sum_{i=1}^{N} x_i \\ \beta \sum_{i=1}^{N} x_i & \alpha + \beta \sum_{i=1}^{N}x_i\\ \end{bmatrix}\\ \mu_{\mathrm{post}} &= \beta \Sigma_{\mathrm{post}}\begin{bmatrix} \sum_{i=1}^{N} y_i \\ \sum_{i=1}^{N} x_i y_i \end{bmatrix}\\ \end{align}

事後分布の観察

事後分布の様子を可視化します.

【コード4の実行結果】

赤い線は真の直線です. 最初はデータ点が1つだけです. この状態で事後分布(の平均と共分散行列)を計算します. 計算した平均と共分散行列をもとに, 正規分布からサンプルします. つまり, 事後分布からのサンプルです. サンプルを10個得ます. これらのサンプルを用いて, 直線を描画したものがグレーの直線です. そして, データ点を1個追加して同じことを繰り返します. つまり, データ点を追加するごとに10本の直線をサンプルします. データ点が追加されるほど, 各サンプルが真の直線に近づきます. つまり, データ点が増えるほど推定の精度が上がります.

事後分布を別の観点から眺めてみましょう.

【コード5の実行結果】

黒い部分は事後密度関数の値の大きい部分を表します. このアニメーションでも, データ点を追加するごとに事後分布を計算します. データ点が増える度に黒い部分の範囲が小さくなります. つまり, データ点が増えるほど, 不確かさが減少します.

予測分布の計算

予測分布を計算します. 予測分布は次式で定義されます.

\[ p(y^*|x^*) = \int p(y|w_0, w_1)p(w_0, w_1|D)dw_0dw_1 \]

これを計算すると, 次のような正規分布になります. 計算過程は省略します.

\begin{align} y^* | x^* &\sim \mathrm{N}(m,\sigma^2) \\ m &= \mu^{\mathrm{T}}\begin{bmatrix}1 \\ x^*\end{bmatrix} \\ \sigma^2 &= \frac{1}{\beta} + \begin{bmatrix}1 \\ x^*\end{bmatrix}^{\mathrm{T}} \Sigma_{\mathrm{post}} \begin{bmatrix}1 \\ x^*\end{bmatrix} \end{align}

x*は新しい入力で, y*はそれに対する予測値です.

予測分布の観察

予測分布を観察しましょう.

【コード6の実行結果】

赤い直線は真の直線です. このアニメーションでも, データ点を追加するごとに予測分布を計算します. 青線は予測分布の平均です. また, 緑色部分は「平均±標準偏差」の幅を表します. データ点がない部分は標準偏差が大きく, データ点が追加されると標準偏差が小さくなります. つまり, データが多くなるほど不確かさが減少します. 平均は真の直線に近づいているので, 予測の精度は段々上がります.

ソースコード

共通のコード

【Juliaコード1; インポート】
using Plots
using LinearAlgebra
using Random
【Juliaコード2; データの作成】
Random.seed!(0);
#サンプル数
n_samps = 10;

#ノイズの精度と事前分布の精度
α = 1.5;
β = 1/0.3^2;

#全データ数
n_data = 20;

#パラメータの真値
w₀_true = -0.1;
w₁_true = 0.6;

#訓練データとテストデータの作成
X_train = 2*rand(n_data).-1 ;
y_train = w₀_true .+ w₁_true * X_train .+ 1/sqrt(β) * randn(n_data);
【Juliaコード3; 関数の作成】
#直線
function line(x::Array{Float64,1}, a::Float64, b::Float64)
    return a*x .+ b
end

#多変量正規分布からのサンプル
function gaussian_sample(mean::Array{Float64,1}, cov::Array{Float64,2}, n_samps::Int64)
    """
    Cholesky分解を用いて, 正規分布N(mean, cov)からサンプルする.
    """
    C = cholesky(cov).L;
    samps = C*randn(2,n_samps) .+ mean;
    return samps
end

#事後パラメータの計算
function calc_post_params(x::Array{Float64,1}, y::Array{Float64,1}, α::Float64, β::Float64)
    #データ数
    N = length(x);

    #和, 二乗和, 積の和
    x_sum = sum(x);
    y_sum = sum(y);
    x_sq_sum = sum(x.^2);
    xy_sum = sum(x.*y);

    #精度行列
    SN_inv = [
        α+β*N β*x_sum;
        β*x_sum α+β*x_sq_sum
    ];

    #平均
    mN = SN_inv\[β*y_sum; β*xy_sum];
    return SN_inv\I(2), mN
end

#予測分布のパラメータ
function calc_pred_params(x::Array{Float64,1}, y::Array{Float64,1}, α::Float64, β::Float64)
    #事後パラメータ
    cov, mN = calc_post_params(x,y,α,β);

    #平均
    pred_mean(x) = dot(mN, [1;x]);
    pred_var(x) = 1/β+dot([1;x], cov*[1;x]);
    return pred_mean, pred_var
end

#正規分布
function normal(x₁::Float64, x₂::Float64, μ::Array{Float64,1}, Σ::Array{Float64,2})
    #ベクトル化
    x = [x₁,x₂];

    return exp(-dot((x.-μ), Σ\(x.-μ))^2/2)
end

アニメーションのコード

【Juliaコード4; 事後分布のアニメーション】
Random.seed!(27);

#アニメーション
anim = Animation();

#データ数を増やしながら逐次的に計算
for N in 1:n_data
    #訓練データの途中までを用いて事後分布を計算する
    xN = X_train[1:N];
    yN = y_train[1:N];

    #真値のプロット
    true_line(x) = line(x, w₁_true, w₀_true);
    y_true = true_line([-1.,1.]);
    true_line_plot = plot([-1.,1.], y_true, color="red", label=["true"], xlabel="x", ylabel="y", ylim=[-1,1]);

    #N番目までのデータから計算する事後パラメータ
    cov_post, mean_post = calc_post_params(xN, yN, α, β);
    try
        w_samps = gaussian_sample(mean_post, cov_post, n_samps);

        #各サンプルから10本の直線をプロット
        for k in 1:n_samps
            #パラメータの取得
            params = w_samps[:,k];

            #傾きと切片
            w₀, w₁ = params[1], params[2];

            #プロット
            sample_line(x) = line(x, w₁, w₀);
            sample_line_plot = plot!([-1.,1.], sample_line([-1.,1.]), color="gray", legend=false);

            #データの散布図
            data_scatter = scatter!(xN, yN, color="blue", legend=false);

            #アニメーションに追加
            frame(anim, sample_line_plot);
            frame(anim, data_scatter);
        end
    catch
        println(N);
    end
end

#保存
#gif(anim, "Bayes-learning.gif", fps=10);
【Juliaコード5; 事後分布のアニメーション】
Random.seed!(27);

#アニメーション
anim = Animation();

x₁_array = -1.0:0.02:1.;
x₂_array = -1.0:0.02:1.;

#真値のプロット
scatter(76, 36, markercolor=:red, legend=false, xlim=[1,50], ylim=[1,50]);

for N in 1:n_data
    #訓練データ
    xN = X_train[1:N];
    yN = y_train[1:N];

    #N番目までのデータから計算される事後パラメータ
    cov_post, mean_post = calc_post_params(xN, yN, α, β);

    #事後分布の確率密度関数の計算
    pdf = [normal(x₁, x₂, mean_post, cov_post) for x₁ in x₁_array, x₂ in x₂_array];
    pdf_heatmap = heatmap(pdf, c=:Greys_9);

    #アニメーションに追加
    frame(anim, pdf_heatmap);
end

#保存
#gif(anim, "Bayes-learning-pdf.gif", fps=5);
【Juliaコード6; 予測分布のアニメーション】
Random.seed!(27);

#アニメーション
anim = Animation();

#データ数を増やしながら逐次的に計算
for N in 1:n_data
    #訓練データの途中までを用いて事後分布を計算する
    xN = X_train[1:N];
    yN = y_train[1:N];

    #真値のプロット
    true_line(x) = line(x, w₁_true, w₀_true);
    y_true = true_line([-1.,1.]);
    true_line_plot = plot([-1.,1.], y_true, color="red", labels="true", xlabel="x", ylabel="y", ylim=[-1,1]);

    #データのプロット
    data_scatter = scatter!(xN, yN, legend=false, labels="data");

    #予測分布
    pred_mean, pred_var = calc_pred_params(xN, yN, α, β);

    #予測平均のプロット
    pred_mean_line = pred_mean.(-1:0.01:1);
    pred_std = sqrt.(pred_var.(-1:0.01:1));
    pred_plot = plot!(-1:0.01:1, pred_mean_line, ribbon=pred_std, linecolor=:blue, labels="predictive mean");

    #アニメーションを保存
    frame(anim, pred_plot);
    frame(anim, data_scatter);
end

#保存
#gif(anim, "Bayes-predictive.gif", fps=10);