集計
多数の因子があるとする。それらは、カテゴリカルな因子である。サンプルについて、それぞれの因子について観測されたデータが得られるとする。
何かの都合で、これらを複合カテゴリ(因子１のカテゴリがｘで因子２のそれがｙで・・・というのを弁別して個別のカテゴリとしたもの)ごとに集計したいとする。複数ＳＮＰのジェノタイプデータも、そのひとつ。
２分探索木を使って格納してみることにした。
リンクしているすべての関数につながっていないので、そのまま使うというわけには行かないので、主に自分のための覚書。
２分探索木に関する日記の記事は、こちらやこちら。
package binaryTree;
/**
* 2分探索木
*
* @version $Revision: 1.6 $, $Date: 2003/04/28 23:23:14 $
*/
import java.io.*;
public class BinarySearchTreeRYMultiCategory {
class Node {
Node left, right;
int id;
int[] category;
double count;
double freq;
double noderk;
double LNexpvalue;
double LNgeomean;
Node(Node left, Node right, int id,int[] category) {
this.left = left; this.right = right;
this.id=id;
this.category = category;
this.count=1;
}
Node(Node left, Node right, int id,int[] category,double freq) {
this.left = left; this.right = right;
this.id=id;
this.category = category;
this.count=freq;
}
Node(Node left, Node right) {
this.left = left; this.right = right;
this.id=-99;
this.category = new int[ncol];
for(int i=0;i<category.length;i++){
this.category[i]=-99;
}
this.count=0;
}
Node(Node left, Node right, int id) {
this.left = left; this.right = right;
this.id=id;
this.category = new int[ncol];
for(int i=0;i<category.length;i++){
this.category[i]=-99;
}
this.count=0;
}
}
public int ncol;
final private Node
zero = new Node(null, null,-999), // 各葉の子(番人用)
rootParent = new Node(null, zero); // {\tt rootParent.right} が根
private Node parent; private boolean isChildLeft; // 検索成功時の親
int numNode;
double N;
public boolean isEmpty() { return (rootParent.right == zero); }
public BinarySearchTreeRYMultiCategory(int ncol_){
ncol=ncol_;
}
public double r2;
public double rk;
public double chi2;
public Node searchNode(int id_) { // 検索（未登録なら {\tt null} を返す）
zero.id=-99; // {\tt zero} は番人
Node child = rootParent;
//int cmp = 1; // 最初の子は {\tt rootParent.right}
double cmp = 1;
boolean loop=true;
do {
parent = child;
if (cmp < 0) { child = child.left; isChildLeft = true; }
else { child = child.right; isChildLeft = false; }
//System.out.println("value="+value);
//System.out.println("child.value="+child.value);
if(child.id==-99){
loop=false;
}else{
cmp=-id_+child.id;
if(cmp==0){
loop=false;
}
}
//} while ((cmp = key.compareTo(child.key)) != 0);
} while (loop);
//System.out.println("childvalue="+child.value);
return (child == zero) ? null : child;
}
public Node insertNode(int id_,int[] category_) { // 挿入（登録）
final Node target = searchNode(id_);
if (target != null){ // すでに登録されている
target.count++;
return null;
}
Node newChild = new Node(zero, zero, id_,category_);
if (isChildLeft) parent.left = newChild;
else parent.right = newChild;
numNode++;
N++;
return newChild;
}
public Node insertNodeWithExp(int id_,int[] category_,double exp) { // 挿入（登録）
final Node target = searchNode(id_);
if (target != null){ // すでに登録されている
target.count++;
return null;
}
Node newChild = new Node(zero, zero, id_,category_);
newChild.LNexpvalue=exp;
newChild.LNgeomean=newChild.LNexpvalue/(double)newChild.category.length;
if (isChildLeft) parent.left = newChild;
else parent.right = newChild;
numNode++;
N++;
return newChild;
}
public Node insertNodeWithFreq(int id_,int[] category_,double freq) { // 挿入（登録）
final Node target = searchNode(id_);
if (target != null){ // すでに登録されている
target.count+=freq;
N+=freq;
return null;
}
Node newChild = new Node(zero, zero, id_,category_,freq);
if (isChildLeft) parent.left = newChild;
else parent.right = newChild;
numNode++;
N+=freq;
return newChild;
}
public Node insertNodeWithFreqExp(int id_,int[] category_,double freq,double exp) { // 挿入（登録）
final Node target = searchNode(id_);
if (target != null){ // すでに登録されている
target.count+=freq;
N+=freq;
return null;
}
Node newChild = new Node(zero, zero, id_,category_,freq);
newChild.LNexpvalue=exp;
newChild.LNgeomean=newChild.LNexpvalue/(double)newChild.category.length;
if (isChildLeft) parent.left = newChild;
else parent.right = newChild;
numNode++;
N+=freq;
return newChild;
}
public boolean deleteNode(int id_) { // 削除できれば {\tt true} を返す
final Node target = searchNode(id_);
if (target == null) return false; // 未登録
Node newChild; // {\tt target} の代わりに {\tt parent} の子になる
if (target.left == zero) newChild = target.right;
else if (target.right == zero) newChild = target.left;
else {
Node s = target.left; // {\tt s} を {\tt target} の次に小さいものに
if (s.right != zero) {
Node p; // {\tt s} の親
do { p = s; s = s.right; } while (s.right != zero);
p.right = s.left; s.left = target.left;
}
s.right = target.right; newChild = s;
}
if (isChildLeft) parent.left = newChild;
else parent.right = newChild;
numNode--;
N--;
return true; // 削除成功, C++ では {\tt delete target;} も必要
}
public boolean decreaseNode(int id_) { // 削除できれば {\tt true} を返す
final Node target = searchNode(id_);
if (target == null) return false; // 未登録
if (target.count>=1){
target.count--;
}else{
deleteNode(id_);
}
return true; // 削除成功, C++ では {\tt delete target;} も必要
}
public void calculate() {
System.out.println("#######");
rk=-1;
if (!isEmpty()) calculate(rootParent.right);
}
private void calculate(Node p) { // 深さ優先探索,中間順でキーを表示
if (p.left != zero) {
depth++; calculate(p.left); depth--;
}
if (p.right != zero) {
depth++; calculate(p.right); depth--;
}
//for (int i = 0; i < depth; i++) System.out.print("\t");
String st="ID\t";
st+=p.id+"\tdata\t";
for(int i=0;i<p.category.length;i++){
st+=p.category[i]+",";
}
st+="\tcount\t";
st+=p.count;
st+="\tfraction\t";
p.freq=p.count/N;
st+=p.freq;
st+="\texpFraction\t"+Math.exp(p.LNexpvalue);
double nodegeomean=Math.exp(p.LNgeomean);
st+="\tGeoMean\t"+nodegeomean;
p.noderk=(Math.log(p.freq)-p.LNgeomean)*(double)(p.category.length)/(double)(p.category.length-1);
p.noderk=Math.exp(p.noderk);
//System.out.println("plus\t"+p.noderk);
rk+=p.noderk;
st+="\trkPerType\t"+p.noderk;
System.out.println(st);
}
public void printNode() {
if (!isEmpty()) printNode(rootParent.right);
}
private int depth = 0;
private void printNode(Node p) { // 深さ優先探索,中間順でキーを表示
if (p.left != zero) {
depth++; printNode(p.left); depth--;
}
for (int i = 0; i < depth; i++) System.out.print("\t");
String st="";
st+=p.id+",";
for(int i=0;i<p.category.length;i++){
st+=p.category[i]+"";
}
st+=",";
//st+=p.count;
//st+=",";
double nodefreq=p.count/N;
st+=nodefreq;
//st+=","+Math.exp(p.LNexpvalue);
double nodegeomean=Math.exp(p.LNgeomean);
//st+=","+nodegeomean;
double noderk=(Math.log(nodefreq)-p.LNgeomean)*(double)(p.category.length)/(double)(p.category.length-1);
noderk=Math.exp(noderk);
//st+=","+noderk;
System.out.println(st);
if (p.right != zero) {
depth++; printNode(p.right); depth--;
}
}
private void giveExpAndGeoMean(){
}
public static void main(String[] args) throws IOException {
int n=10;
double[] data=new double[n];
for(int i=0;i<data.length;i++){
data[i]=Math.random();
}
int[][] d={{0,0,0,0,0,0},{0,1,0,1,0,1},{1,1,1,1,1,1},{0,0,0,0,0,0},{1,1,1,1,1,1}};
int[] id = {0,1,2,0,2};
BinarySearchTreeRYMultiCategory tree = new BinarySearchTreeRYMultiCategory(d[0].length);
/*
System.out.print("命令 Iabc: abcを挿入\n"
+ " Dabc: abcを削除\n"
+ " Sabc: abcを検索\n");
*/
String message;
for (int i=0;i<d.length;i++) {
int value=id[i];
if (tree.insertNode(value,d[i]) != null)
message = "登録しました.";
else message = "登録ずみです.";
//break;
System.out.println(message); tree.printNode();
System.out.println("numNode\t"+tree.numNode);
}
int[] searchID={2};
for(int i=0;i<searchID.length;i++){
if (tree.searchNode(id[searchID[i]]) != null)
message = "登録されています.";
else message = "登録されていません";
System.out.println(message);
}
int valueex=1;
if (tree.searchNode(valueex) != null)
message = "登録されています.";
else message = "登録されていません";
System.out.println(message);
int[] deleteID={0};
for(int i=0;i<deleteID.length;i++){
if (tree.deleteNode(id[deleteID[i]]))
message = "削除しました.";
else message = "登録されていません.";
System.out.println(message);
System.out.println("numNode\t"+tree.numNode);
}
tree.printNode();
/*
for (int i=0;i<data.length;i++) {
String message;
double value=data[i];
String insert="'I'";
switch ('I') {
case 'I':
if (tree.insertNode(value) != null)
message = "登録しました.";
else message = "登録ずみです.";
break;
case 'D':
if (tree.deleteNode(value))
message = "削除しました.";
else message = "登録されていません.";
break;
case 'S':
if (tree.searchNode(value) != null)
message = "登録されています.";
else message = "登録されていません";
break;
default:
message = "使えるのは I, D, S です.";
break;
}
System.out.println(message); tree.printNode();
}
*/
}
}package multiSNPLD;
public class MultiCategoricalChiPresentOnly {
int[][] data;
int nlin;
int ncol;
int[] ncat;
double[][] colfreq;
int[][][] count;
int[][][] typecount;
double N;
double rk;
binaryTree.BinarySearchTreeRYMultiCategory tr;
/**
* @param args
*/
public static void main(String[] args) {
// TODO 自動生成されたメソッド・スタブ
int[][] d0 ={{0,0},{0,1},{1,0},{1,1}};
int[][] d = {{0,0,0},{0,0,1},{0,1,0},{0,1,1},{1,0,0},{1,0,1},{1,1,0},{1,1,1}};
int[][] d2 ={{0,0,0,0,0,0},{1,1,1,1,1,1},{0,0,0,0,0,1}};
int[][] d3 ={
{0,0,0,0,0,0},
{1,1,1,1,1,1},
{0,0,0,1,1,1},
{0,0,1,1,1,1},
{0,1,1,1,1,1},
{0,0,0,0,1,1},
{0,0,0,0,0,1},
{1,1,1,1,1,2},
{1,1,1,1,2,2},
{1,1,1,2,2,2},
{1,1,2,2,2,2},
{1,2,2,2,2,2},
{0,0,0,0,0,1},
{0,0,0,0,1,1},
{0,0,0,1,1,1},
{0,0,1,1,1,1},
{0,1,1,1,1,1},
{2,2,2,2,2,2},
};
double[] freq=new double[d0.length];
for(int i=0;i<freq.length;i++){
freq[i]=Math.random();
}
StatUtilsZ.MiscUtilX.standard(freq);
MultiCategoricalChiPresentOnly mc=new MultiCategoricalChiPresentOnly(d0);
for(int i=0;i<mc.colfreq.length;i++){
for(int j=0;j<mc.colfreq[i].length;j++){
System.out.print(mc.colfreq[i][j]+"\t");
}
System.out.println();
}
System.out.println("mc\t"+mc.tr.rk);
System.out.println("============");
MultiCategoricalChiPresentOnly mc2=new MultiCategoricalChiPresentOnly(d0,freq);
for(int i=0;i<mc.colfreq.length;i++){
for(int j=0;j<mc.colfreq[i].length;j++){
System.out.print(mc.colfreq[i][j]+"\t");
}
System.out.println();
}
System.out.println("mc2\t"+mc2.tr.rk);
/*
for(int i=0;i<mc.count.length;i++){
if(mc.count[i][mc.ncol][0]!=0){
for(int j=0;j<mc.count[i].length;j++){
for(int k=0;k<mc.count[i][j].length;k++){
System.out.print(mc.count[i][j][k]);
}
}
System.out.println();
}
}
System.out.println("rk\t"+mc.rk);
System.out.println("chik\t"+mc.chik);
*/
}
public MultiCategoricalChiPresentOnly(int[][] d){
data=StatUtilsZ.MiscUtilX.DeepCopyInt2(d);
tr = new binaryTree.BinarySearchTreeRYMultiCategory(d[0].length);
nlin=data.length;
ncol=data[0].length;
int[] colmin=new int[ncol];
int[] colmax=new int[ncol];
//int[][] count=new int[ncol+1][1];
N=data.length;
//int nummultcat=1;
for(int i=0;i<data.length;i++){
for(int j=0;j<data[i].length;j++){
if(i==0){
colmin[j]=data[i][j];
colmax[j]=data[i][j];
}else{
if(colmin[j]>data[i][j]){
colmin[j]=data[i][j];
}
if(colmax[j]<data[i][j]){
colmax[j]=data[i][j];
}
}
}
}
colfreq=new double[ncol][0];
for(int i=0;i<colfreq.length;i++){
colfreq[i]=new double[colmax[i]-colmin[i]+1];
//nummultcat*=colfreq[i].length;
}
//System.out.println("nummultcat "+nummultcat);
//count=new int[nummultcat][ncol+1][1];
int[] kurai=new int[ncol];
kurai[0]=1;
for(int i=1;i<kurai.length;i++){
kurai[i]=kurai[i-1]*(colmax[i-1]-colmin[i-1]+1);
//System.out.println("kurai "+kurai[i]);
}
int[][] multtype=new int[data.length][ncol];
int[] id = new int[data.length];
for(int i=0;i<data.length;i++){
for(int j=0;j<data[i].length;j++){
colfreq[j][data[i][j]-colmin[j]]++;
multtype[i][j]=data[i][j]-colmin[j];
}
int serial=0;
for(int j=0;j<multtype[i].length;j++){
//System.out.println("kurai[j]*multtype[j] "+kurai[j]*multtype[j]);
serial+=kurai[j]*multtype[i][j];
}
id[i]=serial;
/*
//if(count[serial][ncol+1][0]==0){
//System.out.println("serial "+serial);
for(int j=0;j<count[serial].length-1;j++){
count[serial][j][0]=multtype[j];
}
//}
count[serial][ncol][0]++;
*/
}
for(int i=0;i<data.length;i++){
double Lnexp=0;
for(int j=0;j<multtype[i].length;j++){
double tmp = colfreq[j][multtype[i][j]]/(double)N;
System.out.println("tmp=\t"+tmp);
Lnexp+=Math.log((double)colfreq[j][multtype[i][j]]/(double)N);
}
tr.insertNodeWithExp(id[i], data[i],Lnexp);
//System.out.println("-------------");
//tr.printNode();
}
tr.printNode();
System.out.println("-------------");
tr.calculate();
/*
rk=-1;chik=0;
double k=1.0/(double)(ncol-1);
for(int i=0;i<count.length;i++){
if(count[i][ncol][0]!=0){
double elem = (1+k)*Math.log((double)count[i][ncol][0]/(double)N);
for(int j=0;j<count[i].length-1;j++){
elem-=k*Math.log((double)colfreq[j][count[i][j][0]]/(double)N);
}
rk+=Math.exp(elem);
}
}
chik=rk*N;
*/
}
public MultiCategoricalChiPresentOnly(int[][] d,double[] freq){
data=StatUtilsZ.MiscUtilX.DeepCopyInt2(d);
tr = new binaryTree.BinarySearchTreeRYMultiCategory(d[0].length);
nlin=data.length;
ncol=data[0].length;
int[] colmin=new int[ncol];
int[] colmax=new int[ncol];
//int[][] count=new int[ncol+1][1];
for(int i=0;i<freq.length;i++){
N+=freq[i];
}
//System.out.println("N="+N);
//N=data.length;
//int nummultcat=1;
for(int i=0;i<data.length;i++){
for(int j=0;j<data[i].length;j++){
if(i==0){
colmin[j]=data[i][j];
colmax[j]=data[i][j];
}else{
if(colmin[j]>data[i][j]){
colmin[j]=data[i][j];
}
if(colmax[j]<data[i][j]){
colmax[j]=data[i][j];
}
}
}
}
colfreq=new double[ncol][0];
for(int i=0;i<colfreq.length;i++){
colfreq[i]=new double[colmax[i]-colmin[i]+1];
//nummultcat*=colfreq[i].length;
}
//System.out.println("nummultcat "+nummultcat);
//count=new int[nummultcat][ncol+1][1];
int[] kurai=new int[ncol];
kurai[0]=1;
for(int i=1;i<kurai.length;i++){
kurai[i]=kurai[i-1]*(colmax[i-1]-colmin[i-1]+1);
//System.out.println("kurai "+kurai[i]);
}
int[][] multtype=new int[data.length][ncol];
int[] id = new int[data.length];
for(int i=0;i<data.length;i++){
for(int j=0;j<data[i].length;j++){
colfreq[j][data[i][j]-colmin[j]]+=freq[i];
multtype[i][j]=data[i][j]-colmin[j];
}
int serial=0;
for(int j=0;j<multtype[i].length;j++){
//System.out.println("kurai[j]*multtype[j] "+kurai[j]*multtype[j]);
serial+=kurai[j]*multtype[i][j];
}
id[i]=serial;
/*
//if(count[serial][ncol+1][0]==0){
//System.out.println("serial "+serial);
for(int j=0;j<count[serial].length-1;j++){
count[serial][j][0]=multtype[j];
}
//}
count[serial][ncol][0]++;
*/
}
for(int i=0;i<data.length;i++){
double Lnexp=0;
for(int j=0;j<multtype[i].length;j++){
double tmp = colfreq[j][multtype[i][j]]/(double)N;
//System.out.println("tmp=\t"+tmp);
Lnexp+=Math.log((double)colfreq[j][multtype[i][j]]/(double)N);
}
tr.insertNodeWithFreqExp(id[i], data[i],freq[i],Lnexp);
//System.out.println("-------------");
//tr.printNode();
}
tr.printNode();
System.out.println("-------------");
tr.calculate();
/*
rk=-1;chik=0;
double k=1.0/(double)(ncol-1);
for(int i=0;i<count.length;i++){
if(count[i][ncol][0]!=0){
double elem = (1+k)*Math.log((double)count[i][ncol][0]/(double)N);
for(int j=0;j<count[i].length-1;j++){
elem-=k*Math.log((double)colfreq[j][count[i][j][0]]/(double)N);
}
rk+=Math.exp(elem);
}
}
chik=rk*N;
*/
}
}
