ベイズ統計なるものに関する質問です

中が見えない箱に玉が２０個入れてあるとします。そのうち５つを取り出したところ、全てが赤い玉でした。このとき２０個の玉の全てが赤い玉である確率はどうなるのでしょうか。
答えを分かり易く誘導記述していただくと、数式にそれほど強くない私としては助かります。

質問者からの補足コメント

• このような問題でも解があるのかな？と思っての質問でした。ベイズ統計についても、何も知らないまま、そんな問題でも答えが出ると何かの本で読んだことを思い出し、専門性の高い人に質問してみようと投稿した次第です。
  お二方の論戦は分からぬなりに大変勉強になりました。ご両名様に厚く感謝いたしながら質問を閉じさせていただきます。。

    補足日時：2019/02/22 07:54

No.22 回答者： ありものがたり 回答日時： 2019/02/21 13:54

←No.19お礼へ

うーん、そのミスを指摘できるのはスゲェとしか。
計算してたとき、6188 と 6586 が近すぎることに
微妙に違和感があったような気はするのだけれど、
検算する気が起こらないほど計算疲れしてた。
やっぱ、ゴリ押しが過ぎたかなあ。
御指摘に従って、真面目に計算してみよう。

使えそうなネタは aＣb = (a-1)Ｃ(b-1) + (a-1)Ｃb
くらいしか思いつかない。これを a=k+1, b=m+1 で使って
kＣm = (k+1)Ｃ(m+1) - kＣ(m+1).
これを k = m+1…n で総和すると、
Σ[k=m…n]kＣm - mＣm = (n+1)Ｃ(m+1) - (m+1)Ｃ(m+1) より
Σ[k=m…n]kＣm = (n+1)Ｃ(m+1).
これは、たぶん周知の式なのだろうが、私は知らなかった。

p(k) = (kＣm)/(nＣm) より、
求めたい確率は {1/(n+1)}p(n)/Σ[k=m…n]{1/(n+1)}p(k)
= (nＣm)/{(n+1)Ｃ(m+1)} = (m+1)/(n+1).
(たぶん鏡の中の)御友達の正解どおりになった。

No.21 回答者： kamiyasiro 回答日時： 2019/02/21 07:01

うんちくの多い#17です。

この問題の類似問題で、企業での実問題を挙げておきます。ご参考まで。

『今日、この工程では200台出荷する予定である。従来、この工程は、平均して0.5%の不良率であるが、正確に0.5%ではなく、その標準偏差は約0.1である。
ところが、今日は設備の調子が悪い。これまで20台生産して（非復元抽出）、5台も廃却不良が出た。
良品200台確保するため、正味何台生産しなければならないか予測して、追加部品を調達せよ。母確率の95％信頼上限を使って安全側で予測せよ。』

こんな問題を自動で解くようにして、生産管理システムのソフト中に入れておけば良いのです。

No.20 回答者： kamiyasiro 回答日時： 2019/02/21 04:37

うんちくの多い#17です。

それもで、かなり食い違いますね。

事前分布を一様分布として、5個非復元試行で観測し、残り15個が赤玉である確率Pは、私はベイズの計算で「一回の試行当たりの赤玉が出る確率p」を求めており、15個は二項分布を想定して解いています。

p＝0.8324
P＝0.06382357

貴殿の29％と私の６％、かなり食い違いが出ています。
何が違うのでしょうか。

No.19 回答者： ありものがたり 回答日時： 2019/02/21 03:26

値が随分食い違うなと思ったら、箱の球は20個じゃないか。

No.13 を訂正。

事前確率分布として、箱の中の赤球が
0, 1, 2, …, 20 個である確率を
各 1/21 づつの等確率と置く。

