お試し並列処理〜細工〜擬似乱数列の管理

パーミューテーションテストでモンテカルロを使うとすると、並列プロセスで用いる擬似乱数列も相互にランダムであることが必要である。シードから開始する擬似乱数列を用いる場合には、採取する擬似乱数を、プロセス数おきに採取させることによって、全プロセスの擬似乱数列が、単一の擬似乱数列からの部分列とすることができる。このようにすることで、擬似乱数列のランダム性を確保するだけで、全プロセスの擬似乱数列のランダム性の精度管理を行うことができる。
擬似乱数列ジェネレータmzを作成し、x=ak+(b-1), k=0,1,2,...1番目の乱数を取る

int seed=xxxxxx;
Utils.MersenneTwisterFast mz=new MersenneTwisterFast(seed);
public static void shuffleDoubleMZ(double[] arr,MersenneTwisterFast mt,int a,int b){

for(int i=0;i<b;i++){
	mt.nextDouble();//空回しb回(剰余分)
}
for(int i=0;i<loop;i++){
	double r=mt.nextDouble();
	for(int j=0;j<a;j++){
		mt.nextDouble();//空回しa回周期
	}
}

お試し並列処理〜小細工〜システム一時ファイルでプロセスの進行制御

  • 管理ファイルの存在・非存在ループ
    • 第1プロセスは管理ファイルを作成する。このファイルは、オリジナル検定セットが終了してから作成する。このファイルが非存在の間は、空ループを回し、k-1プロセスは先へ進まない
String dir="hogedir";
String regulationfile=regulator;
File regulator = new File(dir,p.regulationfile);
 while(!regulator.exists()){
}
    • 第1プロセスは管理ファイルを作成し、パーミュテーション処理の参照データを管理ファイルに書き込む。この情報を読み込んで、適切な参照データが取り込まれるまで、この読み込み作業を繰り返す
 boolean minPOK=true;
 while(minPOK){
  InputStream is = new FileInputStream(regulator.getPath());
  BufferedReader br=new BufferedReader(new InputStreamReader(is);
/*
読み込み
*/
  br.close();
  if(読み込み情報が適切?){
   minPOK=false;
  }
 }
    • パーミュテーション検定セットが終わったら、管理ファイルの読み取り・更新のための待ち行列に並ぶ。
File tmpregulator = new File(p.outdir, tmpfile);
while(tmpregulator.exists()){
/*
システム一時ファイルが存在したら、待つ
*/
}
if(!tmpregulator.exists()){
/*
システム一時ファイルが存在していなければ、ループを抜けて、この分岐に到達する。
??順番待ちがごった返すと、ループを抜けて、到達した後にも衝突が起きてトラブルになることがあるようなので、再度、システム一時ファイルが存在していないことをここで確認している(もしかしたら不要)
*/
boolean success=false;
/*
管理ファイルの専有許可が下りたので、システム一時管理ファイルを作成して、ほかのプロセスの進行をロックする。
ただし、ここでも競合することがあるので、その成否はtry{ }catch(IOException e){}
でくくり、失敗してもよい
失敗の場合は、プロセスを終了させることなく、読み取り・書き込みをしないで、再度、順番待ち行列に戻ることとする(パーミュテーションループの回数を1回分元に戻し、パーミュテーション検定セットを実行しないで、戻ってくる
*/
try{
 success = tmpregulator.createNewFile();
}catch (IOException e){
}
if(success){
/*
管理ファイルの読み取り、更新
*/

/*
システム一時ファイルの確実な消去
このファイルは作ったら、必ず消さないと、すべてのプロセスが止まるので
確実に消す必要がある
失敗したまま先に進むわけにいかない
ここで、単に、システム一時ファイルを消す、とすると
失敗することもあるようなので、以下のような「存在する限り」消す
というループとした
また、これでも、ループの中がファイルの消去コマンドのみだと失敗することが
あるようなので、意味のない処理をはさんだ
ここは、単に時間を稼ぐだけでは、失敗を防げないらしい
*/
 while(tmpregulator.exists()){
  String dummy="0";
  for(int d=0;d<10;d++){
   dummy+=dummy;
  }
  tmpregulator.delete();
 }
}

お試し並列処理〜システム一時ファイルでプロセスの進行制御

参考にしたのは、こちらの記事
javaの1プロセスの中に複数のスレッドを立てたりして制御するには、こんな処理ももちろんあるわけだけれども、今は、複数のノードでの制御なので、Javaの外に制御用存在を置く

  • 想定実行条件
    • (1)今、kxN回のパーミュテーション処理をkプロセスに振り分ける状況を想定する
    • (2)第1プロセスにオリジナル検定セットを実行させ、その終了後に引き続きN回のパーミュテーション検定セットの実行をさせ、残りのk-1プロセスにはN回のパーミュテーション検定セットの実行をさせる。
    • (3)kプロセスはパーミュテーション試行において、オリジナル検定セットの結果を用いるので、オリジナル実行を行わないプロセスは、オリジナル検定セットの結果が出るまで、処理の開始を保留する
    • (4)第1プロセスは、オリジナル検定セットの結果をファイルに書き出し、k-1プロセスはこのオリジナル検定セット結果ファイルを読み取ることで、パーミュテーション試行結果をオリジナルのそれと比較することを可能にする
    • (5)各ノードはパーミュテーション検定セットが終了するたびに、処理の全体の進行と集計を管理するファイル(管理ファイル)を書き換えることによって、全プロセスの処理を統合した結果がファイルにてアップデートされる
    • (6)また、各プロセスは、管理ファイルを読み取ることで、自プロセスの継続・終了の判断を下す
    • (7)管理ファイルは、複数のプロセスが読み書きすることから、この読み書きは常に1プロセスのみに許可されることとする
  • 実装上の要件
    • (1)複数のプロセスは独立したコマンドにて発行する。ただし、そのプロセスは、総プロセス数と、各プロセスのプロセスIDを引数としてとる(通常の並列処理と同様)
    • (2)複数のプロセスは、シェルスクリプトにより、実行ノードが明示的に指定される。また、個々のプロセスは、実行の継続・終了の判断を、全プロセスの結果の集計結果ファイルを参照することによって、自律的に行えるように、個々のプロセスコードを工夫する
    • (3)k-1プロセスは、第1プロセスがオリジナル検定セットの結果が出るまで、開始が保留される。オリジナル検定セットは全体プロセス管理ファイルに決められた書式で書き出される決まりとし、その管理ファイルを参照し、書式に合致するまで空ループを回す
    • (4)パーミュテーション処理は、個々のプロセスで最大回数を指定し、繰り返し実行する。個々のパーミュテーション処理が終わったら、管理ファイルを読み込み、当該処理の終了以前の全プロセス集計結果を読み取り、当該処理の結果を加え、その更新結果を管理ファイルに書き込む
    • (5)複数プロセスの個々のパーミュテーション結果後に行われる、管理ファイルの読み取りと更新は、競合することは、管理ファイルの更新に失敗するという問題と、更新内容に齟齬が生じるという2点から、1度に1プロセスにのみ許可される仕組みが必要である。この仕組みを、システムファイルの作成・消去によって実現する。システム一時ファイルが存在しているときは、あるプロセスが管理ファイルの読み取り・更新をしているときであることとし、このシステム一時ファイルが存在しているときには、処理の進行が保留される。また、システム一時ファイルが存在していないとき、各プロセスは管理ファイルの読み取り・更新を行ってよいことになるので、システム一時ファイルの非存在を確認したら、まず、システム一時ファイルを作成し、そののち、管理ファイルの読み取りと更新を行い、その終了後に、システム一時ファイルの消去を行う。これにて、各プロセスは当該パーミュテーション試行の実行と記録が終了したことになり、次の試行に移行する
    • (6)終了。管理ファイルを読み込み、終了条件に合致していたら、自プロセスを終了する。すべてのプロセスは、順次、終了条件の合致を知ることとなり、最終的に、全プロセスが終了する

お試し並列処理〜パーミュテーションテストで

  • パーミュテーションテスト
    • 処理構成
      • オリジナル検定セットの実行
      • パーミュテーション試行による検定セットの実行
      • パーミュテーション試行による検定セットの結果とオリジナル検定セットのそれとを比較して、集計
    • 並列化
      • 処理の単位は、検定セットの実行
      • 検定セットの実行はかなり大きい処理で、均質
      • 回数を稼ぎたいのも検定セットの実行
      • したがって、並列化は次のようになる
        • (1)オリジナル処理x1
        • (2)パーミュテーション試行による検定セット処理N回(N>1000)を利用可能な並列環境に分配する
        • (3)分配した並列環境での部分集計を取りまとめる
    • 並列化の実際
      • 並列プログラムにより、(2)の分配と(3)の集計とを行うことはもちろん可能である
      • しかしながら、異なる並列機器環境での移植の容易さを実現するために、並列仕様に極力依存しないことが望ましいとする
      • 並列機器間での通信が可能であり、ファイル共有が可能である、というのは、並列環境として、最低限の設定であるとみなせるので、それのみを用いることとする

ケースコントロール関連以外でパーミュテーションテスト

こちらのレビュー(Crit Care. 2004; 8(3): 196–199.)が簡潔にまとまっている。
Statistics review 10: Further nonparametric methods
ちなみにこのレビューシリーズは、全14回(継続中?)
このシリーズ1−14回へのリンクはこちら

  • データ
    • SNPのジェノタイプは3つ
    • それぞれのジェノタイプについて、フェノタイプが「あり か なし か」ではなくて、「量」。
    • 「量」というとき、いわゆるQuantitative traitのようにさまざまな値をとる場合と、順序量の場合とがある。また、さまざまな値をとる場合にも、その値の分布が正規分布やその他の、「きれいな」分布をとる場合と、まったくそうなっていない場合とがある。今、「あり か なし か」で扱えないデータをひとまず一括して扱ってしまいたい、というのが要請であるから、分布を仮定しないで、順序量にも使ってしまう方法を採用する
    • ちなみに、「正規分布を仮定」できるようなQuantitative traitについて、複数のカテゴリ(この場合はジェノタイプ)の間の違いを調べ、そのときに、ジェノタイプ3種類には序列を持ち込まないとしたら、それはANOVA。線形の関係を持ち込むなら、linear regression
    • SNPのジェノタイプに序列を持ち込まずに、各ジェノタイプに観測される、「きれいでない」量的データ(順序量を含む)に対して行う検定は、Kruskal-Wallis。
    • SNPのジェノタイプに序列を持ち込むと、Jonckheere-Terpstra。
  • ケース・コントロール分割表への(お試し)適用
    • 今、あるアレルに着目し、その本数によって、3ジェノタイプの量的情報を、「0,1,2」とする
    • 今、ある質的形質に着目し、そのありなし、を「0,1」とする。
    • 2つの見方ができる。
      • 場合1
        • ケース群とコントロール群について、観測したところ、「0,1,2」という量が観測された。今、ケース群・コントロール群という群分けと「0,1,2」という量の分布とに関係があるかどうかを考える。
        • ケース群・コントロール群という2群なので、群に対して順序を問題にする必要はないが、拡張することを考慮して、順序があるものとする。
          • ケースのAA,Aa,aaがn1=30,n2=40,n3=20人。コントロールがm1=40,m2=30,m3=20人だとする。
            • このデータは2カラム、1カラム目に90行、2カラム目に90行で、各レコードは「0,1,2」で記録される。
        • 「0,1,2」に順位をつけ、Kruskal-Wallis検定が可能。Rでの実行方法と結果は、次の通り。エクセルでは、D22,E22のカラム。自由度は1。
    • -
n1<-30
n2<-40
n3<-20
m1<-40
m2<-30
m3<-20
case<-c(rep(2,n1),rep(1,n2),rep(0,n3))
cont<-c(rep(2,m1),rep(1,m2),rep(0,m3))
kruskal.test(list(case,cont))
        • その結果は
        Kruskal-Wallis rank sum test

data:  list(case, cont) 
Kruskal-Wallis chi-squared = 1.1506, df = 1, p-value = 0.2834
        • このときの、分割表検定の結果(エクセル)は
2x3khai square	2.8571 	0.239651036
Chi-trend	0.952380952	0.329113996
Cochran-Armitage trend	0.947089947	0.330461144
Kruskal-Wallis	1.150596878	0.283424463
      • 場合2
        • AA群、Aa群、aa群の3群について、「ケース=1」「コントロール=0」の観測がなされる。ジェノタイプ3群という群分けと形質の「0,1」という量の分布とに関係があるかどうかを考える。
        • ジェノタイプ3群なので、群に対して順序を問題にしないか、問題にするかの選択が必要としてもよい。。
          • ケースのAA,Aa,aaがn1=30,n2=40,n3=20人。コントロールがm1=40,m2=30,m3=20人だとする。
            • このデータは3カラム、1カラム目に70行、2カラム目に70行、3カラム目に40行で、各レコードは「0,1」で記録される。
        • 「0,1」に順位をつけ、Kruskal-Wallis検定が可能。Rでの実行方法と結果は、次の通り。エクセルでは、D23,D23のカラム。自由度は2。
AA<-c(rep(n1,1),rep(m1,0))
Aa<-c(rep(n2,1),rep(m2,0))
aa<-c(rep(n3,1),rep(m3,0))

kruskal.test(list(AA,Aa,aa))
        Kruskal-Wallis rank sum test

data:  list(AA, Aa, aa) 
Kruskal-Wallis chi-squared = 2.8413, df = 2, p-value = 0.2416
        • Kruskal-Wallisは当然のことながら、ジェノタイプ3タイプについて、順序を考慮していないので、2x3分割表のカイ自乗検定(自由度2)と同じ挙動をする。
        • このときの、分割表検定の結果(エクセル)は
2x3khai square	2.8571 	0.239651036
χsq (1/2)=	1.1429 	0.285049646
χsq (11/12+22)=	2.3377 	0.126278968
χsq (11+12/22)=	0.0000 	1
Chi-trend	0.952380952	0.329113996
Cochran-Armitage trend	0.947089947	0.330461144
Kruskal-Wallis	2.841269841	0.241560596
      • 2x3分割表の各種検定とともにKruskal-Wallisの枠組みでの統計量とそのp値を算出するエクセルはこちら。(アップロード予定)
      • Kruskal-Wallis検定を別途行うエクセルはこちら
    • より面倒くさいデータにあてはめられるように、制約を課した手法を単純なデータにあてはめているので、すこしずつ保守的に寄ってくる。
    • Kruskal-Wallisの検算はRで。

ケースコントロール関連以外でパーミュテーションテスト(メモ)

plinkではQuantitative traitに対する単一SNP検定として、Waldテスト、尤度比検定を用いている。
こちらのレビュー(Crit Care. 2004; 8(3): 196–199.)の方が簡潔にまとまっている。
Statistics review 10: Further nonparametric methods
ちなみにこのレビューシリーズは、全14回(継続中?)
このシリーズ1−14回へのリンクはこちら

自前のパーミュテーションテストルーチンに入れるとして、参考となるサイトは。。。
こちら
ヨンクヒール・タプストラ検定Page's trend test・・・Kendall's tau
ヨンクヒール・タプストラ検定統計量は正規分布漸近近似
Rはこちら

均一確率分布からの最小P値の分布(多面体篇)



2006年12月5日の記事に、この件を書いた。そのときは、FWERの考え方と微分の考え方から、

Exp(min) = ¥int_{0}^{1} (¥alpha ¥times N ¥times (1-¥alpha)^{N-1}) d¥alpha

なる式を示した。

昨日もその別の考え方を書いた

もうひとつの考え方としては、こう。

N=2のとき、xy平面上の1辺の長さが1の正方形を考える。頂点を{0,0,0},{0,1,0},{1,0,0},{1,1,0}とする。このような正方形を底面として、{1,1,1}なる3次元空間の点を頂点とする錐を考える。最小P値の期待値は、この錐の体積である。これをN次元に拡張していくと、N=1のとき、『錐』に相当するのは、{0,0},{1,0},{1,1}を3頂点とする直角二等辺三角形であることがわかるし、多次元における、『最小P値考察用錐』というのは、{0,0,...,0},{1,0,...,0},{0,1,0,...,0},{0,0,1,0,...,0},...,{0,0,...0,1}が作るN-1次元立方体を底面とし、{1,1,1,...,1}を頂点とするN次元錐であって、その体積は¥frac{1}{N+1}である(らしい)ことがわかる。