遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

割り付けの場合を分けてその確率を足し合わせる

R 割り付け
  • k人居て名札を配るときに、「曖昧な記憶に基づいて」配るとしよう
  • 記憶の程度がごく対称的な(異常に対称的な)場合を考えてみる(計算が簡単だから)
  • aさんにaさんの名札を渡す確率はx、aさんにaさん以外の名札を渡す確率は1-xで、aさん以外のだれの名札を渡すかは、均等とする
  • aさんの場合もbさんの場合もすべて同じとする
  • このとき、aさんにだけ着目してaさんが自分の名札を受け取るのはどれくらいの確率になるだろうか
  • 次のようなことに気を付ける
    • その確率はaではない
    • わりつけの場合すべてk!についてx,(1-x)/(k-1)を使って計算して足し合わせると1にならない
  • 色々考えても仕方ないので、全部の通りを地道に計算する、というのも悪くない(ただし、人数が増えると大変だ)
  • ３人くらいならまあよいだろう
  • ただし、せっかく、「すべての場合」を計算するのだから、『記憶の程度』は自由にしてもよいことにする(印象の強かった人は正しい名札をもらう確率が高いし、印象の弱かった人は確率が低いだろう)
library(gtools)

probAllPermutation<-function(P,logarithm=TRUE){
	n<-nrow(P)
	if(!logarithm){
		P<-log(P)
	}
	perms<-permutations(n,n,1:n)
	np<-length(perms[,1])
	ret<-rep(0,np)
	for(i in 1:np){
		ret[i]<-sum(diag(P[perms[i,],]))
	}
	list(perms=perms,probs=ret)
}

CalcProbHx<-function(P,logarithm=TRUE){
	n<-nrow(P)
	p.out<-probAllPermutation(P,logarithm=logarithm)
	M<-matrix(0,n,n)
	for(i in 1:n){
		for(j in 1:n){
			selected<-which(p.out$perms[,i]==j)
			M[i,j]<-sum(exp(p.out$probs[selected]))
		}
	}
	list(perms=p.out$perms,probs=p.out$probs,prob.mat=t(M))
}

n<-3
x<-0.4
y<-0.3
xx<-x^(1/3)
yy<-(y/xx)^(1/2)
P<-matrix(c(xx,yy,yy,yy,xx,yy,yy,yy,xx),byrow=TRUE,n,n)
p.out<-probAllPermutation(P,logarithm=FALSE)
#p.out

calc.out<-CalcProbHx(P,logarithm=FALSE)
calc.out
exp(calc.out$probs)

M1<-(calc.out$prob.mat)/apply(calc.out$prob.mat,1,sum)

M2<-t(calc.out$prob.mat)/apply(calc.out$prob.mat,2,sum)


M1
M2
apply(M1,1,sum)
apply(M1,2,sum)
apply(M2,1,sum)
apply(M2,2,sum)
  • さて、正答に関して「対称性」を仮定して問題を簡単にすることにする
  • その代わり、人数が増えても大丈夫なようにしたい
  • 今日の記事の冒頭で割り付けでの、「正解人数別の場合の数」の計算はしてあるので、それを使う
  • 着目している人について正解する場合については、残りの人に関して、正解人数が0,1,...,(k-1)を考慮すればよい
  • 着目している人について不正解する場合は、次のようにする。残りの人に関して、正解人数が0,1,...,(k-1)を考慮するのは上と同じだが、ただし、着目者で不正解なので、不正解の名札の本人は、正しい名札をもらうことはありえなくなっている。したがって、k-1人のうちi人が「正解」とは言っても、そのi人に「必ず不正解」の人を入れたつもりの時は、「本当の正解の人数は１、減らす」必要があるし、そうでないときは、本当の正解の人数はiのまま。ここで、「正解として『必ず不正解』の人を入れるか入れないか」は、i/(k-1)対(k-1-i)/(k-1)である。場合の数を考えているから、「正解確率x」の多寡は関係ない。これなら計算が簡単だ