箱の中の赤球が k 個という条件下に
5 個取り出した球が全て赤である確率 p(k) は、
k < 4 のとき p(k) = 0,
k ≧ 5 のとき p(k) = (kＣ5)/(20Ｃ5).
表にすると以下のようになる。
p(5) = 1/(20Ｃ5),
p(6) = 6/(20Ｃ5),
p(7) = 21/(20Ｃ5),
p(8) = 56/(20Ｃ5),
p(9) = 126/(20Ｃ5),
p(10) = 252/(20Ｃ5),
p(11) = 462/(20Ｃ5),
p(12) = 792/(20Ｃ5),
p(13) = 1287/(20Ｃ5),
p(14) = 2002/(20Ｃ5),
p(15) = 3003/(20Ｃ5),
p(16) = 4368/(20Ｃ5),
p(17) = 6188/(20Ｃ5),
p(18) = 6586/(20Ｃ5),
p(19) = 11628/(20Ｃ5),
p(20) = 15504/(20Ｃ5).

5 個取り出した球が全て赤だったという条件下に
箱の中の赤球が 20 個である確率は、
(1/21)p(20)/Σ[k=0..20](1/21)p(k)
= 15504/(1+6+21+56+126+252+462+792+1287+2002+3003+4368+6188+6586+11628+15504)
= 15504/52282
≒ 0.2965…

電卓だけでできたが、大変だった。
何かプログラムでも組めばよかったかな？

この回答へのお礼

P(18)の分母ですが　８５６８　を　６５８６　と書き間違われています。　
それで　修正された確率は　１５５０４／５４２６４＝０．２８５７･･･
この答えは（５＋１）／（２０＋１）と同じです。
私の知り合いに数学の得意な人がいて、その人によると、ｎ個の玉のうちｍ個連続で赤が出た場合の、玉全てが赤の確率は　（ｍ＋１）／（ｎ＋１）だよと教えてくれました。
どうしてそうなるのかは、｢お前に説明したところで分からんだろう｣　と言われたので、尤もなことだと結果だけを覚えて帰った次第です。

お礼日時：2019/02/21 07:36

No.18 回答者： kamiyasiro 回答日時： 2019/02/19 06:44

うんちくの多い#17です。

ご質問者様、皆様、お待たせしました。
まず最初に#17の回答を破棄させて下さい。数値が若干異なりました。
また、今回はdhyper関数の使い方が分かったので私も勉強になりました。ありがとうございました。

さて今回は、事前分布の確率密度関数を、
・簡易的な離散確率ではなくキチンと連続分布と考えて、
・本来は積分でやるところをMCMC（マルコフ連鎖モンテカルロ法）を用いて、
正確に計算しました。これに伴って事前分布が一様分布の場合は1/21（離散）ではなくB(1,1)（連続）で与えるようにしました。

MCMCソフトはフリーソフトであるOpenBUGSを使いました。そのスクリプトは以下の通りです。結果は私の危惧通り、簡易的に刻んで計算した場合に対し、試行１回あたりの赤玉出現確率の期待値が次のように違っていました。

・一様分布を想定した場合
0.8928571 → 0.8324
・B(5,5)を想定した場合
0.6889171 → 0.6454

この結果、残り15個が全て赤玉である確率も変わりました。簡易計算はこのような違いが出るので、私は計算手順だけ示して結果を隠していましたが、やっと胸のつかえが取れました。また、他者の回答にあった「単なる条件付き確率の計算なのだった。」は間違いであるとお分かり頂けたと思います。しっかり理解して欲しいものです。

・OpenBUGS（WinBUGS）のスクリプト
・シミュレーション回数は10万回、バーイン区間の除去は1000回です。
・出力には、期待値（mean）だけでなく、信頼区間（パーセンタイル）も表示されています。

#_Red_ball
#
model
{
psi~dbeta(5,5)
x~dhyper(x,m,N,psi)
}
list(x=5,m=5,N=20)
#


①事前分布B(1,1)を設定した時のMCMCの結果（OpenBUGSの出力）

_____mean____sd______MC_error__val2.5pc_median__val97.5pc_start__sample
psi__0.8324__0.1386__6.608E-4__0.4906___0.8677__0.9948____1______100000

上記期待値（mean）を用いて残り15個が赤である確率を求めると（Rの出力）、

dbinom(15,15,0.8324)
[1] 0.06382357

