決断のための分布合算 京大(医)統計遺伝学分野 山田　亮.

Presentation on theme: "決断のための分布合算 京大(医)統計遺伝学分野 山田　亮."— Presentation transcript:

1 決断のための分布合算 京大(医)統計遺伝学分野 山田　亮

2

3 小児～若年発症 ループス 血栓リスクが高そうに思うが… 大急ぎで電子カルテDBを検索した：かき集めても血栓例は10例
N Engl J Med 2011; 365: November 10, 2011DOI: /NEJMp

4 わからなくても決断する あなたは冒険旅行中 分かれ道があって、電光掲示板がある 11例目のあなたは、どちらの道を選ぶか
『右の道を選んだ者、７名あり。４名は幸福に、３名は不幸になった』 『左の道を選んだ者、３名あり。２名は幸福に、１名は不幸になった』 11例目のあなたは、どちらの道を選ぶか 我を過ぐれば憂ひの都あり、 我を過ぐれば永遠の苦患あり、 我を過ぐれば滅亡の民あり

5 ○ × 和 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)

6 期待値で選択することは「悪くない」 方針 「期待値」が大きい方を選ぶ 「期待値」が同じなら、どちらかを選ぶ

7 ゼロからのスタート 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")

8 X Y X期待値 Ｙ期待値 ○ × 0.5 1 2 0.4 3 4 0.25 5 0.2 6 7 8 9 0.375 10 11

9 Y × Y ○ X × X ○

10 Y × Y ○ X × X ○

11 ○率の高いYばかりが選ばれるようになった
… and they lived happily ever after  めでたし、めでたし Y × Y ○ X × X ○

12 ○率の高いYばかりが選ばれるようになった
… and they lived happily ever after  めでたし、めでたし Y × Y ○ X × X ○

13 本当にYばかりが 選ばれるようになるか

14 Ｓｅｌｅ

15

16 Yに落ち着くか Xに落ち着くか Ｙ Ｘ

17 確率的な決断 Multi-armed bandit 問題
複数のスロットマシンがあって、それぞれのマシンには「当たり」の確率が決まっているが、その確率が不明であるという マシンを１つずつ選んでは、勝負をして、各マシンの当否結果を記録しながら、勝負を繰り返すことにする どんなルールで選ぶと、儲けが最大になりやすいか、という問題

18 確率的な決断 Multi-armed bandit 問題
その状況でのThomson samplingとかの方が良い結果が得られることが知られている。 ごく大雑把に言うと、 データを見ても、「100％、どのアームがよいとは言い切れない」から、データから見て、「得策らしくないアームも、ある程度(確率的に)は選ぼう」 いったん、悪い方を選び勝ちになっても、判断を修正するポテンシャルが「確率的な決断」によってもたらされる

19 『確率的な決断にする！』

20 わからなくても決断する 分かれ道があって、電光掲示板がある 11例目のあなたは、どうするか 『確率的な決断が大事である』
『従って、この分かれ道に奇数回目に来た者には、電光掲示板は点灯せず、偶数回目に来た者には、点灯することとする』 11例目のあなたは、どうするか １．適当に選ぶ ２．出直す

21 わからなくても決断する 分かれ道があって、電光掲示板がある 11例目のあなたは、どうするか 『確率的な決断が大事である』
『従って、この分かれ道に奇数回目に来た者には、電光掲示板は点灯せず、偶数回目に来た者には、点灯することとする』 11例目のあなたは、どうするか (１) 適当に選ぶ→こうすればある程度の「確率的決断」が起きる ２．出直す

22 わからなくても決断する 分かれ道があって、電光掲示板がある 11例目のあなたは、どうするか 『確率的な決断が大事である』
『従って、この分かれ道に奇数回目に来た者には、電光掲示板は点灯せず、偶数回目に来た者には、点灯することとする』 11例目のあなたは、どうするか (１) 適当に選ぶ→こうすればある程度の「確率的決断」が起きる (２) 出直す

23 何が問題か 偶数・奇数で不平等

24 何が問題か 偶数・奇数で不平等 平等な何かを作ってみる

25 何を比較する？ 「どちらの道を選ぶと○になる確率が高いのか」 「どちらの道が『○率が高い』のか」 これは○の期待値
「どちらの道の『○の期待値』が高いのか」 「どちらの道が『○率が高い』のか」

