回帰問題の定式化

定式化

問題を整理する

手元にあるデータ\(X=\{\boldsymbol{x}_n\}_{n=1}^N\), \(Y=\{\boldsymbol{y}_n\}_{n=1}^N\)に対して, 回帰問題を以下のように推測問題として定式化します. \(X\)の次元を\(D_x\), \(Y\)の次元を\(D_y\)とします. \(X\), \(Y\)は真の分布\(q(\boldsymbol{x},\boldsymbol{y})\)からのサンプルとします. \(q(\boldsymbol{y}\mid\boldsymbol{x})\)または, 入力から出力への対応\(\mathbb{R}^{D_x}\to\mathbb{R}^{D_y}\)をデータから推測します. モデル\(p(\boldsymbol{y}\mid\boldsymbol{x},\boldsymbol{w})\)を用いて, 何らかの方法で計算した予測分布\(p^*(\boldsymbol{y}\mid\boldsymbol{x})\)によって, 真の分布を推測します. モデルは以下のように書けるとします. ただし, \(g\)は適当な分布の確率(密度)関数とします.

また, 以下の予測分布\(p^*\)の平均によって入力から出力への対応\(\mathbb{R}^{D_x}\to\mathbb{R}^{D_y}\)を推測し, 回帰関数と呼びます.

推測方法

推測方法として, 最尤推測, MAP推測, 平均プラグイン推測, Bayes推測を考えます.

【最尤推測】

---

最尤推測では, 次式で定義される最尤推定量をモデルに代入して予測分布とします.

\begin{equation} \boldsymbol{w}_{\mathrm{ML}} = \underset{w}{\mathrm{argmax}}\prod_{n=1}^Np(\boldsymbol{y}_n\mid\boldsymbol{x}_n,\boldsymbol{w}) = \underset{w}{\mathrm{argmin}}\left\{-\sum_{n=1}^N\log p(\boldsymbol{y}_n\mid\boldsymbol{x}_n,\boldsymbol{w})\right\} \end{equation}

特に, モデル平均が\(f(\boldsymbol{x},\boldsymbol{w})\)ならば\(f^*(\boldsymbol{x})=f(\boldsymbol{x},\boldsymbol{w}_{\mathrm{ML}})\)です. さらに, \(g\)が平均がゼロベクトルの正規分布ならば, 次式が成り立ちます.

\begin{equation} \boldsymbol{w}_{\mathrm{ML}} = \underset{w}{\mathrm{argmin}}\sum_{n=1}^N\|\boldsymbol{y}_n-f(\boldsymbol{x}_n,\boldsymbol{w})\|^2 \end{equation}

これは最小二乗法です.

【MAP推測】

---

MAP推測では, 次式で定義されるMAP推定量をモデルに代入して予測分布とします.

\begin{equation} \boldsymbol{w}_{\mathrm{MAP}} = \underset{w}{\mathrm{argmax}}\ p(\boldsymbol{w}\mid X,Y) = \underset{w}{\mathrm{argmin}}\left\{-\sum_{n=1}^N\log p(\boldsymbol{y}_n\mid\boldsymbol{x}_n\boldsymbol{w})-\log p(\boldsymbol{w})\right\} \end{equation}

特に, モデル平均が\(f(\boldsymbol{x},\boldsymbol{w})\)ならば\(f^*(\boldsymbol{x})=f(\boldsymbol{x},\boldsymbol{w}_{\mathrm{MAP}})\)です. さらに, \(g\)と\(\boldsymbol{w}\)も事前分布が平均ゼロベクトルの正規分布ならば, ある\(\alpha\)が存在して, 次式が成り立ちます.

\begin{equation} \boldsymbol{w}_{\mathrm{MAP}} = \underset{w}{\mathrm{argmin}}\left\{\sum_{n=1}^N \|\boldsymbol{y}_n-f(\boldsymbol{x}_n,\boldsymbol{w})\|^2+\alpha\|\boldsymbol{w}\|^2\right\} \end{equation}

これはridge回帰です.

【平均プラグイン推測】

---

平均プラグイン推測では, 事後平均\(\bar{\boldsymbol{w}}\)をモデルに代入して予測分布とします.

特に, モデル平均が\(f(\boldsymbol{x},\boldsymbol{w})\)ならば\(f^*(\boldsymbol{x})=f(\boldsymbol{x},\bar{\boldsymbol{w}})\)です.

【Bayes推測】

---

Bayes推測では, 次式で予測分布を定めます.

\begin{equation} p^*(\boldsymbol{y}\mid\boldsymbol{x}) = \int p(\boldsymbol{y}\mid\boldsymbol{x},\boldsymbol{w})p(\boldsymbol{w}\mid X,Y)\mathrm{d}\boldsymbol{w} \end{equation}