②事前分布B(5,5)を設定した時のMCMCの結果（OpenBUGSの出力）

_____mean____sd______MC_error__val2.5pc_median__val97.5pc_start__sample
psi__0.6454__0.1214__3.794E-4__0.3949___0.6519__0.8609____1______100000

上記期待値（mean）を用いて残り15個が赤である確率を求めると（Rの出力）、

dbinom(15,15,0.6454)
[1] 0.001404218

ちなみに、企業では類似の技術問題が発生しますが、そのときは、こんな方法（ベイズ）で発生確率を見積もっています。技術者全員が理解してるわけでなく、データサイエンティストという人たちが計算を行います。

No.17 回答者： kamiyasiro 回答日時： 2019/02/18 06:16

#15です。

ご質問者様へ、しばらく閉じないで下さい。
OpenBUGSを使ってMCMCで解いた回答を載せたいと思います。プログラムを書く時間を下さい。ギブアップするかもしれません。

#14様へ、先に投稿したRのスクリプトに書いた関数を利用すると、

事前確率として離散一様分布1/21を想定した場合は、
> N <- 20
> j <- 5
> prob.func(j,N,1)
[1] 0.8928571・・・５個連続して赤が出た時に想定される母確率ｐの期待値
> dbinom(N-5,N-5,prob.func(j,N,1))
[1] 0.1826963・・・残り15個が赤である事象の確率

事前確率としてB(5,5)を想定した場合は、
> prob.func(j,N,2)
[1] 0.6889171・・・５個連続して赤が出た時に想定される母確率ｐの期待値
> dbinom(N-5,N-5,prob.func(j,N,2))
[1] 0.003736838・・・残り15個が赤である事象の確率

となります。グラフでは概算しか読めませんが、数値としてはこんな感じです。なお、企業では期待値を用いず安全側の95％信頼上限を使ったりします。
でも、この計算は積分ではなく刻んで求めていますので、正しいとは言えません。時間はかかりますが、MCMCでの結果を報告させて下さい。

No.16 回答者： ありものがたり 回答日時： 2019/02/18 00:32

答えは出さないの？

No.15 回答者： kamiyasiro 回答日時： 2019/02/17 22:53

#12です。

ご質問者を無視して、長々と細かな話を続けてスミマセンでした。ただ、計算方法は#13とはちょっと違いますので、ご注意ください。ベイズの計算で事後確率の期待値として母確率pを求めたら、20個とも赤である事象の確率Pを別途求めなければなりません。

No.14 回答者： ありものがたり 回答日時： 2019/02/17 18:14

＞お二方の専門性が高く高度な統計議論に

いえいえ、私が書いているのは、素朴な算数です。
少なくとも、専門用語を並べてはいないと思います。

No.13 回答者： ありものがたり 回答日時： 2019/02/17 17:59

ああ、そういうことか。

薀蓄が長いので、話の要点を見落としていた。
なるほど、これは私のほうの言いがかり
だったらしい。

事前確率分布として、箱の中の赤球が
0, 1, 2, …, 10 個である確率を
各 1/11 づつの等確率と置く。

箱の中の赤球が k 個という条件下に
5 個取り出した球が全て赤である確率 p(k) は、
k < 4 のとき p(k) = 0,
k ≧ 5 のとき p(k) = (kＣ5)/(10Ｃ5).
表にすると以下のようになる。
p(5) = 1/252,
p(6) = 6/252,
p(7) = 21/252,
p(8) = 56/252,
p(9) = 126/252,
p(10) = 252/252.

5 個取り出した球が全て赤だったという条件下に
箱の中の赤球が 10 個である確率は、
(1/11)p(10)/Σ[k=0..10](1/11)p(k)
= 252/(1 + 6 + 21 + 56 + 126 + 252)
= 6/11.

ベイズとは言っても、事前確率→事後確率の
順方向なので、単なる条件付き確率の計算なのだった。

この回答へのお礼

お二方の専門性が高く高度な統計議論に（理解が及ばないながら）感心いたしました。

お礼日時：2019/02/17 18:10