最近ベイズがいろんなところで使われていると聞くので少しいじってみようかと。
勝手にこんな状況が起こったらこんなことが使えるんじゃないかと仮定しながら勉強してみます。
生産・製造の現場では機械の状況や出力値の設定によって不良品が出たりすることもあると思います。
その時、なんだか夕方になると不良品が増える気がする。
不良が増える時の値は
(52.27 55.52 47.52 53.05 54.18)
でした。
ひょっとしたら母パラメータが変わっているのかもしれません。
推定してみましょう。
生産にとって重要な機械の出力部分をセンサーで測っていたとします。
購入したてで、朝から晩まで問題なく動いていて、不良品も出していなかった時のパラメータを引っ張り出してきたとします。
既知である正常な時のパラメーター
平均=50
標準偏差=2
今もおそらくこのパラメーターで動いているだろうと仮定します。
（パラメーターがわかっているのがパラメトリックベイズ
パラメーターが不明の時、それらしい分布を当てはめて行うのがノンパラメトリックベイズだそうです。）
ベイズの展開式を今回つかいます。
(これならわかる！ベイズ統計学 ナツメ社)
左式は事後確率です。データDが、原因Hのもとで得られた確率です。
つまり、データDが起こった後、調べて得られた原因Hの確率です。
右式の分子部分はP(D|H)尤度とP(H)事前確率です。
尤度は原因HのもとでデータDが得られるもっともらしい確率です。
事前確率はデータがかかわる前の分析される前の確率です。
分母はそれぞれ起こる確率の総和です。周辺確率と呼ばれています。
とりあえず想定しているパラメーターで分布を書いてみます。
正規分布に従うものとします。とりあえず正規分布が書きたいときはcurve関数を使いますが、今回は事前と事後で比較するのでdnormを使います。
さっそくRで書いてみます。
hist(rnorm(1000,50,2))
1000個のデータを仮定してみると範囲は44から56くらい出てきましたので、
推定の時に考える範囲は40から60くらいを見ておけばいいのではないでしょうか。
ということで
para_mean <- 50
para_sd <- 2
x_axis <- seq(40,60,0.1)
prob_dens <- dnorm(x_axis, mean=para_mean, sd=para_sd)
plot(x_axis, prob_dens, ylim=c(0,0.5), type="l")
既知のパラメーターから想定される確率分布はこのような形になるんですね。
この確率分布が事前確率だと解釈します。
不良時のパラメーターを考えていきます。
(52.27 55.52 47.52 53.05 54.18)
この対象に対して推定を行ってもデータ数が少なすぎてまともな結果にはならなさそうです。
ちなみに
mean(para_defe)
# [1] 52.508
さきほどと同じように不良時のパラメーターについて分布を確かめてみます。
para_defe <- c(52.27, 55.52, 47.52, 53.05, 54.18)
plot(dnorm(para_defe, mean=x_axis, sd=1),type="p")
それぞれのデータについての確率分布を書くことができました。
ちなみに標準偏差については元の1を使用しています。
もし標準偏差が求まる程度に不良品のデータがたくさんあるならば仮定でなくても大丈夫です。
（その場合推定は必要ない気もしますが・・・）
何パターンか試してみてそれっぽいものを見つける必要があるんでしょうか？
こういう時、実際にはどんな値を入れているのか気になってます。
それぞれ別々の確率でしたが、求めるのは不良を出すときの母パラメータでしたので、これらのデータの確率を一つの確率にくっつけるためかけ合わせます。
pro_1 <- dnorm(para_defe[1], mean=x_axis, sd=1)
pro_2 <- dnorm(para_defe[2], mean=x_axis, sd=1)
pro_3 <- dnorm(para_defe[3], mean=x_axis, sd=1)
pro_4 <- dnorm(para_defe[4], mean=x_axis, sd=1)
pro_5 <- dnorm(para_defe[5], mean=x_axis, sd=1)
joint=pro_1*pro_2*pro_3*pro_4*pro_5
plot(x_axis,joint,type="l",col="red")
これが今回の尤度にあたる部分です。
やっていることは同時確率を求めているんです。
同時確率も尤度も数式的には一緒になります。
事前確率と尤度が求まりましたのであとは分子と分母の形にします。
numerator <- prob_dens * joint
denominator <- sum(numerator)
posterior <- numerator / denominator
plot(x_axis, prob_dens, ylim=c(0,0.5), type="l")
points(x_axis, posterior,type="l", col="blue")
青色が事後確率ということになりました。
母パラメーターを求めてみます。
which.max(posterior)
# [1] 125
x_axis[125]
# [1] 52.4
んー
単純に平均した時よりも微妙にかわっていますが、この程度の差は必要なんでしょうか・・・
もしくはここから修正を加えるのか、不良データ発生時の標準偏差を1と仮定したのがまずかったのか・・・
ということで、不良時の母平均がもとまりました。
不良時のデータは母平均からどのくらい離れたところにあるのでしょうか。
緑の線3シグマからはみ出しています。。
見当違い？
ちなみに今回のデータは
rnorm(5,50.5,3)
でなんとなく発生させた値でしたが、ピッタリあてることができませんでした。
標準偏差を3としたときは
3シグマの範囲に収まりました。