30  Bias and Variance

Bias and Variance in Prediction

未知のパラメーターの推定量に焦点を当てた統計的推定とことなり，教師あり機械学習の予測では $x$ が与えられたときに $y$ を良く予測できる関数 $h$

$$\begin{array}{rl}& y=h(x)+ϵ\\ & \text{where }\mathbb{E}[ϵ|x]=0\end{array}$$

を学習することを目指します．予測関数の良さについては決定理論(decision theory)で議論されるものですが，一旦L2距離で良さを判断する場合，良い関数は $\mathbb{E}[y|x]$ となります．この関数の学習を目指して，

$$D=\{x^{(i)},y^{(i)}{\}}_{i=1}^n$$

をtraining datasetとして，${\hat{f}}_D$ という関数を構築したとします．母集団からランダムサンプリング された $S$ にはノイズ ${ϵ}_i$ が含まれており，$S$ を用いて学習した ${\hat{f}}_D$ は確率的に変動する側面があります．その上で同じ母集団からi.i.dでサンプリングされた未知のデータセット $\{x^{\ast },y^{\ast }\}$ においても $x$ を特徴量として用いることで 精度良く $y$ を予測することができること = 汎化誤差(generalization error)が小さい ${\hat{f}}_D$ が望ましい予測関数となります．

汎化誤差を $y$ と ${\hat{f}}_D$ のL2距離と考え，$x^{\ast }$ で条件づけた場合，

$$\begin{array}{}\text{(30.1)}& \begin{array}{rl}& \mathrm{MSE}({\hat{f}}_D)\\ & ={\mathbb{E}}_{D,y^{\ast }}[(y^{\ast }-{\hat{f}}_D(x^{\ast }){)}^2]\\ & ={\mathbb{E}}_{D,y^{\ast }}[(h(x^{\ast })+ϵ-{\hat{f}}_D(x^{\ast }){)}^2]\\ & ={\mathbb{E}}_{y^{\ast }}[{ϵ}^2]+{\mathbb{E}}_D[(h(x^{\ast })-{\hat{f}}_D(x^{\ast }){)}^2]\\ & ={\mathbb{E}}_D[{ϵ}^2]+\underset{=\text{bias and constant}}{\underbrace{{\mathbb{E}}_D[h(x^{\ast })-{\hat{f}}_D(x^{\ast }){]}^2}}+\mathrm{Var}(h(x^{\ast })-{\hat{f}}_D(x^{\ast }))\\ & =\underset{\text{irreducible error}}{\underbrace{\mathrm{Var}(ϵ)}}+\text{Bias}({\hat{f}}_D{)}^2+\mathrm{Var}({\hat{f}}_D)\end{array}\end{array}$$

上式の第３項のVarianceは予測関数 $f$ の特定のtraing data集合 $D$ の選び方に関する敏感さを表していると解釈されます．

▶  regession with gauss basis function

上記の議論を踏まえると，Bias, Varianceをバランス良く小さくするモデルがMSEの意味で良い予測モデルということになります．

経験則的に，柔軟性のある複雑なモデルはBiasは小さいがVarianceが大きく，逆に柔軟性の低いモデルはBiasが大きいがVarianceが小さいことが知られています．

以下の例では，

$$\begin{array}{rl}{y}_i& =\sin (2\pi x_i)+{ϵ}_i\\ {ϵ}_i& \stackrel{\mathrm{iid}}{\sim }N(0,{0.5}^2)\\ x_i& \stackrel{\mathrm{iid}}{\sim }\mathrm{Uniform}(0,1)\end{array}$$

というData Generating Processのもとで $N_{train}=25$ からなるsampleを生成し，そのsample生成を100回繰り返す処理を考えます．つまり，

• それぞれのsample sizeは $N_{train}=25$
• sample数は $L=100$

ということです．24個のガウス基底を基底関数と用いたRidge Linear Regressionをそれぞれのsampleに対して実行し，その結果得られた各予測モデル(100個の予測モデル)の予測結果についてPythonで確認してみました．

なお，Ridge penaltyは $[0.1,1.5,15]$ と恣意的に設定しています．なおTest datasetのsample sizeは $N_{test}=1000$ としています

import numpy as np
from sklearn.linear_model import Ridge
import matplotlib.pyplot as plt