26 何を比較する？ 「どちらの道を選ぶと○になる確率が高いのか」 「どちらの道が『○率が高い』のか」 これは○の期待値
「どちらの道の『○の期待値』が高いのか」 「どちらの道が『○率が高い』のか」 違う設問だから答えも違う

27 何を比較する？ 「どちらの道を選ぶと○になる確率が高いのか」 「どちらの道が『○率が高い』のか」 これは○の期待値
「どちらの道の『○の期待値』が高いのか」 「どちらの道が『○率が高い』のか」 違う設問だから答えも違う Xの方が高いかもしれないし、Ｙの方が高いかもしれない

28 ○ × 和 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)

29 ○ × 和 X 4 3 7 Y 2 1 期待値　(4+1)/(7+2) 期待値(2+1)/(3+2)

30 「どちらの道の『○の期待値』が高いのか」
× 和 X 4 3 7 Y 2 1 期待値　(4+1)/(7+2) 期待値(2+1)/(3+2) 「どちらの道の『○の期待値』が高いのか」

31 ○ × 和 X 4 3 7 Y 2 1 (4+1)/(7+2) 「どちらの道が『○率が高い』のか」

32 ○ × 和 X 4 3 7 Y 2 1 (4+1)/(7+2) (2+1)/(3+2) 「どちらの道が『○率が高い』のか」

33 計算できる(式の導出は省略)

34 ○率が高い確率応じて X,Yを「確率的に」選択してみよう

35

36 ○率の高いYばかりが選ばれるようになった
… and they lived happily ever after  めでたし、めでたし Y × Y ○ X × X ○

37 『いいことを教えてやろう』

38 『右の道を選んだ者、７名あり。４名は幸福に、３名は不幸になった』
その内訳は 男　５名。２名は幸福に、３名は不幸に 女　２名。２名とも幸福に 『左の道を選んだ者、３名あり。２名は幸福に、１名は不幸になった』 男　３名。２名は幸福に、１名は不幸に 女　は左の道を選んでおらぬ

39

40 男女合算 と 女のみ 女のみ 男女合算 女のみ 男女合算 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)) 男女合算

41 「男女に違いなし」なら 「男女に違いあり」なら 男女合算の情報を使った方が正確 男女合算の情報に基づいて集計した方が、早く、収束する
男女別々の情報を使った方が正確

42 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] > Decision_beta.2(c(2,0,0,0)+1) [1] 0.75 0.42 0.75

43 男女合算 女のみ 道の選択確率が異なる 『道 Ｘ、道 Ｙ、どちらにしよう？』 『男女合算、女のみ、どちらにしよう？』 確率的に選んだ
男女合算　女のみ 道の選択確率が異なる 『道　Ｘ、道　Ｙ、どちらにしよう？』 確率的に選んだ 『男女合算、女のみ、どちらにしよう？』 確率的に選んでみる

44 男女合算 女のみ 道の選択確率が異なる 『道 Ｘ vs. 道 Ｙ、どちらにしよう？』 『男女合算 vs. 女のみ、どちらにしよう？』
男女合算　女のみ 道の選択確率が異なる 『道　Ｘ vs. 道　Ｙ、どちらにしよう？』 確率的に選んだ 『男女合算 vs. 女のみ、どちらにしよう？』 確率的に選んでみる

45 男女合算 道X,道Yのベータ分布 Y X X

46 男女合算→道X 道X→男,道Y→女のベータ分布

47 X ○ × 和 男 2 3 5 女 道X 男,女のベータ分布 女 X 男

48 X ○ × 和 男 2 3 5 女 道X 男,女のベータ分布 (p_m=p, p_f=p) 女 X 男

49 X ○ × 和 男 2 3 5 女 道X 男, 女のベータ分布 (p_m,p_f) (p_m=p, p_f=p) 女 X 男

50 (p_m,p_f) 『男女に違いがあってもよい』という立場 (p_m=p, p_f=p) 『男女に違いがない』という場合 女 X 男

51 X ○ × 和 男 2 3 5 女 女 X 男

52 X ○ × 和 男 2 3 5 女 女 X 男

53 X ○ × 和 男 2 3 5 女 女 この積分値を活用する X 男

54 『男女差がある場合』 サンプル数が０：積分値１ サンプル数が増える：小さくなる 『男女差がない場合』より 小さくなる程度は急速