# 返り値のp1は着目者を正解する確率、p2は不正解のそれ
# ただし、足して１にならないことに注意
Seigo<-function(x,k){
	y<-(1-x)/(k-1)
	w<-Waritsuke3(k-1)
	P1<-P2<-0
	for(i in 0:(k-1)){
		P1<-P1+w[k,i+1]*y^i*x^(k-1-i)
		P2<-P2+w[k,i+1]*(y^i*x^(k-1-i)*i/(k-1) + y^(i+1)*x^(k-1-i-1)*(k-1-i)/(k-1) )
	}
	p1<-x*P1
	p2<-y*P2
	list(p1=p1,p2=p2)
}
Seigo2<-function(x,y,k){
	#y<-(1-x)/(k-1)
	w<-Waritsuke3(k-1)
	P1<-P2<-0
	for(i in 0:(k-1)){
		P1<-P1+w[k,i+1]*y^i*x^(k-1-i)
		P2<-P2+w[k,i+1]*(y^i*x^(k-1-i)*i/(k-1) + y^(i+1)*x^(k-1-i-1)*(k-1-i)/(k-1) )
	}
	p1<-x*P1
	p2<-y*P2
	list(p1=p1,p2=p2)
}
  • k=3の場合について、地道法と比べてみる
Seigo(x,k)


xx<-x
yy<-y
P<-matrix(c(xx,yy,yy,yy,xx,yy,yy,yy,xx),byrow=TRUE,n,n)
p.out<-probAllPermutation(P,logarithm=FALSE)
#p.out

calc.out<-CalcProbHx(P,logarithm=FALSE)
calc.out

> Seigo(x,k)
$p1
[1] 0.73125

$p2
[1] 0.002375

> 
> 
> xx<-x
> yy<-y
> P<-matrix(c(xx,yy,yy,yy,xx,yy,yy,yy,xx),byrow=TRUE,n,n)
> p.out<-probAllPermutation(P,logarithm=FALSE)
> #p.out
> 
> calc.out<-CalcProbHx(P,logarithm=FALSE)
> calc.out
$perms
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    1    3    2
[3,]    2    1    3
[4,]    2    3    1
[5,]    3    1    2
[6,]    3    2    1

$probs
[1] -0.3160815 -6.0968251 -6.0968251 -8.9871968 -8.9871968 -6.0968251

$prob.mat
         [,1]     [,2]     [,3]
[1,] 0.731250 0.002375 0.002375
[2,] 0.002375 0.731250 0.002375
[3,] 0.002375 0.002375 0.731250
  • 超対称な場合、「ある個人」について正解する確率と不正解する確率はxとkのみの関数になるので、どうなっているか見てみることができる
ks<-2:10
ts<-seq(from=0.01,to=1-0.01,by=0.01)

P1s<-P2s<-matrix(0,length(ks),length(ts))
for(i in 1:length(ks)){
	for(j in 1:length(ts)){
		tmp<-Seigo(ts[j],ks[i])
		P1s[i,j]<-tmp$p1
		P2s[i,j]<-tmp$p2
	}
}

image(P1s/P2s)
 matplot(t(P1s/P2s),type="l")
ks<-2:10
xy.ratios<-1.2^(10:-10)

xy.ratios<-seq(from=1,to=10,length=100)
preV<-postV<-P1s<-P2s<-P2multis<-Xvalue<-matrix(0,length(ks),length(xy.ratios))
for(i in 1:length(ks)){
	for(j in 1:length(xy.ratios)){
		rr<-xy.ratios[j]
		Xvalue[i,j]<-1/(1+rr*(ks[i]-1))
		tmp<-Seigo(Xvalue[i,j],ks[i])
		P1s[i,j]<-tmp$p1
		P2s[i,j]<-tmp$p2
		P2multis[i,j]<-P2s[i,j]*(ks[i]-1)
		preV[i,j]<-Xvalue[i,j]/(Xvalue[i,j]*(1+(ks[i]-1)*rr))
		postV[i,j]<-P1s[i,j]/(P1s[i,j]+(ks[i]-1)*P2s[i,j])
	}
}

image(P1s/P2s)
 matplot(1/t(P1s/P2s),type="l")
 matplot(1/t(P1s/P2multis),type="l")
 
 matplot((P1s/P2s),type="l")
 matplot((P1s),type="l")
 matplot((P2s),type="l")
matplot(t(P1s),type="l")
matplot(t(P2s),type="l")

plot(c(Xvalue),c(P1s/P2s))