消しゴム事件 続く

• 一件落着に思えた消しゴム事件、さらなる相談が入った
  • 消しゴムに「●ちゃん♥」と書いてあったのは本当らしい
  • 消しゴムのタイプというのは、音楽の先生からの情報らしいが、音楽の先生は消しゴムの種類について、ちょっとあやふや(キャラクター商品なのだが、キャラクターを識別できていないようだ)という
  • また、△の消しゴムタイプに関する情報も、消しゴムを配った先生からのもので、この先生のキャラクター識別力にも問題がありそうだ、という
• 考えよう
  • まず、消しゴムのタイプが本当は1,2,...,kであるので、それぞれの場合について考える
    • 本当はタイプiであるときに消しゴムを配った先生が"j=1,2,...,kである"と教えてくれる確率がある
    • 本当はタイプiであるときに音楽の先生が"j=1,2,...,kである"と判定結果を示す確率がある
  • 消しゴム配布先生が、音楽の先生がと教えてくれたとき、本当のタイプがである尤度は
    • である
    • したがって、本当のタイプがであるのはと考えるのがよいだろう
  • 消しゴム配布先生が、音楽の先生がと教えてくれたとき、本当のタイプがiである尤度がである
    • このとき、その本当のタイプi=1,2,...,kのそれぞれについて、の場合には、「△が落とし主だ」と判断する可能性がある
    • の場合は、その判断は適切(ただし、確率的に)
    • の場合には、たまたま、「配布先生と音楽の先生との発言が一致した」だけなので、このような場合には、「△が落とし主だ」と判断するのは不適切である
  • 「配布先生と音楽の先生」の言葉は本当は信用ならない可能性があるけれども、「鵜呑み」にする場合
    • のすべてでが「△が持ち主である」という判定を許すだけの尤度をもたらすとすると
    • の確率で、「△が持ち主である」という判定を下すことになる
    • しかしながらの分は「適切」だが
    • の分は「不適切」である
    • 誤判断の混入が問題となる
  • 「配布先生と音楽の先生」の言葉は本当は信用ならない可能性があるのでその分を割り引く場合
    • 適切であるのみで、「△が持ち主である」という判定を許すだけの尤度がもたらされるとすると、結局、「ご判断」する点では同じである
    • の高さによって、エラー情報率が高いときには、が十分に下がるようにして、「誤判断」が起きないようにすることが求められる

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)
		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)
		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))
}

SimpleLR<-function(q,p,N){
	q/(q+p*N)
}

# 消しゴムのタイプ数
N.type<-3

# typeの実数
Np<-c(5,100,150)
# typeの比率
p<-Np/sum(Np)

# 落し物消しゴムの真のタイプ
o.type<-1:N.type
# 花輪君消しゴムの真のタイプ
h.type<-1:N.type

# タイプコールはどうなるか
# 正答率
P1<-0.9
P2<-0.9

Keshigomuya<-diag(rep(P1,N.type))
Ongaku<-diag(rep(P2,N.type))

for(i in 1:N.type){
	Keshigomuya[i,which(Keshigomuya[i,]==0)]<-(1-sum(Keshigomuya[i,]))/length(which(Keshigomuya[i,]==0))
	Ongaku[i,which(Ongaku[i,]==0)]<-(1-sum(Ongaku[i,]))/length(which(Ongaku[i,]==0))
}

Keshigomuya
Ongaku

# 
# 消しゴムのタイプ数
N.type<-2

# typeの実数
Np<-c(10,100000)
# typeの比率
p<-Np/sum(Np)

# 落し物消しゴムの真のタイプ
o.type<-1:N.type
# 花輪君消しゴムの真のタイプ
h.type<-1:N.type

# タイプコールはどうなるか
# 正答率
P1<-0.9
P2<-0.9

Keshigomuya<-diag(rep(P1,N.type))
Ongaku<-diag(rep(P2,N.type))

for(i in 1:N.type){
	Keshigomuya[i,which(Keshigomuya[i,]==0)]<-(1-sum(Keshigomuya[i,]))/length(which(Keshigomuya[i,]==0))
	Ongaku[i,which(Ongaku[i,]==0)]<-(1-sum(Ongaku[i,]))/length(which(Ongaku[i,]==0))
}

Keshigomuya
Ongaku




