決断のための分布合算 京大(医)統計遺伝学分野 山田 亮
小児～若年発症 ループス 血栓リスクが高そうに思うが… 大急ぎで電子カルテDBを検索した：かき集めても血栓例は10例 N Engl J Med 2011; 365:1758-1759November 10, 2011DOI: 10.1056/NEJMp1108726
わからなくても決断する あなたは冒険旅行中 分かれ道があって、電光掲示板がある 11例目のあなたは、どちらの道を選ぶか 『右の道を選んだ者、７名あり。４名は幸福に、３名は不幸になった』 『左の道を選んだ者、３名あり。２名は幸福に、１名は不幸になった』 11例目のあなたは、どちらの道を選ぶか 我を過ぐれば憂ひの都あり、 我を過ぐれば永遠の苦患あり、 我を過ぐれば滅亡の民あり
○ × 和 X 4 3 7 Y 2 1 ○率をベータ分布で推定 (4+1)/(7+2) 4/7 道X 道Y 期待値 最頻値 2/3 ○率 p <- seq(from=0,to=1,length=100) X <- c(4,3) Y <- c(2,1) dx <- dbeta(p,X[1]+1,X[2]+1) dy <- dbeta(p,Y[1]+1,Y[2]+1) Mean.x <- (X[1]+1)/(sum(X)+2) Mean.y <- (Y[1]+1)/(sum(Y)+2) Mode.x <- X[1]/sum(X) Mode.y <- Y[1]/sum(Y) Mean.x Mean.y Mode.x Mode.y matplot(p,cbind(dx,dy),type="l") abline(v=Mean.x,col=1) abline(v=Mean.y,col=2) abline(v=Mode.x,col=1,lty=2) abline(v=Mode.y,col=2,lty=2) 最頻値 2/3 ○率 (2+1)/(3+2)
期待値で選択することは「悪くない」 方針 「期待値」が大きい方を選ぶ 「期待値」が同じなら、どちらかを選ぶ
ゼロからのスタート Xの○率 0.2 Yの○率 0.25 n.iter <- 20 N <- 1000 res <- list() Suc.x <- Suc.y <- Select.X <- matrix(0,n.iter,N) ps <- c(0.6,0.65) for(i in 1:n.iter){ res[[i]] <- matrix(0,N,4) for(j in 2:N){ m1 <- (res[[i]][j-1,1]+1)/(res[[i]][j-1,1]+res[[i]][j-1,2]+1) m2 <- (res[[i]][j-1,3]+1)/(res[[i]][j-1,3]+res[[i]][j-1,4]+1) if(m1==m2){ s <- sample(1:2,1) }else if(m1>m2){ s <- 1 }else{ s <- 2 } res[[i]][j,] <- res[[i]][j-1,] x <- sample(1:0,1,prob=c(1-ps[s],ps[s])) tmp <- 2*(s-1)+x+1 res[[i]][j,tmp] <- res[[i]][j,tmp] +1 Suc.x[i,] <- res[[i]][,2]/(res[[i]][,1]+res[[i]][,2]) Suc.y[i,] <- res[[i]][,4]/(res[[i]][,3]+res[[i]][,4]) Select.X[i,] <- res[[i]][,1]+res[[i]][,2] matplot(2:N,t(Select.X[,2:N])/(2:N),type="l",xlab="No.samples",main="Fraction of X")
X Y X期待値 Ｙ期待値 ○ × 0.5 1 0.333333 2 0.4 3 4 0.285714 0.25 5 0.2 6 0.222222 7 0.272727 8 0.307692 9 0.375 10 0.352941 11
Y × Y ○ X × X ○
Y × Y ○ X × X ○
○率の高いYばかりが選ばれるようになった … and they lived happily ever after めでたし、めでたし Y × Y ○ X × X ○
○率の高いYばかりが選ばれるようになった … and they lived happily ever after めでたし、めでたし Y × Y ○ X × X ○
本当にYばかりが 選ばれるようになるか
Ｓｅｌｅ
Yに落ち着くか Xに落ち着くか Ｙ Ｘ
確率的な決断 Multi-armed bandit 問題 複数のスロットマシンがあって、それぞれのマシンには「当たり」の確率が決まっているが、その確率が不明であるという マシンを１つずつ選んでは、勝負をして、各マシンの当否結果を記録しながら、勝負を繰り返すことにする どんなルールで選ぶと、儲けが最大になりやすいか、という問題
確率的な決断 Multi-armed bandit 問題 その状況でのThomson samplingとかの方が良い結果が得られることが知られている。 ごく大雑把に言うと、 データを見ても、「100％、どのアームがよいとは言い切れない」から、データから見て、「得策らしくないアームも、ある程度(確率的に)は選ぼう」 いったん、悪い方を選び勝ちになっても、判断を修正するポテンシャルが「確率的な決断」によってもたらされる
『確率的な決断にする！』
わからなくても決断する 分かれ道があって、電光掲示板がある 11例目のあなたは、どうするか 『確率的な決断が大事である』 『従って、この分かれ道に奇数回目に来た者には、電光掲示板は点灯せず、偶数回目に来た者には、点灯することとする』 11例目のあなたは、どうするか １．適当に選ぶ ２．出直す
わからなくても決断する 分かれ道があって、電光掲示板がある 11例目のあなたは、どうするか 『確率的な決断が大事である』 『従って、この分かれ道に奇数回目に来た者には、電光掲示板は点灯せず、偶数回目に来た者には、点灯することとする』 11例目のあなたは、どうするか (１) 適当に選ぶ→こうすればある程度の「確率的決断」が起きる ２．出直す
わからなくても決断する 分かれ道があって、電光掲示板がある 11例目のあなたは、どうするか 『確率的な決断が大事である』 『従って、この分かれ道に奇数回目に来た者には、電光掲示板は点灯せず、偶数回目に来た者には、点灯することとする』 11例目のあなたは、どうするか (１) 適当に選ぶ→こうすればある程度の「確率的決断」が起きる (２) 出直す
何が問題か 偶数・奇数で不平等
何が問題か 偶数・奇数で不平等 平等な何かを作ってみる
何を比較する？ 「どちらの道を選ぶと○になる確率が高いのか」 「どちらの道が『○率が高い』のか」 これは○の期待値 「どちらの道の『○の期待値』が高いのか」 「どちらの道が『○率が高い』のか」
何を比較する？ 「どちらの道を選ぶと○になる確率が高いのか」 「どちらの道が『○率が高い』のか」 これは○の期待値 「どちらの道の『○の期待値』が高いのか」 「どちらの道が『○率が高い』のか」 違う設問だから答えも違う
何を比較する？ 「どちらの道を選ぶと○になる確率が高いのか」 「どちらの道が『○率が高い』のか」 これは○の期待値 「どちらの道の『○の期待値』が高いのか」 「どちらの道が『○率が高い』のか」 違う設問だから答えも違う Xの方が高いかもしれないし、Ｙの方が高いかもしれない
○ × 和 X 4 3 7 Y 2 1 ○率をベータ分布で推定 (4+1)/(7+2) 4/7 道X 道Y 期待値 最頻値 2/3 ○率 p <- seq(from=0,to=1,length=100) X <- c(4,3) Y <- c(2,1) dx <- dbeta(p,X[1]+1,X[2]+1) dy <- dbeta(p,Y[1]+1,Y[2]+1) Mean.x <- (X[1]+1)/(sum(X)+2) Mean.y <- (Y[1]+1)/(sum(Y)+2) Mode.x <- X[1]/sum(X) Mode.y <- Y[1]/sum(Y) Mean.x Mean.y Mode.x Mode.y matplot(p,cbind(dx,dy),type="l") abline(v=Mean.x,col=1) abline(v=Mean.y,col=2) abline(v=Mode.x,col=1,lty=2) abline(v=Mode.y,col=2,lty=2) 最頻値 2/3 ○率 (2+1)/(3+2)
○ × 和 X 4 3 7 Y 2 1 期待値 (4+1)/(7+2) 期待値(2+1)/(3+2)
「どちらの道の『○の期待値』が高いのか」 × 和 X 4 3 7 Y 2 1 期待値 (4+1)/(7+2) 期待値(2+1)/(3+2) 「どちらの道の『○の期待値』が高いのか」
○ × 和 X 4 3 7 Y 2 1 (4+1)/(7+2) 「どちらの道が『○率が高い』のか」
○ × 和 X 4 3 7 Y 2 1 (4+1)/(7+2) (2+1)/(3+2) 「どちらの道が『○率が高い』のか」
計算できる(式の導出は省略)
○率が高い確率応じて X,Yを「確率的に」選択してみよう
○率の高いYばかりが選ばれるようになった … and they lived happily ever after めでたし、めでたし Y × Y ○ X × X ○
『いいことを教えてやろう』
『右の道を選んだ者、７名あり。４名は幸福に、３名は不幸になった』 その内訳は 男 ５名。２名は幸福に、３名は不幸に 女 ２名。２名とも幸福に 『左の道を選んだ者、３名あり。２名は幸福に、１名は不幸になった』 男 ３名。２名は幸福に、１名は不幸に 女 は左の道を選んでおらぬ
男女合算 と 女のみ 女のみ 男女合算 女のみ 男女合算 p <- seq(from=0,to=1,length=100) 男女合算 と 女のみ 女のみ 男女合算 女のみ p <- seq(from=0,to=1,length=100) #y <- cbind(dbeta(p,4+1,3+1),dbeta(p,2+1,1+1),dbeta(p,2+1,3+1),dbeta(p,2+1,0+1),dbeta(p,2+1,1+1),dbeta(p,0+1,0+1)) y <- cbind(dbeta(p,4+1,3+1),dbeta(p,2+1,1+1),dbeta(p,2+1,0+1),dbeta(p,0+1,0+1)) y1 <- cbind(dbeta(p,4+1,3+1),dbeta(p,2+1,0+1)) y2 <- cbind(dbeta(p,2+1,1+1),dbeta(p,0+1,0+1)) par(mfcol=c(1,2)) matplot(p,y1,type="l",ylab="X",ylim=c(0,3)) matplot(p,y2,type="l",ylab="Y",ylim=c(0,3)) par(mfcol=c(1,1)) 男女合算
「男女に違いなし」なら 「男女に違いあり」なら 男女合算の情報を使った方が正確 男女合算の情報に基づいて集計した方が、早く、収束する 男女別々の情報を使った方が正確
0.58 0.25 0.42 0.75 女のみ 男女合算 p <- seq(from=0,to=1,length=100) #y <- cbind(dbeta(p,4+1,3+1),dbeta(p,2+1,1+1),dbeta(p,2+1,3+1),dbeta(p,2+1,0+1),dbeta(p,2+1,1+1),dbeta(p,0+1,0+1)) y <- cbind(dbeta(p,4+1,3+1),dbeta(p,2+1,1+1),dbeta(p,2+1,0+1),dbeta(p,0+1,0+1)) y1 <- cbind(dbeta(p,4+1,3+1),dbeta(p,2+1,0+1)) y2 <- cbind(dbeta(p,2+1,1+1),dbeta(p,0+1,0+1)) yall <- cbind(dbeta(p,4+1,3+1),dbeta(p,2+1,1+1)) yfemale <- cbind(dbeta(p,2+1,0+1),dbeta(p,0+1,0+1)) par(mfcol=c(1,2)) matplot(p,y1,type="l",ylab="X",ylim=c(0,3)) matplot(p,y2,type="l",ylab="Y",ylim=c(0,3)) par(mfcol=c(1,1)) dx <- dbeta(p,4+1,3+1) dy <- dbeta(p,2+1,1+1) dxy <- outer(dx,dy,"*") dx.female <- dbeta(p,2+1,0+1) dy.female <- dbeta(p,0+1,0+1) dxy.female <- outer(dx.female,dy.female,"*") par(mfcol=c(2,2)) matplot(p,yall,type="l",ylab="X",ylim=c(0,3)) image(dxy,xlim=c(0,1),ylim=c(0,1),xlab="X",ylab="Y") contour(dxy,add=TRUE) abline(0,1) matplot(p,yfemale,type="l",ylab="Y",ylim=c(0,3)) image(dxy.female,xlim=c(0,1),ylim=c(0,1),xlab="X",ylab="Y") contour(dxy.female,add=TRUE) Decision_beta.2(c(4,3,2,1)+1) Decision_beta.2(c(2,0,0,0)+1) > Decision_beta.2(c(4,3,2,1)+1) [1] 0.4242424 > Decision_beta.2(c(2,0,0,0)+1) [1] 0.75 0.42 0.75
男女合算 女のみ 道の選択確率が異なる 『道 Ｘ、道 Ｙ、どちらにしよう？』 『男女合算、女のみ、どちらにしよう？』 確率的に選んだ 男女合算 女のみ 道の選択確率が異なる 『道 Ｘ、道 Ｙ、どちらにしよう？』 確率的に選んだ 『男女合算、女のみ、どちらにしよう？』 確率的に選んでみる
男女合算 女のみ 道の選択確率が異なる 『道 Ｘ vs. 道 Ｙ、どちらにしよう？』 『男女合算 vs. 女のみ、どちらにしよう？』 男女合算 女のみ 道の選択確率が異なる 『道 Ｘ vs. 道 Ｙ、どちらにしよう？』 確率的に選んだ 『男女合算 vs. 女のみ、どちらにしよう？』 確率的に選んでみる
男女合算 道X,道Yのベータ分布 Y X X
男女合算→道X 道X→男,道Y→女のベータ分布
X ○ × 和 男 2 3 5 女 道X 男,女のベータ分布 女 X 男
X ○ × 和 男 2 3 5 女 道X 男,女のベータ分布 (p_m=p, p_f=p) 女 X 男
X ○ × 和 男 2 3 5 女 道X 男, 女のベータ分布 (p_m,p_f) (p_m=p, p_f=p) 女 X 男
(p_m,p_f) 『男女に違いがあってもよい』という立場 (p_m=p, p_f=p) 『男女に違いがない』という場合 女 X 男
X ○ × 和 男 2 3 5 女 女 X 男
X ○ × 和 男 2 3 5 女 女 X 男
X ○ × 和 男 2 3 5 女 女 この積分値を活用する X 男
『男女差がある場合』 サンプル数が０：積分値１ サンプル数が増える：小さくなる 『男女差がない場合』より 小さくなる程度は急速
小さくはなるが、対立仮説が作る同時分布のピークが対角線のごく近傍に来るので、小さくなっても困らない 『男女差がない場合』 サンプル数が０：積分値１ サンプル数が増える：小さくなる 小さくはなるが、対立仮説が作る同時分布のピークが対角線のごく近傍に来るので、小さくなっても困らない
X ○ × 和 男 2 3 5 女 どう活用する？ X
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差あり』 の同時分布 正方形部分の積分は１
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差あり』 の同時分布 正方形部分の積分は１
『男女に差あり』 の同時分布 正方形部分の積分は こちらも１ 『男女に差あり』 の同時分布 正方形部分の積分は１ X ○ × 和 男 女 X 女 X ○ × 和 男 2 3 5 女 『男女に差あり』 の同時分布 正方形部分の積分は こちらも１ 『男女に差あり』 の同時分布 正方形部分の積分は１
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差なし』 部分
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差なし』 部分
X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差なし』 部分
１ １ １
１ １ 積分 １
？ ？ 事前確率 事後確率 ？ ？
0.5 ？ 事前確率 事後確率 0.5 ？
１ １ 積分 １
0.5 1/(1+r) 事前確率 事後確率 0.5 r/(1+r) =r
仮説の比率が決まれば ２つのベータ分布の 重みづけ混合分布 1/(1+r) Beta(a+c+1,b+d+1) + r/(1+r) Beta(c+1,d+1)
道X、道Y、 それぞれの混合ベータ分布 に基づき 『〇率が高くなる確率』 による確率的選択 モンテカルロに 『○率が高くなる確率』を算出
やってみよう
女 男 200 800 1→1000人 8割男、2割女 男 X:0.2, Y:0.4 女 X:0.2, Y:0.4 0.2 vs. 0.25 10000人 女 n.pt <- 1000 X <- matrix(sample(1:2,n.pt,replace=TRUE,prob=c(0.8,0.2)),ncol=1) Prob.Vec <- list() tmp.Pr <- (X-1)*0.2+0.2 Prob.Vec[[1]] <- cbind(tmp.Pr,1-tmp.Pr) Prob.Vec[[2]] <-cbind(tmp.Pr,1-tmp.Pr) better.1 <- which(X==1) better.2 <- which(X==2) n.iter <- 10 out.mat.1 <- matrix(0,n.iter,length(better.1)) out.mat.2 <- matrix(0,n.iter,length(better.2)) for(i in 1:n.iter){ model.2.out <- my.simulate.model.2(model.2,X,Prob.Vec,n.iter=n.pt) better.1.selection <- (-1)*(model.2.out$selection[better.1]-2) better.2.selection <- model.2.out$selection[better.2]-1 out.mat.1[i,] <- cumsum(better.1.selection)/(1:length(better.1)) out.mat.2[i,] <- cumsum(better.2.selection)/(1:length(better.2)) } matplot(t(out.mat.1),type="l") matplot(t(out.mat.2),type="l") 男 200 800
1→1000人 8割男、2割女 男 X:0.2, Y:0.2 女 X:0.4, Y:0.4 女 n.pt <- 1000 X <- matrix(sample(1:2,n.pt,replace=TRUE,prob=c(0.8,0.2)),ncol=1) Prob.Vec <- list() tmp.Pr <- (X-1)*0.2+0.2 Prob.Vec[[1]] <- cbind(tmp.Pr,1-tmp.Pr) Prob.Vec[[2]] <-cbind(tmp.Pr,1-tmp.Pr) better.1 <- which(X==1) better.2 <- which(X==2) n.iter <- 10 out.mat.1 <- matrix(0,n.iter,length(better.1)) out.mat.2 <- matrix(0,n.iter,length(better.2)) for(i in 1:n.iter){ model.2.out <- my.simulate.model.2(model.2,X,Prob.Vec,n.iter=n.pt) better.1.selection <- (-1)*(model.2.out$selection[better.1]-2) better.2.selection <- model.2.out$selection[better.2]-1 out.mat.1[i,] <- cumsum(better.1.selection)/(1:length(better.1)) out.mat.2[i,] <- cumsum(better.2.selection)/(1:length(better.2)) } matplot(t(out.mat.1),type="l") matplot(t(out.mat.2),type="l") 男 800 200
女 男 800 200 1→1000人 8割男、2割女 男 X:0.2, Y:0.4 女 X:0.4, Y:0.2 0.2 vs. 0.25 10000人 女 model.1 <- my.make.model(1,1) model.2 <- my.make.model(1,2) model.3 <- my.make.model(1,0) n.pt <- 1000 X <- matrix(sample(1:2,n.pt,replace=TRUE,prob=c(0.8,0.2)),ncol=1) Prob.Vec <- list() tmp.Pr <- (X-1)*0.2+0.2 Prob.Vec[[1]] <- cbind(tmp.Pr,1-tmp.Pr) tmp.Pr <- (X-1)*(-0.2)+0.4 Prob.Vec[[2]] <-cbind(tmp.Pr,1-tmp.Pr) better.1 <- which(Prob.Vec[[1]][,1]>Prob.Vec[[2]][,1]) better.2 <- which(Prob.Vec[[1]][,1]<Prob.Vec[[2]][,1]) n.iter <- 10 out.mat.1 <- matrix(0,n.iter,length(better.1)) out.mat.2 <- matrix(0,n.iter,length(better.2)) for(i in 1:n.iter){ model.2.out <- my.simulate.model.2(model.2,X,Prob.Vec,n.iter=n.pt) better.1.selection <- (-1)*(model.2.out$selection[better.1]-2) better.2.selection <- model.2.out$selection[better.2]-1 out.mat.1[i,] <- cumsum(better.1.selection)/(1:length(better.1)) out.mat.2[i,] <- cumsum(better.2.selection)/(1:length(better.2)) } matplot(t(out.mat.1),type="l") matplot(t(out.mat.2),type="l") 男 800 200
1→1000人 8割男、2割女 男 X:0.2, Y:0.2 女 X:0.4, Y:0.2 女 n.pt <- 1000 X <- matrix(sample(1:2,n.pt,replace=TRUE,prob=c(0.8,0.2)),ncol=1) Prob.Vec <- list() tmp.Pr <- (X-1)*0.2+0.2 Prob.Vec[[1]] <- cbind(tmp.Pr,1-tmp.Pr) tmp.Pr <- (X-1)*0+0.2 Prob.Vec[[2]] <-cbind(tmp.Pr,1-tmp.Pr) better.1 <- which(X==1) better.2 <- which(X==2) n.iter <- 10 out.mat.1 <- matrix(0,n.iter,length(better.1)) out.mat.2 <- matrix(0,n.iter,length(better.2)) for(i in 1:n.iter){ model.2.out <- my.simulate.model.2(model.2,X,Prob.Vec,n.iter=n.pt) better.1.selection <- (-1)*(model.2.out$selection[better.1]-2) better.2.selection <- model.2.out$selection[better.2]-1 out.mat.1[i,] <- cumsum(better.1.selection)/(1:length(better.1)) out.mat.2[i,] <- cumsum(better.2.selection)/(1:length(better.2)) } matplot(t(out.mat.1),type="l") matplot(t(out.mat.2),type="l") 男 800 200
悪くない…（？） それ以外の選択戦略との比較が必要 ・・・割愛
やれやれ できた！
○率の高いYばかりが選ばれるようになった … and they lived happily ever after めでたし、めでたし Y × Y ○ ほんとに よかったのぉ、 おばあさん X × よかったですね、 おじいさん X ○
『いいことを教えてやろう』
『右の道を選んだ者、７名あり。４名は幸福に、３名は不幸になった』 幸福になった者の体重は67,53,86,71kg、不幸になった者の体重は48,52,51kgであった 『左の道を選んだ者、３名あり。２名は幸福に、１名は不幸になった』 幸福になった者の体重は41,53,49kg、不幸になった者の体重は88,68,64kgであった
やってくる人の「体重の分布」
Xが得 X,Yの成功率 Yが得
帰結ごとにカーネル推定 説明変数(体重)における「みなし観測度数」を推定 「みなし観測度数」に基づく「みなしベータ分布」 「みなしベータ分布」をX,Y道間で比較 「めずらしい体重」の人の場合「みなし観測度数」は小さくなり、「ありふれた体重」の人の場合は「みなし度数」が大きくなる
Xが良いはずの人 Ｙが良いはずの人
全10人 Xが良いはずの人 Ｙが良いはずの人
全50人 Xが良いはずの人 Ｙが良いはずの人
全125人 Xが良いはずの人 Ｙが良いはずの人
全250人 Xが良いはずの人 Ｙが良いはずの人
全500人 Xが良いはずの人 Ｙが良いはずの人
全1000人 Xが良いはずの人 Ｙが良いはずの人
Xが良いはずの人 Ｙが良いはずの人
悪くない…（？）
0.5 1/(1+r) 事前確率 事後確率 0.5 r/(1+r) 今回は「帰無仮説」は入れていない… 検討中
『いいことを教えてやろう』
『いいことを教えてやろう』 もう、いい！
『いいことを教えてやろう』 もう、いい！ というわけにもいかず
帰無・対立の事前・事後確率 男女 + その他複数の名義尺度 複数の量的変数 名義尺度と量的変数の組合せ
いくつかのこと 量的変数・多次元 多名義尺度における「仮説数」のハンドリング 帰結のカーネル分布推定が効かなくなる k-NN (k-nearest Neighbors)で代用できる？？ 多名義尺度における「仮説数」のハンドリング 2^k : k=10くらいまでは力技でも？？