ベイジアン鶴亀算/t分布モデル

鶴亀算を解こう。

ある庭に何羽かの鶴と何匹かの亀がいる。
20人の観測者が頭の数と足の数を数えたところ次のようなデータが得られた。

observed_heads observed_feet

96 275

96 272

94 274

95 272

...中略... ...中略...

96 273

庭にいる鶴の数と亀の数を求めよ。

観測者によって頭の数と足の数がバラバラだ。これは各観測者が数え間違いをしている (かもしれない) からだ。

ベイズの考えを使って解く

典型的な鶴亀算の問題は、与えられた「頭の数」と「足の数」を真の値として決定論的に計算する。
今回の問題では「頭の数」と「足の数」に観測誤差を含むので、誤差を加味して「鶴の数」と「亀の数」を評価しなければいけない。

このような問題にはベイズ推定の考え方が使える。隠れたパラメータがあり、パラメータに応じて確率的にデータが観測されていると考える。ベイズの定理 (逆確率の定理) を使うと、データから逆算してパラメータを推定できる。

いま「鶴の数」「亀の数」は未知だが、庭にいる「鶴の数」「亀の数」に応じて「頭の数」「足の数」が観測されているはずだ。ベイズ推定で、観測された「頭の数」「足の数」から「鶴の数」「亀の数」を推定できる。

実行環境・ツール

今回も Stan を使って推定する。Stan は Python から呼び出すことにする。Python の実行環境には Google Colab を使う。
というわけで完全に動作する Google Colab Notebook がここにある。

colab.research.google.com

モデル定義

確率モデルを定義しよう。
ここでの確率モデルとは、「頭の数」「足の数」がどのような分布から生成されているかを仮定するものだと考えればいい。

真の鶴の数 T* と真の亀の数 K* に対して、真の頭の数 H*, 真の足の数 F* は、

H* = T* + K*
F* = 2 T* + 4 K*

こうなるはずだ。

よって観測された頭の数 H と観測された足の数 F は、それぞれ T* + K* と 2 T* + 4 K* を中心とした山なりの分布から生成されていると考えたい。
山なりの分布なので中心である T* + K* や 2 T* + 4 K* の値が観測されやすい。しかし確率分布からの標本なので、もちろん中心から外れてしまうこともある。

今回は H, F が t分布に従うと仮定する。
t分布は3つのパラメータ (母数) によって形が決まる。 $\nu$ , $\mu$ , $\sigma$ だ。

$H \sim StudentT(h | \nu = nu, \mu = T + K, \sigma = sigma\_H)$
$F \sim StudentT(h | \nu = nu, \mu = 2T + 4K, \sigma = sigma\_F)$

このような仮定をおいて T, K, sigma_H, sigma_F, nu を推定対象としよう。

Stan コード

モデルを Stan のコードで書き下すとこうなる。
このコードでは H = T + K, F = 2 T + 4 K という知識を反映しているものの。T = ~, K = ~ の形に解いた式は入っていない。
H, F の値が尤もらしくなるように推定すれば、自然に T, K の尤もらしい値も求まるはずだ。

data {
  int<lower=1> N;      // データ数
  vector[N] H;         // 観測された頭の数
  vector[N] F;         // 観測された足の数
  real H_bar;          // 頭の数の標本平均
}
parameters {
  real<lower=0> T;         // 鶴の数（未知）
  real<lower=0> K;         // 亀の数（未知）
  real<lower=0> sigma_H;   // 観測誤差(観測者の能力): 頭の数
  real<lower=0> sigma_F;   // 観測誤差(観測者の能力): 足の数
  real<lower=2> nu;        // t自由度
}
transformed parameters {
  vector[N] mu_H = rep_vector(T + K, N);        // E[H] = T + K
  vector[N] mu_F = rep_vector(2*T + 4*K, N);    // E[F] = 2T + 4K
}
model {
  (T + K) ~ normal(H_bar, 10);    // 総数は観測された頭の数に近いはず
  sigma_H ~ student_t(3, 0, 10);
  sigma_F ~ student_t(3, 0, 10);
  nu ~ gamma(2, 0.1);             // 平均≈20のロバスト化

  H ~ student_t(nu, mu_H, sigma_H);  // 頭の数 H は mu_H := T + K を中心とした分布に従う
  F ~ student_t(nu, mu_F, sigma_F);  // 足の数 F は mu_F := 2*T + 4*K を中心とした分布に従う
}

データ

冒頭の問題文では "中略" となっていたデータの全貌をお見せしよう。

observed_heads	observed_feet
96	275
96	272
94	274
95	272
97	273
102	266
98	267
100	277
94	268
97	273
95	270
97	272
97	268
96	271
94	265
98	281
95	272
92	275
98	274
96	273

全部で20件ある。

サンプリングと推定

というわけで事後分布を求めて各パラメータを推定した。

	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
T	56.630	54.294	58.854	1.159	6817.0	1.0
K	39.659	38.261	41.033	0.704	6813.0	1.0
sigma_H	2.221	1.423	3.145	0.446	9629.0	1.0
sigma_F	3.836	2.488	5.357	0.756	9835.0	1.0
nu	20.302	2.333	47.689	13.800	9759.0	1.0

鶴の数 T は、平均が 56.630、おおよその範囲が 55\~59 となった。鶴の数 K は、平均が 39.659、おおよその範囲が 38\~41 となった。

「T ≒ 56.630 → 鶴は57羽」のように点推定するだけでなく、確信度に応じて幅を持たせて推定できるのがベイズ推定のうれしいところだ。

また sigma_H, sigma_F の値にも着目してほしい。この値は H, F を観測したときの標準偏差 (ばらつき) を表している。
つまり sigma_H, sigma_F が大きいほど「ばらつきが大きい」=「観測者の能力が低い」と解釈できる。

データを増やして再推定する

追加で20件のデータを取ったとして、推定結果がどう変わるか見てみよう。
追加のデータはこれだ。

observed_heads	observed_feet
95	279
95	273
95	267
94	274
95	271
90	277
95	274
97	277
94	273
95	268
92	270
96	269
94	270
93	277
100	275
95	265
96	270
94	277
97	272
96	272

最初の20件と合わせて合計40件のデータを、先ほどと同じモデルに渡して再度推定してみる。

	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
T	55.071	53.526	56.579	0.776	6995.0	1.0
K	40.513	39.576	41.430	0.474	7003.0	1.0
sigma_H	2.097	1.554	2.723	0.299	9471.0	1.0
sigma_F	3.699	2.810	4.664	0.481	10348.0	1.0
nu	20.076	2.697	45.919	13.147	9186.0	1.0

T, F の平均値はさほど変わっていないが、幅を持たせて推定したときの幅が狭くなっているのに着目してほしい。これはデータを増やしたことで、推定結果により確信が持てたということだ。

	データ20件	データ40件
T の平均	56.630	55.071
T の範囲	54\~59	54\~57
F の平均	39.659	40.513
F の範囲	38\~41	40\~41

答え合わせ

今回使ったデータは私が適当に生成したものなので、実は真の値を知っている。
真の値は T=56, K=40 だ。

	真の値	点推定	区間推定
T	56	55.071	54\~57
F	40	40.513	40\~41

まぁ、いい感じ...？

まとめ

ベイズの考え方を使って、方程式を解くことなく鶴亀算に取り組んでみた。

私からは以上です。

無駄と文化

実用的ブログ