2015年10月17日土曜日

ベイズの基礎:正規分布による推定

正規分布は最初にあらず からの続き。

課題:平日の大学生の睡眠時間を推定する。ただし標準偏差を 0.5 時間とする。観測データは {6.5 時間} の 一つを用いる。

「平日の大学生の睡眠時間」は概ね正規分布に従うと仮定できるだろう。したがって、推定するのは正規分布で、標準偏差が分かっているので、求めるのは平均値。

 注意:現実的には、標準偏差は不明、課題用に便宜的に提供されている。


左のように3つのモデル(6, 7, 8 時間睡眠)に同じ事前確率を割り当てる。

尤度は正規分布関数より

LIKELIHOOD = e−z2/2


観測データが1件

左の正規分布関数から、以下のように事後確率 post を求める

> M <- c(6,7,8)
> z <- (6.5 - M)/0.5
> like <- exp(-z^2/2)
> prod <- like*(1/3)
> post <- prod/sum(prod)
[1] 0.495462643 0.495462643 0.009074715

事後確率の算出の流れは左の通り。

ほぼ6時間と7時間で「半々の確率」という結果。これではダメだと容易に分かる。観測データが一件しかないため、こういう結果になった。


モデルが不十分

観測データが 28 件で、同様に事後確率を求める。

x <- c(8.0,6.0, 6.5, 7.0, 7.0, 6.5, 6.5, 6.0,7.0,7.0,6.0,6.0,7.5,6.0,7.0,7.0,6.5,6.0,6.5,6.5,9.0,9.0,7.0,5.0,8.0,6.5,6.0,6.0)

z は左の式で求める。「(横棒上付)x」は、観測値の平均。

> M <- c(6,7,8)
> z2 <-(sqrt(length(x))*(mean(x)-M))/0.5
> like2 <- exp(-z2^2/2)
> prod2 <- like2*(1/3)
> post2 <- prod2/sum(prod2)
> post2
[1] 6.914400e-13 1.000000e+00 3.305701e-37

事後確率の算出の流れは左の通り。

7時間が「100%」という更にダメダメな結果になりました。


モデルを細分化

モデルを「6.0 から 7.5 まで 0.1時間刻み」に変更。

M2 <- seq(6,7.5,by=0.1)
z3 <- (sqrt(length(x))*(mean(x)-M2))/0.5
like3 <- exp(-z3^2/2)
prod3 <- like3*(1/3)
post3 <- prod3/sum(prod3)
plot(M2, post3, type="h", xlab="M", ylab="Prob", main="Post for Mean Sleeping Time")

6.66.9 時間の確率:
> post3[7:10]
[1] 0.1197590 0.3670438 0.3670438 0.1197590
> sum(post3[7:10])
[1] 0.9736055

97% 、左の正規分布の中央から4本の合計値。

課題18-5,10 正規分布なの? に続く。

0 件のコメント:

コメントを投稿