ベイズの基礎：課題16-4 夫は妻より年上か？

課題16-3 コインが表の確率分布 からの続き。


Activity 16-4: Marriage ages

課題：大阪市と京都市に在住の夫婦を無作為に24組抽出したところ「同い年：夫が年上：妻が年上 = 2 : 16 : 6」となった。夫が年上である割合を推定せよ。


単純な発想として

　16 / (16+6) ≈ 0.727

はどうだろうか？

答えは「No」ですね。理由は、同い年を除いて（年齢が違う夫婦のみで）考えているから。

では次はどうか？

　16 / 24 ≈ 0.667

これって「単にサンプルの平均」ですよね？


膨大なモデル：ありがとう R

今回、これまでの課題と格段に異なるのは、モデルの数が 101 個という多さ。夫が年上である確率（予測モデル）を 0% から 0.01 刻みの 101 個を想定する。

p = 0.00, 0.01, 0.02, ..., 0.98, 0.99, 1

各モデルの事前確率は「一様」として 1/101 。

尤度は、予測モデルが 0.2 の場合、「16 組が夫年上、6 組が妻年上」のデータから以下のように求める。

尤度 = 0.2¹⁶ × (1 - 0.2)⁶

事前確率と尤度から、全予測モデルの事後確率を算出して、グラフ化した。

p = seq(0,1,by = 0.01) # 0.01 刻みの確率モデル
plot(seq(0, 1, by=0.01),
  (p^16*(1-p)^6*(1/101))/sum(p^16*(1-p)^6*(1/101)),
  type="h", ylab = "Prob.", xlab="p Model",
  main="The Proportion of Couples in which the Husband is Older")

初めは OpenOffice のスプレッドシートで算出しようとしたが、10 分程度格闘して断念。R を使ったら 3 分未満で完了。モデル数が 1,000 件, 10万件だろうと大丈夫。つくづく、スプレッドシートで統計学をやっている人を気の毒に思う。そして Thanks, R!!


確率の推定

post = (p^16*(1-p)^6*(1/101))/sum(p^16*(1-p)^6*(1/101)) #事後確率
l <- list(p=seq(0,1,by=0.01),post=post) # {モデル, 事後確率} リスト作成

以下のように「モデル0」から「モデル1.0」の順に事後確率が格納される。

> l
$p
  [1] 0.00 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10 0.11 0.12
 [14] 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.20 0.21 0.22 0.23 0.24 0.25
 [27] 0.26 ...
 [92] 0.91 0.92 0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.00

$post
  [1] 0.0e+00 1.6e-28 1.0e-23 6.2e-21 5.8e-19 1.9e-17 3.3e-16 3.7e-15
  [9] 2.9e-14 1.8e-13 9.1e-13 3.9e-12 1.5e-11 5.0e-11 1.5e-10 4.3e-10
 [17] 1.1e-09 ...
 [97] 3.7e-05 7.7e-06 7.9e-07 1.5e-08 0.0e+00

「夫が年上の確率が 50% より高い確率」は「0.50 < 予測モデル ≦ 1.00」の事後確率の合計値

> sum(l$post[52:101])
[1] 0.9805412

つまり、98% の確率で「夫が年上の確率が 50% より高い」。

事後確率が高い順に並べ替える

> p2 <- l$p[sort.list(l$post, decreasing = T)] 
> post2 <- l$post[sort.list(l$post, decreasing = T)]
> l2 <- list(p=p2, post=post2)

> cumsum(l2$post) #事後確率の累積過程

  [1] 0.043 0.086 0.129 0.172 0.214 0.255 0.296 0.336 0.375 0.414 0.450

 [12] 0.487 0.521 0.555 0.587 0.618 0.647 0.675 0.702 0.727 0.751 0.773

 [23] 0.794 0.814 0.832 0.850 0.865 0.879 0.893 0.905 0.916 0.926 0.935

 [34] 0.944 0.951 0.957 0.963 0.968 0.973 0.977 0.980 0.984 0.986 0.989

 [45] 0.991 0.992 0.994 0.995 0.996 0.997 0.997 0.998 0.998 0.999 0.999

 ...

[100] 1.000 1.000

50% を超えるのが「0.521」の 13 番目、90% を超えるのが「0.905」の 30 番目。この二つの地点を、予測モデルで表示したのが以下。

> l2$p

$p

  [1] 0.73 0.72 0.74 0.71 0.75 0.70 0.76 0.69 0.77 0.68 0.67 0.78 0.66

 [14] 0.79 0.65 0.80 0.64 0.81 0.63 0.62 0.82 0.61 0.83 0.60 0.84 0.59

 [27] 0.58 0.85 0.57 0.56 0.86 0.55 0.87 0.54 0.53 0.88 0.52 0.51 0.89

 [92] 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0.00 1.00

テキストに沿って推定結果を表すと

You have just constructed a probability interval. The proportion p lies between 0.66 and 0.78 with probability 0.521 - that is, [0.66, 0.79] is a 52.1% probability interval.

予測モデル p は 0.66 と 0.78 の区間に 0.521 の確率で納まる

or

[0.66, 0.78] は52.1%確率区間

課題16-5 超能力あんのか？ に続く。