ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

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

2016-08-25 ベルマン方程式

[][][][]ベルマン方程式

  • 調べものをしましたが、日本語サイトの中でこちらが一番、わかりやすく参考にさせていただきました
  • ベルマン方程式というものがある。Wikipediaではベルマン方程式(リチャード・E・ベルマン)、もしくはBellman_equation
  • 選択肢があるときに、どれを選ぶとよいかを迷いながら選ぶような状況があるとする
  • どれを選ぶかの基準は、「選ぶことで得られる利益の総和」を大きくしたい、という基準
  • この基準に照らして、もし最適な選択肢を取ることができるのであれば、その利益の総和はどうなるだろうか、というようなことを表した方程式
  • 選択肢を選んだときに、どのくらいの利益が生じるかは確率的に決まる(ある時は利益が出るし、あるときは出ない、とか)ので、常に、どの選択肢を選んだものか…と悩ましい状態に置かれるが、それに関する式だということ
  • この論文では、複数の選択肢があって、選ぶたびに利益は0か1かになるような単純な選択課題であるが、そこに出てくる式(3.1)を理解したい
  • 式は以下のような形をしている
    • V(x_{1,t},x_{2,t},...,x_{K,t})=
    •   max_{k} ¥{(profit_{k,t} +
    • d(q_{k,t} ¥times V(x_{1,t},...,x_{k,t}+(1,0),...) + (1-q_{k,t}) ¥times V(x_{i,t},...,x_{k,t}+(0,1),...)¥}
  • 各選択肢について、0回以上の試行結果の成否情報がある。それがx_{i,t}。成否の回数からなる長さ2の整数ベクトル。ただし、iは第i番目の選択肢。tは、これまでの全部でt試行がなされていることを意味する(複数の選択肢があるので、それぞれが実行された回数の総和になる)。
  • 成否情報があったときに、最適選択をするとすると、その先で得られる利益の総和(の期待値)がV(x_{1,t},...)。全試行数がtのときのそれ
  • このような状態(x_{i,t})にあって、最善選択肢が第k番だとなっている、とする。
  • そうすると、全試行数がt+1に増えた段階では、第k番選択肢を選ぶことで得られる利益の期待値proit_{k,t}が利益となることはわかる。ただし、この追加利益は未知であるから、tの時点でのその推定値を使うことになる。二値成否の場合には、ベータ分布から知られる期待値を使うのがよいのでprofit_{k,t} = ¥frac{s_{k,t}+1}{s_{k,t}+f_{k,t}+2}、というように成否回数s,fと下駄に当たる、1、2を使った期待値を使っているのが、リンク先の論文
  • この第k番選択肢が選ばれることがわかっているものとして、t+1以降にどんな利益が続くのかは不明だが、それをVで表してしまっている。ただし、このt+1以降で考慮するべきなのは、tにて第k番選択肢を選んで、成功する確率qと失敗する確率(1-q)とがあって、それぞれの場合について、起きる確率の重みを考慮して加える必要がある
  • また、このqについては、ベータ分布にて成功確率を推定するとすれば、q= ¥frac{s_{k,t}+1}{s_{k,t}+f_{k,t}+2}
  • さて、こんな風なVに関する確率コミの漸化式ができたので、これを用いて、推定することを考える
  • 全体で最適選択をしたら、どんな総利益が出るかを一発で出すことはできないけれど、漸化式をちまちま解いていくと最終的には答えがでるね、このように全体の答えを解くのに、ちまちまと部分部分で考えるのが動的計画法(こちら)
  • ちなみにdは、繰り返して実行するにしても、そのうち止めるという現実的な条件を入れるためのたで、継続する確率をdとしている。無限に続けるならd=1
  • こちらにベルヌーイ事象のGittins indexのGlazebrook Weberアルゴリズム(ベタなアルゴリズム)のMatlabコードがある。それのベタアルゴリズムの方のR化したもの(バグが潜んでいる可能性あり)。一応、答え合わせは合っている感じ…
  • ぐるぐるとループが3重に回っているので遅いです
  • 0-1の値についてグリッド探索すると得られることが知られているらしい。ペイパーのAppendixの書いてあるのだが…、以下のソースがそれを実装しているのかどうか、ちょっと心もとないが、結果があっているのでよしとする
  • 成功率は不明であるが、それがわかりさえすれば、繰り返して実施したときの総報酬期待値はその成功率に関して、1/(1-d)倍したものになる
  • なので、その成功率をグリッドで定め多ものに基づく総報酬期待値について、小さい方から、試してみる。で、試した値が、漸化式の値に近いことが解れば、その総報酬期待値を「当該観測成否回数」における総報酬期待値とする
  • 漸化式の値に近いことは、「成功率グリッド値を小さい方から試して行って、初めて、漸化式の値を超えたあたり」という判定方法とする
  • そのようにしたときの成功率グリッド値(グリッドの幅を考慮して少し加減した値)をGittins indexとして登録する
my.gittins <- function(d,Mh,grid){
	G <- matrix(0,Mh-1,Mh-1)
	Gittins <- matrix(0,Mh-1,Mh-1)
	
	for(s in 1:(Mh-1)){
		G[s,Mh-s] <- s/Mh
	}
	ps <- seq(from = grid/2, to =1, by=grid)
	for(ii in 1:length(ps)){
		p <- ps[ii]
		# 1/(1-d) = 1 + d + d^2 + .....
		# p/(1-d)はdの確率で継続するときの…、全報酬の期待値
		known <- p/(1-d)
		for(i in (Mh-1):2){
			for(s in 1:(i-1)){
				experimental <- s/i * (1+d*G[s+1,i-s]) + (i-s)/i * d * G[s,i-s+1]
				if(Gittins[s,i-s] ==0 & (known > experimental)){
					Gittins[s,i-s] <- p - grid/2
				}
				G[s,i-s] <- max(known,experimental)
			}
		}
	}
	

	
	
	return(list(G=G,Gitting=Gittins))
}

Git09 <- my.gittins(0.9,200,0.001)
Git09[[2]][1:8,1:8]
  • 全試行回数を有限に設定し、そこから遡る形で計算
  • 結果
> Git09[[2]][1:8,1:8]
      [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
[1,] 0.703 0.500 0.380 0.302 0.249 0.210 0.181 0.159
[2,] 0.800 0.635 0.516 0.434 0.372 0.324 0.287 0.257
[3,] 0.845 0.707 0.601 0.518 0.456 0.406 0.365 0.331
[4,] 0.872 0.754 0.658 0.581 0.518 0.468 0.426 0.390
[5,] 0.891 0.787 0.700 0.628 0.568 0.517 0.475 0.439
[6,] 0.904 0.811 0.732 0.664 0.607 0.558 0.516 0.479
[7,] 0.914 0.831 0.757 0.694 0.640 0.592 0.551 0.514
[8,] 0.922 0.846 0.778 0.719 0.667 0.621 0.581 0.545
  • Gittins indexとして論文に掲載されているもの
g1 <- c(0.7029,0.8001,0.8453,0.8723,0.8905,0.9039,0.9141,0.9221)
g2 <- c(0.5001,0.6346,0.7072,0.7539,0.7869,0.8115,0.8307,0.8461)
g3 <- c(0.3796,0.5163,0.6010,0.6579,0.6996,0.7318,0.7573,0.7782)
g4 <- c(0.3021,0.4342,0.5184,0.5809,0.6276,0.6642,0.6940,0.7187)
g5 <- c(0.2488,0.3702,0.4561,0.5179,0.5676,0.6071,0.6395,0.6666)
g6 <- c(0.2103,0.3245,0.4058,0.4677,0.5168,0.5581,0.5923,0.6212)
g7 <- c(0.1815,0.2871,0.3647,0.4257,0.4748,0.5156,0.5510,0.5811)
g8 <- c(0.1591,0.2569,0.3308,0.3900,0.4387,0.4795,0.5144,0.5454)

G.mat <- rbind(g1,g2,g3,g4,g5,g6,g7,g8)
> t(G.mat)
         g1     g2     g3     g4     g5     g6     g7     g8
[1,] 0.7029 0.5001 0.3796 0.3021 0.2488 0.2103 0.1815 0.1591
[2,] 0.8001 0.6346 0.5163 0.4342 0.3702 0.3245 0.2871 0.2569
[3,] 0.8453 0.7072 0.6010 0.5184 0.4561 0.4058 0.3647 0.3308
[4,] 0.8723 0.7539 0.6579 0.5809 0.5179 0.4677 0.4257 0.3900
[5,] 0.8905 0.7869 0.6996 0.6276 0.5676 0.5168 0.4748 0.4387
[6,] 0.9039 0.8115 0.7318 0.6642 0.6071 0.5581 0.5156 0.4795
[7,] 0.9141 0.8307 0.7573 0.6940 0.6395 0.5923 0.5510 0.5144
[8,] 0.9221 0.8461 0.7782 0.7187 0.6666 0.6212 0.5811 0.5454

2016-08-18 ExAC

[][][][]ExAC

  • Analysis of protein-coding genetic variation in 60,706 humans
  • 91000エクソーム→フィルタリング・QC→60706人
  • ヨーロッパ人を大多数に、東アジア、南アジア、アフリカ、Latinoを主要サンプルとする
  • 10195872(1千万)候補バリアントを見つけ、確度の高い7404909(七百万)を高品質バリアントとしてレポート(うち、317381はins/dels)
  • 平均密度1バリアント/8塩基対(エクソーム標的)
  • マイナーアレル頻度:1%未満が99%超、シングルトンが54%
  • 機能性との関係(アミノ酸置換なし vs. ありでの違い。ins/delsにおけるフレームシフトなし vs. フレームシフトあり)
  • バリアントの生じやすさとの関係:CpG(C->T置換)は多い
  • Trialleleic (7.9 %) はポアソンモデルからそれほど遠くない
  • Infinite site仮説は本当ではなくて、同じ場所に同じタイプの変異が人類史上、複数回起きていることの確認も取れた
  • 近接バリアントによる機能レスキューも検出
  • X染色体遺伝子配列の保存性の高さ(常染色体遺伝子に比べ)は追認
  • 機能喪失型バリアントの入り方により機能パスウェイの『必須度』を測ると、「必須度が高い」のは、spliceosome, ribosome,proteasome)、反対に「必須度が低い」のはolfactory receptors(追認)
  • この機能上必須度の高い遺伝子が病気の原因バリアントを持つわけではないことも確認(ある程度の多様性を持ちうるのは必須度が『そこそこ』な遺伝子…)
  • 既報のメンデル型遺伝病原因バリアントは、フェノタイプフリーな個人にも多数認められており、このフェノタイプフリーが浸透率の低さのせいなのか、病原性のなさのせいなのかなどにつき、今後、原因バリアントの整理が必要
  • Protein-truncating variants(PTV)が7404909バリアントのうち179774(十八万)個。平均すると一人当たり、PTVをホモで35遺伝子に、ヘテロで85遺伝子に持つことになる

2016-08-09 ぱらぱらめくる『逆問題の数学』

[][][][][]ぱらぱらめくる『逆問題の数学』

逆問題の数学

逆問題の数学

  • 本の前提
    • リッジ回帰は、y=X¥betaのようなデータについて¥betaの値が説明変数にばらけるような割り付けをするタイプの回帰
    • そのときに||y-X¥beta||2 + ¥alpha||¥beta||_2の最小化をする
  • チホノフ正則化、Tikhonov正則化は、||.||のノルムの定義をもう少し一般化したもの(のようだ)
    • 一般化してしまうと、どんな内積・ノルムを仮定するかが、俄然自由になってしまって、広い「汎関数空間」の探索作業になってしまうTikhonov一派は、「経験的に?」どういうアルファを設定するのがよいか、とかそういうことについての提唱をしている(らしい)
  • まえがきを読むとだいたい、この本はわかる
    • 逆問題と言うのがある。データから、世界がどうなっているのかを理解しようとする、そんな問題設定
    • 非線形だったりして、やっかいだし、ill-posed problemと呼ばれるような、解が不安定で解きにくい(解けない場合)もある
  • その逆問題だが、実用的な課題でもあるため、工学諸分野でそれぞれ問題設定などがなされている(生物分野もそうだろう)
  • それとは別に、数学としての逆問題というものある
  • 数学から言えば、逆問題は、「関数解析」として記述されるよ、となる
  • チホノフ正則化、と言うのが、ある一定の役割を持つ
  • チホノフ正則化を含めた正則化は、「正則化してその前提で解を得よう」という立場だから、「解がないことは頭にない」
  • 逆問題の解の存在を示すことはそもそも難しく、『解がある』前提で解くので、結果を信じるかどうかは、別途考えるべきこと
  • 一般的な逆問題の他に、時系列データで差分作用素を考えるとき、ヤコビアンを使うが、ここでの逆問題では、解の存在と一意性があるという。ただし、非線形なので、解が安定かどうかはわからない。関数解析ではあるけれど、線形代数で押し通せる(ようだ)
  • そして、いくつかの物理・工学領域での逆問題を数学的スタンスで書いている章が続く

2016-08-08 ぱらぱらめくる『はじめての統計データ分析』5.6.

[][]5.6 複数要因・組み合わせ、独立だったり順序があったり ぱらぱらめくる『はじめての統計データ分析』

  • stanファイル似て指定する
  • モデルを書き換えて、知りたいことを生成量にすれば、なんでもできる
  • 要因が増えるなら、モデル式の項数を増やす
  • 複数カテゴリに順序を入れてそれらが守られているかどうかは、事後分布から発生した標本のうちそれが満足されるものを評価すれは良いので、generated quantities 生成量 を指定する

2016-08-07 ぱらぱらめくる『はじめての統計データ分析』3.4.

[][]3.4. 2群のデータ ぱらぱらめくる『はじめての統計データ分析』

  • t検定が、等分散性の仮定の有無を区別するように、また、対応のありなしを区別するように、事後分布推定もそれぞれの仮定の下で行う
  • 仮定はmodelとして指定する
    • 等分散性の有無は、2群の分散を1つのパラメタで表すか否かで書き分ける
    • 対応があることは、分散共分散行列の係数推定を入れることで表す
data {
  int<lower=0> n1;     int<lower=0> n2;       #データ数
  real x1[n1];         real x2[n2];           #データ
  real EQU;                                   #SD、1 - 共通、0 - 異なる
  real mL;real mH;real sL;real sH;            #事前分布
}
parameters {
  vector<lower=mL,upper=mH> [2] mu;          #平均(範囲指定)
  real<lower=sL,upper=sH> sigma1;            #標準偏差(範囲指定)
  real<lower=sL,upper=sH> dummy;             #ダミーのSD(範囲指定)
}
transformed parameters {
  real<lower=0>           sigma2;             #標準偏差2
  sigma2<-if_else(EQU>0.5,sigma1,dummy);
}
model {
  x1 ~ normal(mu[1],sigma1);                  #正規分布1
  x2 ~ normal(mu[2],sigma2);                  #正規分布2
}
generated quantities{
  real xaste[2];
  real log_lik;
  xaste[1]<- normal_rng(mu[1],sigma1);        #予測分布1
  xaste[2]<- normal_rng(mu[2],sigma2);        #予測分布2
  log_lik<-normal_log(x1,mu[1],sigma1)+
           normal_log(x2,mu[2],sigma2);       #対数尤度
}
data { 
  int<lower=0> n;                             #データ数   
  vector[2] x[n];                             #データ     
  real EQU;                                   #SD、1 - 共通、0 - 異なる
  real mL;real mH;real sL;real sH;            #事前分布
}
parameters {
  vector<lower=mL,upper=mH> [2] mu;          #平均(範囲指定)
  real<lower=sL,upper=sH> sigma1;            #標準偏差(範囲指定)
  real<lower=sL,upper=sH> dummy;             #ダミーのSD(範囲指定)
  real<lower=-1,upper=1>  rho;                #相関
}
transformed parameters {
  real<lower=0>           sigma2;             #標準偏差2
  cov_matrix[2]           Sigma;
  sigma2<-if_else(EQU>0.5,sigma1,dummy);
  Sigma[1,1] <- pow(sigma1,2);                #共分散行列
  Sigma[2,2] <- pow(sigma2,2);
  Sigma[1,2] <- sigma1*sigma2*rho;
  Sigma[2,1] <- Sigma[1,2];
}
model {
  for(i in 1:n){x[i]~multi_normal(mu,Sigma);} #2変量正規分布
}
generated quantities{
  vector[2] xaste;
  real log_lik;
  xaste <-  multi_normal_rng(mu,Sigma);       #予測分布
  log_lik <- multi_normal_log(x, mu, Sigma);  #対数尤度
}