第7章 ノンパラメトリック組み合わせ解析の例 7.2 Multivariate paired Observation



  • 7.2 Multivariate paired Observation
    • 2群に属するn個体についてq種類の変数をk時刻に観測する、というようなデータ構造
    • 各時刻において、2群の値に大小の傾向がないことを検定
    • 教科書に記載されている個別時刻のP値が少しずれているように思えるので少々不安だが。。。このソースで、6章までと整合性がとれていることから、擬似乱数列のぶれによる値の乖離なのかもしれない。ちなみに、個別仮説1,2,3に対するPは0.15605 0.01785 0.01155 、Combining fucntiosとしてFisher,LiptakLogit,Tippettのぞれぞれにて

corrP Fisher 0.0114

corrPLiptakLogit 0.0114

corrP Tippett 0.0354


public class TensorDataVec {
public int rank;

public int[] dimension;

public double[][] datavec;
public double[] data;
public int[][] index;

public static void main(String[] args) {

int rank_= 4;
int[] dim_ = {3,10,2,1};
//double [] data_ = {8.5,8.6,9.8,9.4,14.9,14.1,6.1,6.2,7.7,7.3,8.0,7.4,9.8,9.6,11.5,11.0,13.0,12.9,
//14.5,14.1,16.7,15.9,19.0,18.6,6.1,6.3,6.4,6.5,7.7,7.9,5.1,4.9,8.7,7.6,9.1,8.7,6.4,6.7,9.2,8.9,13.3,12.6,
//7.9,7.6,7.7,8.0,10.2,10.5,5.1,5.1,6.3,6.5,9.6,9.2,10.7,10.2,12.8,12.1,14.7,13.8};
double[] data_ = {8.5,8.6,6.1,6.2,9.8,9.6,14.5,14.1,6.1,6.3,5.1,4.9,6.4,6.7,7.9,7.6,5.1,5.1,10.7,10.2,
9.8,9.4,7.7,7.3,11.5,11.0,16.7,15.9,6.4,6.5,8.7,7.6,9.2,8.9,7.7,8.0,6.3,6.5,12.8,12.1,
14.9,14.1,8.0,7.4,13.0,12.9,19.0,18.6,7.7,7.9,9.1,8.7,13.3,12.6,10.2,10.5,9.6,9.2,14.7,13.8};

TensorDataVec t = new TensorDataVec(rank_,dim_,data_);




//Create Tensor Object and print out its features
//System.out.println("Tensor object is consisted of int rank, int[] dimension, double[] data" +
//" and int[][] index");

//PrintTensor(t,"\t","\n");



/*
* 個々の仮説は、第1次元のそれぞれについて、第3次元をパーミュテーションしても
*/
double DEV= 0.00000001;
int perm =10000;
//異なるdelta値での統計量を算出
//double[] delta = {-1000,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,1000};
double[] delta = {0};
double[][][] permT = new double[t.dimension[0]][2][perm];
String mval = "";
String mval2 ="";
int c = 1000000;
int seed = (int)(Math.random()*c);
//int seed =100;
for(int m=0;m<delta.length;m++){
System.out.println("\n\n=================\nd value_" + m + "\n=================\n");
double[][] x = new double[t.dimension[0]][t.dimension[1]];
// int seed = 10000;
MersenneTwisterFast mt = new MersenneTwisterFast(seed);
MersenneTwisterFast mt2 = new MersenneTwisterFast(seed);
for(int i=0;i<x.length;i++){
for(int j=0;j<x[i].length;j++){
int[] a = {i,j,0,0};
int[] b = {i,j,1,0};

double aElem = ReturnElem(t,a);
double bElem = ReturnElem(t,b);
//System.out.println("aElem bElem delta " + aElem + " " + bElem + " " + delta[m]);
x[i][j] = aElem-bElem+delta[m];
//System.out.println("x " + x[i][j]);
}


}
double[][] to = new double[x.length][2];
//System.out.println("to.len " + to.length);
for(int i=0;i<to.length;i++){

to[i][0] = signPerm.sum(x[i]);
to[i][1] = signPerm.K(x[i]);

}


//複数の統計量を比較できるように
permT = signPerm.SumAndKAndWilcoxonSignPermVector(x,perm,mt);

//permutation P(個別仮説)の算出
double[] Pvalue = new double[t.dimension[0]];
for(int i=0;i<Pvalue.length;i++){
Pvalue[i] = signPerm.PermTestP(to[i][0],permT[i][0]);
}

//permutation P(個別仮説)をオリジナルデータと全パーミュテーション試行に算出
double[][] PAllPerm = new double[t.dimension[0]][perm];
for(int i=0;i<PAllPerm.length;i++){
PAllPerm[i] = signPerm.PermTestPPerm(permT[i][0]);
}





//PermTestP(to[0],permT[0]);
String StpermP = "PermP\t";
for(int i=0;i<Pvalue.length;i++){
StpermP += Pvalue[i] + "\t";
}
System.out.println(StpermP);

// T''算出
//Direct Combining Function

double tppDirect =0;
for(int i=0;i<to.length;i++){
tppDirect+=to[i][0];
}
double[] tppDirectPerm = new double[perm];
for(int i=0;i<perm;i++){
tppDirectPerm[i]=0;
for(int j=0;j<permT.length;j++){
tppDirectPerm[i]+=permT[j][0][i];
}
}
double corrPDirect = signPerm.PermTestP(tppDirect,tppDirectPerm);
//Fisher's omnibus combining function
double tppFisher = CombFxFisher(Pvalue);
double[] tppFisherPerm = new double[perm];
// LiptakLogit's omnibus combining function
double tppLiptakLogit = CombFxLiptakLogit(Pvalue);
double[] tppLiptakLogitPerm = new double[perm];
// Tippett's omnibus combining function
double tppTippett = CombFxTippett(Pvalue);
double[] tppTippettPerm = new double[perm];

for(int i=0;i<tppFisherPerm.length;i++){
double[] tmp = new double[PAllPerm.length];
for(int j=0;j<tmp.length;j++){
tmp[j]=PAllPerm[j][i];
}
tppFisherPerm[i] = CombFxFisher(tmp);
tppLiptakLogitPerm[i] = CombFxLiptakLogit(tmp);
tppTippettPerm[i] = CombFxTippett(tmp);
}

double corrPFisher = signPerm.PermTestP(tppFisher,tppFisherPerm);
System.out.println("corrP Fisher " + corrPFisher);
double corrPLiptakLogit = signPerm.PermTestP(tppLiptakLogit,tppLiptakLogitPerm);
System.out.println("corrPLiptakLogit " + corrPLiptakLogit);
double corrPTippett = signPerm.PermTestP(tppTippett,tppTippettPerm);
System.out.println("corrP Tippett " + corrPTippett);

System.out.println("corrP Direct " + corrPDirect);

/*
for(int i=0;i<permT.length;i++){
// 度数分布・信頼区間など
int[] h=signPerm.hist(permT[i][0],20);
double[][] hdouble = signPerm.histAccum(permT[i][0],1000);
// 統計量の信頼区間を算出表示
double confA = 0.000001;
double[] confInterval = signPerm.confInterval(hdouble,confA);

System.out.println("confident interval\t" + confInterval[0] + "\t" + confInterval[1]);

// 度数分布を標準出力
// System.out.println(histSt);
System.out.println("Histogram\n