2007-10-20

多重検定を重層的に繰り返したとき)(2)

この日の記事の問題点は
『統計量 $P_j=\prod_{i=1}^M p_{ij}$ を考えたこと』
独立事象について、生起確率を掛け合わせることと、P値(生起確率の累積)を掛け合わせることとは別物であるから。

以下の記述は、この点について問題があることに留意しつつ、積分その他については、メモとして使うので、残す。
この留意は以降、『複数の多テストスタディの累積』のシリーズの(６)、2007/10/26まで適用されるべきである。

Variable inflationがあるとき、累積化により増幅される。
補正するには、個々のセットでの補正によるy=x化がよい。
それができないときに、累積化後にy=xに乗るかどうかだけが知りたいときには、pをランクに変えてやれば、その用(詳細、略)には立つ。
累積にあたって、「片側的評価」をすることと、一部のテストが一部のセットのみで行われたときの、補完についての構成が今後必要。

import java.io.IOException;
import java.io.BufferedWriter;
import java.io.FileWriter;


public class MultSet {

	public int N;//No. tests per set
	public int M;//No. sets
	public double[][] p;//p ncol=M,nrow=N
	public double[] multP;
	public double[] expMultP;
	public double[] rankMultP;
	/**
	 * @param args
	 */
	public static void main(String[] args) throws IOException{
		// TODO 自動生成されたメソッド・スタブ
		int N=1000;
		int M=2;
		double[][] p= new double[N][M];
		BufferedWriter[] bw = new BufferedWriter[4];
		String out = "ChiUniform10";
		String outrank = out+"_rank";
		out+=".txt";
		outrank+=".txt";
		
		String corOut="cor"+out;
		String corOutRank="cor"+outrank;
		bw[0] = new BufferedWriter(new FileWriter(out));
		bw[1] = new BufferedWriter(new FileWriter(outrank));
		bw[2] = new BufferedWriter(new FileWriter(corOut));
		bw[3] = new BufferedWriter(new FileWriter(corOutRank));
		double[][] correctedP=new double[N][M];
		
		int x = 10000000;
		int seed = (int)(Math.random()*x);
		Utils.MersenneTwisterFast mz = new Utils.MersenneTwisterFast(seed);
		double lambda = 1;
		int df1=1;//自由度
		org.apache.commons.math.distribution.ChiSquaredDistributionImpl chidf1 = 
			new org.apache.commons.math.distribution.ChiSquaredDistributionImpl(df1);
		try{
			double uniformP=mz.nextDouble();
			uniformP=1;
			for(int i=0;i<p.length;i++){
				double tmp2 = mz.nextDouble();
				tmp2=0.5;
				double randomlambda = 1+(lambda-1)*2*tmp2;
				double tmpX = mz.nextDouble();
				for(int j=0;j<p[i].length;j++){
					double tmp = mz.nextDouble();
					//tmpX=tmp;
					tmp = chidf1.inverseCumulativeProbability(tmpX);
					//tmp*=lambda;
					tmp*=randomlambda;
					tmp = 1-chidf1.cumulativeProbability(tmp);
					p[i][j]=tmp;
					if(j>0){
						p[i][j]=uniformP;
					}
					//p[i][j]=Math.abs(mz.nextGaussian())*mz.nextDouble();
					//p[i][j]=Math.random();
				}
			}
			
			for(int i=0;i<M;i++){
				double[] tmp = new double[N];
				for(int j=0;j<tmp.length;j++){
					tmp[j]=p[j][i];
				}
				BBPPhenotypeVector3GC.gC gc = new BBPPhenotypeVector3GC.gC(tmp);
				for(int j=0;j<tmp.length;j++){
					correctedP[j][i]=1-chidf1.cumulativeProbability(chidf1.inverseCumulativeProbability(1-tmp[j])/lambda);
				}
			}
		} catch (org.apache.commons.math.MathException ex) {
			
		}
		
		MultSet ms = new MultSet(p);
		MultSet.PrintThomas(bw[0],ms.expMultP);
		ms.RankMultP();
		MultSet.PrintThomas(bw[1],ms.rankMultP);
		MultSet ms2 = new MultSet(correctedP);
		MultSet.PrintThomas(bw[2],ms2.expMultP);
		ms2.RankMultP();
		MultSet.PrintThomas(bw[3],ms2.rankMultP);

	}
	
	public MultSet(double[][] p_){
		N=p_.length;
		M=p_[0].length;
		p=StatUtilsX.MiscUtil.DeepCopyDouble2(p_);
		multP=new double[N];
		expMultP=new double[N];
		for(int i=0;i<N;i++){
			for(int j=0;j<M;j++){
				multP[i]+=Math.log(p[i][j]);
			}
			expMultP[i]=ExpMultPLog(multP[i],M);
		}
		
	}
	
	public static double ExpMultP(double p,int M){
		double ret=0;
		double sl[] = StatUtilsX.Fisher.serialLogfact(M);
		for(int i=0;i<M;i++){
			ret += Math.pow((-1),(i-1))/Math.exp(sl[i-1])*Math.pow(Math.log(p),(i-1));
		}
		ret *=p;
		return ret;
	}
	public static double ExpMultPLog(double logp,int M){
		double ret=0;
		double sl[] = StatUtilsX.Fisher.serialLogfact(M);
		for(int i=1;i<=M;i++){
			double tmp =Math.log(-logp);
			//System.out.println(tmp);
			//ret += Math.pow((-1),(i-1))/Math.exp(sl[i-1])*Math.pow(logp,(i-1));
			ret += 1/Math.exp(sl[i-1])*Math.exp(tmp*(i-1));
			
			//ret += Math.pow((-1),(i-1))/Math.exp(sl[i-1])*Math.pow(-1, i-1)*Math.exp(tmp*(i-1));
		}
		ret *=Math.exp(logp);
		return ret;
	}
	public static void PrintThomas(BufferedWriter bw,double[] p)throws IOException{
		double[] exp= StatUtilsX.MiscUtil.DeepCopyDouble1(p);
		StatUtilsX.MiscUtil.quickSortDouble(exp, 0, exp.length-1);
		for(int i=0;i<exp.length;i++){
			double tmp = (double)(i+1)/(double)(exp.length+1);
			BioBankPerm2.Run.out3File(bw, exp[i]+"\t"+tmp+"\n");
		}
		bw.close();
	}
	public void RankMultP(){
		double[][] rank = new double[N][M];
		for(int i=0;i<M;i++){
			double[] tmp = new double[N];
			for(int j=0;j<N;j++){
				tmp[j]=p[j][i];
			}
			int[][] tmprank=StatUtilsX.MiscUtil.bublSortRetInt2(tmp);
			for(int j=0;j<N;j++){
				//System.out.println("t\t"+tmprank[j][1]);
				
				rank[j][i]=(double)(tmprank[j][1]+1)/(double)(N+1);
			}
		}
		MultSet ms = new MultSet(rank);
		rankMultP = StatUtilsX.MiscUtil.DeepCopyDouble1(ms.expMultP);
	}

}