import regmonkey_style.stylewizard as sw

sw.set_templates("regmonkey_line")


## create data function
def true_func(x):
    return np.sin(2 * np.pi * x)


def generate_data(func, sample_size, std, domain=[0, 1]):
    x = np.random.uniform(domain[0], domain[1], sample_size)
    t = func(x) + np.random.normal(scale=std, size=x.shape)
    return x, t


def create_gauss_basis(x, gauss_mean, gauss_var):
    if x.ndim == 1:
        x = x[:, None]
    else:
        assert x.ndim == 2
    basis = [np.ones(len(x))]
    for m in gauss_mean:
        gauss_transfomred = np.exp(-0.5 * np.square(x - m) / gauss_var)
        basis.append(gauss_transfomred.flatten())
    return np.asarray(basis).transpose()


## Generate Data
np.random.seed(1212)
gauss_basis_mean = np.linspace(0, np.pi / 2, 24)
gauss_basis_var = 0.01
penalty_list = [15, 1.5, 0.1]

x_test = np.linspace(0, 1, 1000)
X_model_test = create_gauss_basis(x_test, gauss_basis_mean, gauss_basis_var)
y_test = true_func(x_test)

x_train, y_train = generate_data(true_func, 25 * 100, 0.5)
x_train, y_train = x_train.reshape(100, 25), y_train.reshape(100, 25)


for penalty in penalty_list:
    y_list = []
    fig, axes = plt.subplots(1, 2, figsize=(9, 3))
    for i in range(100):
        X_model_train = create_gauss_basis(
            x_train[i, :], gauss_basis_mean, gauss_basis_var
        )
        clf = Ridge(alpha=penalty)
        clf.fit(X_model_train, y_train[i, :])
        y_pred = clf.predict(X_model_test)
        y_list.append(y_pred)

    for y in y_list[:10]:
        axes[0].plot(x_test, y, alpha=0.8, linewidth=1, label=f"penalty: {penalty:.2f}")
    axes[0].set_title(f"model regression with penalty: {penalty:.1f}", fontsize=10)
    for ax in axes:
        ax.set_ylim(-1.9, 1.9)
    axes[1].plot(x_test, y_test, label="true")
    axes[1].plot(
        x_test, np.asarray(y_list).mean(axis=0), linestyle="--", label="pred-mean"
    )
    axes[1].legend()
plt.show()

Bias-Variance Decompositionをシミュレーション標本で計算する場合は Equation 30.1 をシミュレーション標本と対応させれば良いので

$$\begin{array}{rl}\bar{f}(x)& =\frac{1}{L}\sum _l^L{\hat{f}}^{(l)}(x)\\ {\mathrm{Bias}}^2& =\frac{1}{N_{test}}\sum _n\{\bar{f}(x_n)-h(x_n){\}}^2\\ \mathrm{Variance}& =\frac{1}{N_{test}}\sum _n\frac{1}{L}\sum _l{\{{\hat{f}}^{(l)}(x_n)-\bar{f}(x_n)\}}^2\end{array}$$

この考えに基づいて ridge penalty に対するBias, Variance, MSEをそれぞれplotすると以下のようにトレードオフ関係が確認できます

"""
 * Computes the bias-variance decomposition for a given list of penalty values.
 * 
 * @param penalty_list - A list of penalty values.
 * @returns An array of bias and variance values corresponding to each penalty value.
"""
from multiprocessing import Pool

penalty_list = np.linspace(0.1, 8, 500)
bias_list = []
variance_list = []
h_test = true_func(x_test)

x_train, y_train = generate_data(true_func, 25 * 100, 0.5)
x_train, y_train = x_train.reshape(100, 25), y_train.reshape(100, 25)


def compute_bias_variance(penalty):
    y_list = []
    for i in range(100):
        X_model_train = create_gauss_basis(
            x_train[i, :], gauss_basis_mean, gauss_basis_var
        )
        clf = Ridge(alpha=penalty)
        clf.fit(X_model_train, y_train[i, :])
        y_pred = clf.predict(X_model_test)
        y_list.append(y_pred)
    y_list = np.asanyarray(y_list)
    bias_term = np.mean((h_test - np.mean(y_list, axis=0)) ** 2)
    variance_term = np.mean(np.mean((y_list - np.mean(y_list, axis=0)) ** 2, axis=0))

    return bias_term, variance_term


