2012-02-04

消しゴム事件

尤度推定

事件の概要
- こんな事件の相談を受けた
- 相談はこちらの「話題１」
- 全校の男子生徒数N人の学校で、全男子生徒にk種類の消しゴムが配られた
- k種類の消しゴムの個数の内訳はn1,n2,...,nkという
- その後、その消しゴムの落し物があった
- 「●ちゃん（女子）♥」と書かれていた
- 見つかった場所は音楽室
- その消しゴムは、ある特定の日に落とされたことが確実で、その日に音楽室に出入りした男子はm人
- 落し物の消しゴムの種類が「タイプi」であるとの情報が入った
- ■ちゃん（女子）は意中の男子△の消しゴムの種類が何なのかの情報を入手することにした
- 次の２通りが考えられる
  - 「タイプi」だった
  - 「タイプj( $j\ne i$ )」だった
- さて、△は「●ちゃん♥」なのか、そうでないのか！という問題
考えよう
「タイプi」だった場合
- 音楽室に出入りして、消しゴムを落としたかもしれない男子m人のうち、タイプiの消しゴムを持っていたのは、t=1,2,...,人のいずれか
  - １人は落としているから0人の可能性はない
- t=tである確率は $Pr(t)=\;_m C_t p^t (1-p)^{m-t};p=\frac{n_i}{\sum n_s}$
- 「正確な確率」は $PrEx(t) = \; \frac{n_{1.}!n_{2.}!n_{.1}!n_{.2}!}{n_{..}!t!(n_{1.}-t)!(n_{.1}-t)!(n_{..}-n_{1.}-n_{.1}+t)!}$
- 実際、これは、 $t=t$ である確率だが、「タイプi」の落し物をしやすいtとそうでないtとがあって、それを考慮すると、尤度は $Pr(t),PrEx(t)\times t$ になる
- さて、tの値別に、△が落とし主である確率は
  - t=1の場合、△が落とし主である確率は1
  - t=2の場合、△が落とし主である確率は1/2
  - t=tの場合、△が落とし主である確率は1/t
- したがって $\frac{\sum_{t=1}^{t_{MAX}} {Pr(t)/t}}{\sum_{t=1}^{t_{MAX}} Pr(t)}$
「タイプj()」だった場合
- △が落とし主である確率は0

mySuspect<-function(p,N){
	ks<-0:N
	if(p>=1){
		ks<-N
		Prk.log<-lgamma(N+1)-lgamma(ks+1)-lgamma(N-ks+1)+ks*log(p)+log(ks)
		PrHanawa.log<-Prk.log-log(ks)
		LikeHanawa<-0
		LikeAll<-0
		Pr.Hanawa<-0
		LR.Hanawa<-0
	}else{
		Prk.log<-lgamma(N+1)-lgamma(ks+1)-lgamma(N-ks+1)+ks*log(p)+(N-ks)*log(1-p)+log(ks)
		PrHanawa.log<-Prk.log-log(ks)
		LikeHanawa<-sum(exp(PrHanawa.log[2:length(ks)]))
		LikeAll<-sum(exp(Prk.log)[2:length(ks)])
		Pr.Hanawa<-LikeHanawa/LikeAll
		LR.Hanawa<-LikeHanawa/(LikeAll-LikeHanawa)
	}
	
	return(list(pr.log.per.k=Prk.log,pr.log.suspect.per.k=PrHanawa.log,like.all=LikeAll,like.suspect=LikeHanawa,pr.suspect=Pr.Hanawa,LR.suspect=LR.Hanawa,p=p,N=N,ks=ks))
}
mySuspect.Exact<-function(m1,m2){
	N<-sum(m1)
	ks<-0:min(m1[1],m2[1])
	M<-cbind(ks,m1[1]-ks,m2[1]-ks,m1[2]-m2[1]+ks)
	minM<-apply(M,1,min)
	selected<-which(minM>=0)
	M<-M[selected,]
	selected2<-which(M[,1]>0)
	M<-M[selected2,]
	tmp<-sum(lgamma(m1+1),lgamma(m2+1))-lgamma(N+1)
	Prk.log<--apply(lgamma(M+1),1,sum)+tmp+log(M[,1])
	PrHanawa.log<-Prk.log-log(M[,1])
	LikeHanawa<-sum(exp(PrHanawa.log))
	LikeAll<-sum(exp(Prk.log))
	Pr.Hanawa<-LikeHanawa/LikeAll
	LR.Hanawa<-LikeHanawa/(LikeAll-LikeHanawa)

	
	return(list(pr.log.per.k=Prk.log,pr.log.suspect.per.k=PrHanawa.log,like.all=LikeAll,like.suspect=LikeHanawa,pr.suspect=Pr.Hanawa,LR.suspect=LR.Hanawa,p=m1[1]/N,N=N,ks=M[,1]))
}