消しゴム事件

  • 事件の概要
    • こんな事件の相談を受けた
    • 相談はこちらの「話題1」
    • 全校の男子生徒数N人の学校で、全男子生徒にk種類の消しゴムが配られた
    • k種類の消しゴムの個数の内訳はn1,n2,...,nkという
    • その後、その消しゴムの落し物があった
    • 「●ちゃん(女子)♥」と書かれていた
    • 見つかった場所は音楽室
    • その消しゴムは、ある特定の日に落とされたことが確実で、その日に音楽室に出入りした男子はm人
    • 落し物の消しゴムの種類が「タイプi」であるとの情報が入った
    • ■ちゃん(女子)は意中の男子△の消しゴムの種類が何なのかの情報を入手することにした
    • 次の2通りが考えられる
      • 「タイプi」だった
      • 「タイプj(j\ne i)」だった
    • さて、△は「●ちゃん♥」なのか、そうでないのか!という問題
  • 考えよう
  • 「タイプi」だった場合
    • 音楽室に出入りして、消しゴムを落としたかもしれない男子m人のうち、タイプiの消しゴムを持っていたのは、t=1,2,...,t_{MAX}=\text{min}(n_i,m)人のいずれか
      • 1人は落としているから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(j \ne i)」だった場合
    • △が落とし主である確率は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]))
}