「Probabilistic Programming and Bayesian Methods for Hackers」を読んだ - (5章)

第5章は損失関数について。損失関数の計算に真の値を使わずにどう計算するのか気になっていたが事後分布を使うと知って納得。

損失関数

損失関数は予測した分布がどれくらい悪いのかを評価する関数。計算に便利で導入されたらしいが、コンピュータがあるので今使われる理由は「robust to application」だかららしい。よく使われるような有名なものを以下に書いた。

$$ L(\theta, \hat{\theta}) = \begin{cases} &(\theta - \hat{\theta})^{2} &: &\text{平方損失}\\ &|\theta - \hat{\theta}| &: &\text{絶対損失}\\ &1_{\hat{\theta} \neq \theta} &: &\text{0-1損失}\\ &-\hat{\theta}log(\theta) - (1-\hat{\theta})log(1-\theta) &: &\text{log損失 } (\theta, \hat{\theta} \in [0,1]) \end{cases} $$ 平方損失は線形回帰に使われる。0/1損失は分類アルゴリズム、log損失は期外学習に使われるらしい。

この章の例もそうだがこれら以外で自分で簡単に定義できるので他にも考えることができて、、などとすることもできる。

頻度主義の場合はのに推定値を代入して使うと思われる。ベイズ主義ではは予測された事後分布になり、を計算して代用する。

例1 「The Price is Right」

• 「The Price is Right」というテレビ番組は2人の出演者が商品の値段を予測比べをして、勝てばその商品がもらえる番組らしい。ただし予測が実際の値段を超えず、250以上下回っていない場合のみはいけない。

• 事後分布の計算

この例では価格の事前分布としてこれまでの番組での価格分布を正規分布にして用いる。そしてここからサンプリングして出てきた結果を元に事後分布を計算するのがいつもの流れだが、今回はまだ価格を知ってはおらず出演者が予測した価格分布しか情報がない。そこで事前分布の平均と分散の逆数を、として価格がxになる確率をとして事後分布を計算する。

ここで分布を修正するときに@pymc.potentialを使う。@pymc.potentialは出来上がったモデルに特に親子関係がない条件を付け加えたいときに使われるデコレータで、potentialになった変数はpymc.modelに付け加える必要がない。

@pm.potential
def error(true_price, ...):

などとするといい。

黒い曲線が元の分布で赤色の分布が事後分布。

• 損失関数の計算

import scipy.optimize as sop

def showdown_loss(guess, true_price, risk=80000):
    loss = np.zeros_like(true_price)
    loss[guess > true_price] = risk
    ix = [(guess <= true_price) * (guess > true_price - 250)]
    loss[ix] = -2 * np.abs(true_price[ix])
    ix = [guess <= true_price - 250]
    loss[ix] = np.abs(true_price - guess - 250)

expected_loss = lambda guess, risk: shutdown_loss(guess, mcmc.trace("true_price")[:], risk)
for _p in risks:
...
    min_results = sop.fmin(expected_loss, 15000, args=(_p,), disp=False)

scipy.optimize.fminは関数の最小値を探してくれる。

損失関数が最小なのは事後分布の最大値の20000前後からかなり小さくなっている。

例2: 金利の推測

$$ R = \alpha + \beta x + \epsilon $$ というモデルについて考える。 とは平均0、精度0.0001の正規分布、には一様分布なstdの2乗の逆数を使う。

import pymc as pm
from pymc.Matplot import plot as mcplot

std = pm.Uniform("std", 0, 100, trace=False)  # this needs to be explained.


@pm.deterministic
def prec(U=std):
    return 1.0 / (U) ** 2

beta = pm.Normal("beta", 0, 0.0001)
alpha = pm.Normal("alpha", 0, 0.0001)


@pm.deterministic
def mean(X=X, alpha=alpha, beta=beta):
    return alpha + beta * X

obs = pm.Normal("obs", mean, prec, value=Y, observed=True)
mcmc = pm.MCMC([obs, beta, alpha, std, prec])

mcmc.sample(100000, 80000)
mcplot(mcmc)

収束していたらOK。

損失関数は飽きたのでパス。

例3: ダークマターによるハロー探し

TODO