with Pool() as pool:
    results = pool.map(compute_bias_variance, penalty_list)

bias_list, variance_list = zip(*results)

import plotly.graph_objects as go

# Plot bias vs penalty
fig = go.Figure()
fig.add_trace(go.Scatter(x=penalty_list, y=bias_list, mode="lines", name="Bias"))
fig.add_trace(
    go.Scatter(x=penalty_list, y=variance_list, mode="lines", name="Variance")
)
fig.add_trace(
    go.Scatter(
        x=penalty_list,
        y=np.asarray(bias_list) + np.asarray(variance_list),
        mode="lines",
        name="Total",
    )
)
fig.update_layout(
    title="Bias-Variance decomposition", xaxis_title="Penalty", yaxis_title="Value"
)
fig.show()

▶  Example 30.1 線形クラスにおけるBias and Variance

True functionを

$$f(x)=\theta {\theta }^{T}x$$

と定義し，$\theta \theta$ は未知であるとします．Training dataset $D$ を用いて，

$${\hat{f}}_D(x)={\widehat{\theta \theta }}_D^{T}x$$

を学習したとします．このとき，未知のデータセット $SS=\{(x^{\ast },y^{\ast })\}$ において，

$$\begin{array}{rl}\text{Bias}({\hat{f}}_D)& =\mathbb{E}[{\hat{f}}_D(x)-f(x)|SS]\\ & =\mathbb{E}[{\widehat{\theta \theta }}_D^{T}x-\theta {\theta }^{T}x|SS]\\ & =\mathbb{E}[{\widehat{\theta \theta }}_D-\theta \theta |SS{]}^{T}x\\ & =\text{Bias}({\widehat{\theta \theta }}_D{)}^{T}x\end{array}$$

同様に

$$\begin{array}{rl}& \mathrm{Var}({\hat{f}}_D)\\ & =\mathbb{E}[({\hat{f}}_D(x)-f(x){)}^2|SS]\\ & =\mathbb{E}[({\widehat{\theta \theta }}_D^{T}x-\theta {\theta }^{T}x{)}^2|SS]\\ & =\mathbb{E}[(({\widehat{\theta \theta }}_D^T-\theta {\theta }^T)x{)}^T(({\widehat{\theta \theta }}_D^T-\theta {\theta }^T)x)|SS]\\ & =\mathbb{E}[x^T({\widehat{\theta \theta }}_D-\theta \theta )({\widehat{\theta \theta }}_D-\theta \theta {)}^{T}x|SS]\\ & =x^T\mathrm{Var}(\widehat{\theta \theta })x\end{array}$$

なお，最後の式変形は $a$ を定数としたときに $\mathrm{Var}(a-X)=\mathrm{Var}(X)$ を用いています．

🍵 Green Tea Break

予測におけるBias-Variance分解は Overfitting，Underfitting という結びつけて考えることができます．

▶  Overfitting

training dataにおける予測精度はとても良いがtest dataにおける予測精度がとても悪い予測モデルに対してOverfittingという評価をします．これは予測モデルがtraining dataのノイズを拾いすぎてしまったときに発生します．

Overfittingは表現力(capacity)の高いモデルを用いたときに発生する可能性が高いです．実際に発生したときは，

• high bias and high variance
• low bias and high variance

いずれも考えることができますが，基本的にはhigh varianceを疑いアクションをとるという対処方針となります，対処方法として，

• 正則化
• より大きいデータセットの活用
• inputeとして用いる特徴量の数を少なくする
• より表現力が狭められたモデルを学習モデルに用いる
• BoostingよりBagging

といったことをトライすることが推奨されます．

▶  Underfitting

Underfittingとは，training dataにおいても十分な学習ができていない状況を指す概念です． モデルの表現力が低い場合に起きることが多く，またhigh biasの状況と関連しているケースが多いです．この場合，traing dataset自体にfitできていないので，さらなるデータ収集することは費用対効果に見合いません．

• 罰則項の緩和
• 特徴量の拡充
• より表現力のあるモデルのトライ
• Boosting

が推奨されます．

▶  Rule of thumbs