Accuracy<-function(Obs1,Obs2,thres,p,N){
	Ongaku<-Obs1
	Keshigomuya<-Obs2
	N.type<-length(Ongaku[1,])
	# 落し物消しゴムの真のタイプ
	o.type<-1:N.type
	# 花輪君消しゴムの真のタイプ
	h.type<-1:N.type


	True.ij<-list()
	Call.mat<-list()
	Stat.Naive.mat<-list()
	cnt<-1
	for(i in 1:length(o.type)){
		o.call<-Ongaku[i,]
		for(j in 1:length(h.type)){
			True.ij[[cnt]]<-c(i,j)
			h.call<-Keshigomuya[j,]
			Call.mat[[cnt]]<-outer(o.call,h.call,FUN="*")
			Stat.Naive.mat[[cnt]]<-matrix(0,N.type,N.type)
			for(k1 in 1:N.type){
				for(k2 in 1:N.type){
					if(k1==k2){
						Stat.Naive.mat[[cnt]][k1,k2]<-mySuspect(p[k1],N)$pr.suspect
						#Stat.Naive.mat[[cnt]][k1,k2]<-SimpleLR(1,p[k1],N)
					}
				}
			}
			cnt<-cnt+1
		}
	}
	
	TrueFalse<-list()
	cnt<-1
	for(i in 1:N.type){
		for(j in 1:N.type){
			TrueFalse[[cnt]]<-c(0,0)
			if(i==j){
				for(k1 in 1:N.type){
					for(k2 in 1:N.type){
						if(k1==k2){
							if(Stat.Naive.mat[[cnt]][k1,k2]>thres){
								if(i==k1){
									TrueFalse[[cnt]][1]<-TrueFalse[[cnt]][1]+Call.mat[[cnt]][k1,k2]*Stat.Naive.mat[[cnt]][k1,k2]
										TrueFalse[[cnt]][2]<-TrueFalse[[cnt]][2]+Call.mat[[cnt]][k1,k2]*(1-Stat.Naive.mat[[cnt]][k1,k2])
								}else{
									TrueFalse[[cnt]][2]<-TrueFalse[[cnt]][2]+Call.mat[[cnt]][k1,k2]
								}
							}
						}
						
					}
				}
			}else{
				
			}
			cnt<-cnt+1
		}
	}

	list(TrueFalse=TrueFalse,N.type=N.type,Obs1=Obs1,Obs2=Obs2,thres=thres,N=N,
	True.ij=True.ij,Call.mat=Call.mat,Stat.Naive.mat=Stat.Naive.mat)

}
# 犯人扱いにする閾値
thres<-0.99
# 犯人候補者数
N<-10

Accuracy(Obs1=Ongaku,Obs2=Keshigomuya,thres=thres,p=p,N)

# マーカー数
N.m<-5
p.cuts<-seq(from=0.9,to=1,by=0.01)
g<-0.05
N.type<-2
Niter<-1
for(xx in 1:Niter){
	# マーカーの信用度
	ps1<-rep(1,N.m)-runif(N.m)*0.1
	ps2<-rep(1,N.m)-runif(N.m)*0.1
	worseps<-ps1
	for(i in 1:N.m){
		if(ps1[i]>ps2[i])worseps[i]<-ps2[i]
	}
	for(yy in 1:length(p.cuts)){
		p.cut<-p.cuts[yy]
		selected<-which(worseps>p.cut)
		prod.p1<-prod(ps1[selected])
		prod.p2<-prod(ps2[selected])
		right.call<-prod.p1*prod.p2
		n.selected<-length(selected)
		g.f<-g^n.selected
		p<-c(g.f,1-g.f)

		Keshigomuya<-diag(rep(prod.p1,N.type))
		Ongaku<-diag(rep(prod.p2,N.type))

		for(i in 1:N.type){
			Keshigomuya[i,which(Keshigomuya[i,]==0)]<-(1-sum(Keshigomuya[i,]))/length(which(Keshigomuya[i,]==0))
			Ongaku[i,which(Ongaku[i,]==0)]<-(1-sum(Ongaku[i,]))/length(which(Ongaku[i,]==0))
		}
		if(n.selected>0){
			print(Accuracy(Obs1=Ongaku,Obs2=Keshigomuya,p=p,thres=thres,N))
		}
		
	}

}