特に, モデル平均が\(f(\boldsymbol{x},\boldsymbol{w})\)ならば, 回帰関数は以下のようになります.

\begin{equation} f^*(\boldsymbol{x}) = \int f(\boldsymbol{x},\boldsymbol{w})p(\boldsymbol{w}\mid X,Y)\mathrm{d}\boldsymbol{w} \end{equation}

実験

\(D_x=D_y=1\)とします. \(g\)を平均をゼロベクトル, 共分散行列を\(\sigma^2I\)とする正規分布とします. 関数列\(\{\phi_i\}_{i=1}^{d_w}\)を次式で定めます.

このとき, \(f\)を次のように定めます.

事前分布は平均がゼロベクトル, 精度が\(\lambda_wI\)の正規分布とします. また, \(\Phi\)を計画行列とします. このとき, 各推定値は次式を求めます. \(\alpha=\sigma^2\lambda_w\)とします.

Bayes推測の回帰関数は, 以下の予測平均です.

実はこの回帰関数は, MAP推測と平均プラグイン推測の回帰関数と一致します.

次のように, データ点が\(N=25\)点あるとします.

【コード2の実行結果】

4つの手法を実行した結果を下図に示します. MAP推測, 平均プラグイン推測, Bayes推測は一致しています.

【コード4の実行結果】

コード

#mathematics
using LinearAlgebra

#statistics
using Random

#visualize
using Plots
pyplot()

#macros
using UnPack

#sample X and Y
Random.seed!(42)
F(x) = 0.7*(sin(3*x) + cos(2*x) + sin(5*x))
N = 25
X = 6*(rand(N).-0.5)
Y = F.(X) + 0.5*randn(N)

#plot the data
fig1 = plot(X, Y, st=:scatter, xlabel="x", ylabel="y", title="data", label="data", 
    markershape=:x, markerstrokewidth=1, markersize=10, markercolor=:blue)
savefig(fig1, "figs-RP/fig1.png")

#basis function 
ϕi(x,μi,σi) = exp(-(x-μi)^2/σi/σi/2)
function ϕ(x,params)
    @unpack μs,σs,dw = params
    vals = ones(dw)
    for i in 1:dw-1
        vals[i+1] = ϕi(x,μs[i],σs[i])
    end
    return vals
end

#design matrix
function design_matrix(X,params)
    N = length(X)
    mat = zeros(N,dw)
    for n in 1:N
        mat[n,:] = ϕ(X[n],params)
    end
    return mat
end

#model
f(x,wvec,params) = dot(wvec, ϕ(x,params))

#define model
σ = 1
λw = 0.01
dw = 8
μs = [-3,-2,-1,0,1,2,3]
σs = ones(dw-1)
params = (σ=σ, λw=λw, μs=μs, σs=σs, dw=dw, N=N)

#matrix, vector
Φ = design_matrix(X,params)
ΦTΦ = Φ' * Φ
ΦTy = Φ' * Y

#maximum liklihood
wML = ΦTΦ\ΦTy

#MAP estimate
α = σ^2*λw
wMAP = (ΦTΦ+α*I(dw))\ΦTy

#mean plugin
wmean = (ΦTΦ/σ/σ+λw*I(dw))\(ΦTy/σ/σ)

#Bayes estimate
function μpred(x,params)
    @unpack σ,λw,dw = params
    tmp1 = (ΦTΦ+ϕ(x,params)*ϕ(x,params)'+σ^2*λw*I(dw))\ΦTy
    tmp2 = (ΦTΦ+ϕ(x,params)*ϕ(x,params)'+σ^2*λw*I(dw))\ϕ(x,params)
    return dot(ϕ(x,params), tmp1)/(1-dot(ϕ(x,params), tmp2))
end

#visualize
fig2 = plot(X,Y,st=:scatter,xlabel="x",ylabel="y",title="fitting",label="data",ylim=[-2,2],legend=:bottomright,
    markershape=:x, markerstrokewidth=1, markersize=10, markercolor=:blue)
plot!(-3:0.1:3, x->f(x,wML,params), label="ML", ls=:dot, lw=2)
plot!(-3:0.1:3, x->f(x,wMAP,params), label="MAP", ls=:dash, lw=2)
plot!(-3:0.1:3, x->f(x,wmean,params), label="plug-in", ls=:solid, lw=1, alpha=0.5)
plot!(-3:0.1:3, x->μpred(x,params), label="Bayes", ls=:dashdot, lw=2)
savefig(fig2, "figs-RP/fig2.png")