汎化性能が悪い場合に，まずbias-varianceのどちらの要因が重きを占めているのか？を考えることが重要です．その判別方法として経験則的に知られているステップとして，

• Training data errorが大きい場合は，もしモデルがトレーニングデータ自体にうまくフィットできないということになるので，そのモデルはhigh biasを持っている可能性が高いと判断する
• Cross validation errorとprediction erro in training datasetの差は、モデルや推定量の分散として扱うことができるので，Cross validation errorが高い場合は overfitting の状態と判断する

Training dataset errorに比べCross validation errorが高いときに，線形モデルからneural networkモデルへ変更することはあまり意味のない（むしろ状況を悪化させる）となるケースが多いというアクションになります．

Decomposition of Expected Test Error

上では未知のデータセット $\{x^{\ast },y^{\ast }\}$ について $\{x^{\ast }\}$ コントロールしたときのBias-Variance Decompositionを確認しましたが 未知のデータセット $(xx,y)$ とtraining dataset $D$ についてのexptected decompositonを確認していきます．

なお，Expected Prediction functionについて以下のnotationを導入します

$$\overline{f}=E_D\left[{\hat{f}}_D\right]=\underset{D}{\int }{\hat{f}}_{D}Pr(D)\partial D$$

$\overline{f}$ は $D$ の確率で重みづけたそれぞれの予測関数のweighted averageと理解できます．

$$\begin{array}{}\text{(30.2)}& \begin{array}{rl}& {\mathbb{E}}_{xx,y,D}[({\hat{f}}_D(xx)-y{)}^2]\\ & ={\mathbb{E}}_{xx,y,D}[({\hat{f}}_D(xx)-\overline{f}(xx)+\overline{f}(xx)-y{)}^2]\\ & =\underset{\text{Variance}}{\underbrace{{\mathbb{E}}_{xx,D}[({\hat{f}}_D(xx)-\overline{f}(xx){)}^2]}}+{\mathbb{E}}_{xx,y}[(\overline{f}(xx)-y{)}^2]\end{array}\end{array}$$

なお途中の式変形において次の性質を利用しています

$$\begin{array}{rl}& {\mathbb{E}}_{\mathbf{x},y,D}\left[({\hat{f}}_D(\mathbf{x})-\overline{f}(\mathbf{x}))(\overline{f}(\mathbf{x})-y)\right]\\ & ={\mathbb{E}}_{\mathbf{x},y}\left[E_D[{\hat{f}}_D(\mathbf{x})-\overline{f}(\mathbf{x})](\overline{f}(\mathbf{x})-y)\right]\because \text{LIE}\\ & ={\mathbb{E}}_{\mathbf{x},y}\left[(E_D[{\hat{f}}_D(\mathbf{x})]-\overline{f}(\mathbf{x}))(\overline{f}(\mathbf{x})-y)\right]\\ & ={\mathbb{E}}_{\mathbf{x},y}\left[(\overline{f}(\mathbf{x})-\overline{f}(\mathbf{x}))(\overline{f}(\mathbf{x})-y)\right]\\ & ={\mathbb{E}}_{\mathbf{x},y}\left[0\right]\\ & =0\end{array}$$

次に，Equation 30.2 の第二項について考えます．

$$\begin{array}{rl}& {\mathbb{E}}_{xx,y}[(\overline{f}(xx)-y{)}^2]\\ & ={\mathbb{E}}_{xx,y}[(\overline{f}(xx)-h(xx)+h(xx)-y{)}^2]\\ & ={\mathbb{E}}_{xx}[(\overline{f}(xx)-h(xx){)}^2]+{\mathbb{E}}_{xx,y}[(h(xx)-y{)}^2]\\ & ={\mathbb{E}}_{xx}[{\text{Bias}}^2(xx)]+{\mathbb{E}}_{xx}[\mathbb{E}[{ϵ}^2|xx]]\\ & =\mathbb{E}[{\text{Bias}}^2]+\mathbb{E}[{ϵ}^2]\end{array}$$

以上より

$${\mathbb{E}}_{xx,y,D}[({\hat{f}}_D(xx)-y{)}^2]=\mathbb{E}[\text{Variance of }{\hat{f}}_D]+\mathbb{E}[{\text{Bias}}^2]+\mathbb{E}[{ϵ}^2]$$