55 小さくはなるが、対立仮説が作る同時分布のピークが対角線のごく近傍に来るので、小さくなっても困らない
『男女差がない場合』 サンプル数が０：積分値１ サンプル数が増える：小さくなる 小さくはなるが、対立仮説が作る同時分布のピークが対角線のごく近傍に来るので、小さくなっても困らない

56 X ○ × 和 男 2 3 5 女 どう活用する？ X

57 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女

58 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女

59 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女

60 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差あり』 の同時分布 正方形部分の積分は１

61 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差あり』 の同時分布 正方形部分の積分は１

62 『男女に差あり』 の同時分布 正方形部分の積分は こちらも１ 『男女に差あり』 の同時分布 正方形部分の積分は１ X ○ × 和 男 女 X
女 X ○ × 和 男 2 3 5 女 『男女に差あり』 の同時分布 正方形部分の積分は こちらも１ 『男女に差あり』 の同時分布 正方形部分の積分は１

63 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女

64 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差なし』 部分

65 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差なし』 部分

66 X ○ × 和 男 女 X ○ × 和 男 2 3 5 女 『男女に差なし』 部分

67 １ １ １

68 １ １ 積分 １

69 ？ ？ 事前確率 事後確率 ？ ？

70 0.5 ？ 事前確率 事後確率 0.5 ？

71 １ １ 積分 １

72 0.5 1/(1+r) 事前確率 事後確率 0.5 r/(1+r) =r

73 仮説の比率が決まれば ２つのベータ分布の 重みづけ混合分布 1/(1+r) Beta(a+c+1,b+d+1) + r/(1+r) Beta(c+1,d+1)

74 道X、道Y、 それぞれの混合ベータ分布 に基づき 『〇率が高くなる確率』 による確率的選択 モンテカルロに 『○率が高くなる確率』を算出

75 やってみよう

76 女 男 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)* 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

77 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)* 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

78 女 男 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)* 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

79 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)* 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

80 悪くない…（？） それ以外の選択戦略との比較が必要 ・・・割愛

81 やれやれ できた！

82 ○率の高いYばかりが選ばれるようになった
… and they lived happily ever after  めでたし、めでたし Y × Y ○ ほんとに よかったのぉ、 おばあさん X × よかったですね、 おじいさん X ○

83 『いいことを教えてやろう』

84 『右の道を選んだ者、７名あり。４名は幸福に、３名は不幸になった』
幸福になった者の体重は67,53,86,71kg、不幸になった者の体重は48,52,51kgであった 『左の道を選んだ者、３名あり。２名は幸福に、１名は不幸になった』 幸福になった者の体重は41,53,49kg、不幸になった者の体重は88,68,64kgであった

85 やってくる人の「体重の分布」

86 Xが得 X,Yの成功率 Yが得

87 帰結ごとにカーネル推定 説明変数(体重)における「みなし観測度数」を推定 「みなし観測度数」に基づく「みなしベータ分布」
「みなしベータ分布」をX,Y道間で比較 「めずらしい体重」の人の場合「みなし観測度数」は小さくなり、「ありふれた体重」の人の場合は「みなし度数」が大きくなる

88 Xが良いはずの人 Ｙが良いはずの人

89 全10人 Xが良いはずの人 Ｙが良いはずの人

90 全50人 Xが良いはずの人 Ｙが良いはずの人

91 全125人 Xが良いはずの人 Ｙが良いはずの人

92 全250人 Xが良いはずの人 Ｙが良いはずの人

93 全500人 Xが良いはずの人 Ｙが良いはずの人

94 全1000人 Xが良いはずの人 Ｙが良いはずの人

95 Xが良いはずの人 Ｙが良いはずの人

96 悪くない…（？）

97 0.5 1/(1+r) 事前確率 事後確率 0.5 r/(1+r) 今回は「帰無仮説」は入れていない… 検討中

98 『いいことを教えてやろう』

99 『いいことを教えてやろう』 もう、いい！

100 『いいことを教えてやろう』 もう、いい！ というわけにもいかず

101 帰無・対立の事前・事後確率 男女 + その他複数の名義尺度 複数の量的変数 名義尺度と量的変数の組合せ

102 いくつかのこと 量的変数・多次元 多名義尺度における「仮説数」のハンドリング 帰結のカーネル分布推定が効かなくなる
k-NN (k-nearest Neighbors)で代用できる？？ 多名義尺度における「仮説数」のハンドリング 2^k : k=10くらいまでは力技でも？？