正味の重量分布を予測する（後半：ベイズ推定）
はじめに
チョコボールの重量分布を予測します。
前回の記事では、 チョコボールの重量分布を予測するための方法として最尤推定を行いました。
しかし、最尤推定には課題がありました。 最尤推定は点推定であるため、推定した値がどのくらい確信を持っているのかを 評価できないのでした。
そこで今回は、 ベイズ推定を使ってパラメータの予測結果の分布を推定します。
何を予測するのか
チョコボールの正味の重量分布を予測します。 モデルは、前回最尤推定した場合と同じで、重量は正規分布すると仮定し、 正規分布の平均と分散を推定します。
パラメータ（平均と分散）をとすると、
データ
が与えられた条件での事後分布
を推定します。
これまでの計測で以下のような計測データが得られています。
ベイズ推定によってこのデータからどのように分布が推定されるのか、
次章以降で計算してみます。
ベイズ推定の基本
ベイズの定理は以下の式で定義されます。
ここで、分母のは、データXが得られる確率です。
は、
が
について周辺化されたものだと解釈できます。
ここで、は可能な
の集合を表します。
この式からわかるように、
は全体を積分すると１にするための
正規化定数です。
ということで、ベイズの式は次のように表せられます。
この式は、
データが与えられた条件での
の確率（事後分布
）は、
尤度関数
と事前分布
の積に比例するということをあわらしています。
事前分布というのは、
「分析者が、求めたいパラメータ
がどのような分布をしていると思っているのか？」を表現するものです。
主観確率とも呼ばれますが、まさに分析者の主観です。
主観ってなんぞと思うでしょうが、
問題によっては予め事前知識を持っていることが少なくなく、
その事前知識を盛り込んでデータの分析をしようということです。
例えば今回の例では、チョコボールの仕様上、重量は28gということになっています。
ということは、仮に何もデータが無くても28gからそう遠くない重量になっているだろうと思えるわけです。
あたりがつけられないような場合には、一様分布を当てはめれば良いのです（無情報事前分布）。
あとは、実際に得られたデータを当てはめて式を解けばよいのですが、 ざっくりいうと、素直にベイズの式を展開するには積分がいっぱい出てきて大変です。 大変というか解析的に計算できない場合もあります。
そこで、計算のテクニックとして２つの方法があります。
- 事前分布と尤度関数の組み合わせを特別な組み合わせにすることで、比較的容易に計算できる事後分布にする（共役事前分布）
- 複雑な関数をコンピュータシミュレーションで近似計算する（MCMC法）
今回は、後者のMCMC法を利用して計算してみます。
MCMCによる重量分布の予測
MCMC法というのは、すごーくざっくり説明すると、 複雑な形状の確率分布の大きさに依存してデータの生成をシミュレーションします（下図）。
このシミュレーションでサンプルしたデータの密度で確率分布の形状を表現して いろいろな統計計算に利用するというものです。
事前準備
pythonでpymcを使います。 また、Anacondaを使います。
pymcのインストールは以下のコマンドを入力するだけです。
conda install pymc
コード
jupyter notebookで作ってみました。
５番目のセルでモデルの設定をしています。
# 精度$\tau$(分散の逆数)
tau = 1.0 / pm.Uniform("std", 0, 100)**2
# 分布の中心
## 事前分布の中心として、製品仕様の28を、tauを0.01(分散=100)
center = pm.Normal("center", 28, 0.01)
# 観測とモデルを結びつける
observations = pm.Normal("obs", center, tau, value=data['net_weight'].values, observed=True)
# Modelオブジェクト生成
model = pm.Model([tau, center])
求めたいパラメータは、重量分布の平均と分散です。 分散の変数として"tau"を、平均の変数を"center"として、 それぞれ一様分布と正規分布を初期分布として設定しています。
また、観測モデルをobservationという名前で定義しています。 観測は平均と分散を中心にして正規分布で得られると設定しています。 観測値は得られるものなので、observed=Trueとしています。
最後に、pymcのModelを作成しています。
その次のセルでMCMCでのパラメータ推定を実行しています。 今回は、収束性を見るためにも、バーンイン期間無しでシミュレーションしています。
mcmc = pm.MCMC(model) mcmc.sample(30000) # バーンイン無し
MCMCシミュレーションでは、最初の数千点は探索をするために大きくブレるため、 最終的な分布には含めません。 この最初の期間をバーンインと言って、捨てた方が良いデータです。
以下のようにすると、バーンイン期間10000点を結果に含めないという意味になります。
mcmc.sample(30000, 10000)
以下の図は、今回のMCMCシミュレーションの軌跡をプロットした図です。
上段が平均の軌跡で、下段が標準偏差の軌跡です。 左側が全ての軌跡をプロットしたもので、右側がバーンイン期間10000点を除いた軌跡です。 平均も標準偏差も早々に収束していることが見えます。
推定結果
30回分の計測データを使ってベイズ推定による平均と分散の推定をしたところ、予測分布は以下のようになりました。 下図の左が平均の予測分布で右が標準偏差の予測分布です。
この図から、重量分布の平均は29.5付近にあるだろうと強く確信していることを示しています。 平均の分布の緑の線は90%信用区間を表しており、 約29.39〜約29.57の範囲に収まることを90%確信しています。 仕様よりも1g以上多いと予測しています。
それぞれの分布の期待値を平均と分散の推定値として分布を描いてみます。
紫色の線がベイズ推定による予測分布で、 青線が最尤推定による予測分布です。 結果としてはどちらもほぼ同じ分布を予測できています。
データ数毎の推定を見てみる
ベイズ推定の良い点は、データが少なく、予測があやふやなときには 予測があやふやであるということが分布の広がりから言える点にありました。 そこで、データを計測していく毎にパラメータの推定分布がどのように変化していったのか見てみます。
この図は、左上からデータ数が[1, 3, 5, 10, 20, 30]点での平均と標準偏差の推定分布を表しています。 それぞれの推定分布は左側が平均の分布で、右側が標準偏差の分布です。 平均の分布はx軸のレンジを揃えています。 また、平均の分布には90%信用区間を緑の線で表しています。
この図から、データが増える毎に予測の確信度が上昇していることが見て取れます。
補足
今回の設定は、尤度関数と事前分布には正規分布を使いました。 この場合、共役事前分布の関係になっているので、事後分布も正規分布になります。 そのため、わざわざMCMCで繰り返しシミュレーションしなくても、 解析的に計算することができます。
ですが、MCMCを使いたかったので今回はMCMCを使いました。 共役事前分布を使って解析的に求める方法については別の機会に書くかもしれません。
終わりに
ということで、ベイズ推定を使ってチョコボールの重量予測をやってみました。 最終的な結果としては最尤推定と同じように、平均29.5gくらいの重量と予測できました。
良心的ですね！森永さん！
参考文献
とりあえず細かい理論は置いておいても使ってみたいという時に。
- 出版社/メーカー: 森永製菓
- 発売日: 2016/03/01
- メディア: 食品&飲料
- この商品を含むブログを見る
SIMERST 改良版 携帯タイプはかり ポケットデジタルスケール(秤) 0.001g-100g精密 業務用(プロ用) デジタルスケール 電子天秤
- 出版社/メーカー: Simerst
- メディア:
- この商品を含むブログを見る
