Multiple correlation coefficient関連ソース

• Multiple correlation coefficientの記事はこちら
• Utilesパッケージ(自分の)に以下の６ソース
  
  • Calculator
  • InOutUtils
  • InvMat
  • MiscUtil
  • MultRegression
  • TestRun
• DiscreteMathToolsパッケージ(自分の)以下のソースも呼び出している
  
  • MatrixExec
• Utils.TestRunを実行するとひとまず、群馬大青木先生のページ(こちら)の例題が解かれる
• 以下、ソース
  
  • Utils.Calculator

package Utils;
import DiscreteMathTools.MatrixExec;
public class Calculator {

	public static double mean(double[] v){

		/*

		 * 平均は、データ列の０を中心とした１次のモーメント

		 */

		double ret = MiscUtil.momentX(v,0,1);

		

		return ret;

	}

	

	public static double[] meanV(double[][] m){

		double[] ret = new double[m[0].length];

		for(int i=0;i<ret.length;i++){

			double[] tmpv = new double[m.length];

			for(int j=0;j<tmpv.length;j++){

				tmpv[j]= m[j][i];

			}

			ret[i] = mean(tmpv);

			//System.out.println("mean " + mean[i]);

		}

		return ret;

	}

	public static double ss(double[] a){

		double ret =var(a)*a.length;

		return ret;

		

	}

	public static double var(double[] a){

		/*

		 * 分散は、データ列の平均を中心とした

		 * ２次のモーメント(２次のcentral moment)

		 */

		double mean = MiscUtil.momentX(a,0,1);

		double sigma = MiscUtil.momentX(a,mean,2);

		//double tmp =0;

		//for(int i=0;i<a.length;i++){

			//tmp += (a[i]-mean)*(a[i]-mean);

		//}

		//tmp /= a.length;

		//System.out.println("tmp " + tmp);

		//System.out.println("sigma " + sigma);

		return sigma;

	}

	public static double[] varV(double[][] a){

		double[] ret = new double[a[0].length];

		for(int i=0;i<ret.length;i++){

			double[] tmpv = new double[a.length];

			for(int j=0;j<ret.length;j++){

				tmpv[j]= a[j][i];

			}

			ret[i] = var(tmpv);

			//System.out.println("mean " + mean[i]);

		}

		return ret;

	}

	

	public static double max(double[] v){

		double ret = v[0];

		for(int i=0;i<v.length;i++){

			if(ret<v[i]){

				ret = v[i];

			}

		}

		return ret;

	}

	public static double min(double[] v){

		double ret = v[0];

		for(int i=0;i<v.length;i++){

			if(ret>v[i]){

				ret = v[i];

			}

		}

		return ret;

	}
	public static double[] sortUp(double[] v){

		double[] ret = MiscUtil.DeepCopyDouble1(v);

		MiscUtil.bublSortDouble(ret);

		return ret;

	}

	public static double[] sortDown(double[] v){

		double[] ret = sortUp(v);

		double[] ret2 = new double[ret.length];

		for(int i=0;i<ret2.length;i++){

			ret2[i]=ret[ret.length-i-1];

		}

		return ret2;

	}

	public static double[] mode(double[] v){

		double[] sort = sortUp(v);

		int numtie = 1;

		double[] modelist = {sort[0]};

		int counter=1;

		double current = sort[0];

		for(int i=1;i<sort.length;i++){

			System.out.println("sort[" + i +"] " + sort[i]);

			if(sort[i]==current){

				counter++;

				System.out.println("counter " + counter);

			}else{

				System.out.println("numtie " + numtie);

				if(counter>numtie){

					modelist = new double[1];

					modelist[0] = sort[i-1];

					numtie = counter;

				}else if(counter==numtie){

					modelist = MiscUtil.AddElemDouble1(modelist,sort[i-1]);

				}else{

					

				}

				counter=1;

				current=sort[i];

			}

		}

		return modelist;

	}

	public static double median(double[] v){

		double[] sort = sortUp(v);

		double ret = 0;

		if (v.length % 2 == 0){

			int i2 = v.length/2;

			int i1 = i2-1;

			ret = (sort[i1]+sort[i2])/2;

		}else{

			int i = v.length/2;

			ret = sort[i];

		}

		return ret;

	}

	public static double covar(double[] a, double[] b){

		double ret = 0;

		double meana = mean(a);

		double meanb = mean(b);

		for(int i=0;i<a.length;i++){

			ret += (a[i]-meana)*(b[i]-meanb);

		}

		

		return ret;

	}

	

	public static double[][] covarmatrix(double[][] m){

		double[] mean = meanV(m);

		/*

		double[] mean = new double[m[0].length];

		for(int i=0;i<mean.length;i++){

			double[] tmpv = new double[m.length];

			for(int j=0;j<tmpv.length;j++){

				tmpv[j]= m[j][i];

			}

			mean[i] = mean(tmpv);

			//System.out.println("mean " + mean[i]);

		}

		*/

		double[][] ret = new double[m[0].length][m[0].length];

		for(int i=0;i<m[0].length;i++){

			for(int j=0;j<m[0].length;j++){

				ret[i][j]=0;

				

				for(int k=0;k<m.length;k++){

					ret[i][j] += (m[k][i]-mean[i])*(m[k][j]-mean[j]);

				}

				

			}

		}

		return ret;

	}

	public static double[][] cov(double[] y, double[][] m){

		double meany = mean(y);

		double[] mean = new double[m[0].length];

		for(int i=0;i<mean.length;i++){

			double[] tmpv = new double[m.length];

			for(int j=0;j<tmpv.length;j++){

				tmpv[j]= m[j][i];

			}

			mean[i] = mean(tmpv);

			//System.out.println("mean " + mean[i]);

		}

		double[][] ret = new double[m[0].length][1];

		for(int i=0;i<m[0].length;i++){

			ret[i][0]=0;

			for(int j=0;j<y.length;j++){

				ret[i][0] += (y[j]-meany)*(m[j][i]-mean[i]);

			}

		}

		return ret;

	}

	public static void MultRegression(double[] y, double[][] mat){

		MultRegression mr = new MultRegression(y,mat);

		MultRegression.PrintMultRegression(mr,"\t","\n");

		

	}

	public static void MultRegressionPartial1(double[] y, double[][] mat,int[] x_){

		MultRegression mr = new MultRegression(y,mat,x_);

		MultRegression.PrintMultRegression(mr,"\t","\n");

		

	}

	public static void MultRegressionPartial2(double[] y, double[][] mat,int[] x_,int[] z_){

		MultRegression mr = new MultRegression(y,mat,x_,z_);

		MultRegression.PrintMultRegression(mr,"\t","\n");

		

	}

	public static double partialcorrcoef(double[] y, double[][] mat, int[] x,int[] z){

		MultRegression mrx = new MultRegression(y,mat,x);

		int[] xANDz = MiscUtil.DeepCopyInt1(x);

		for(int i=0;i<z.length;i++){

			xANDz = MiscUtil.AddElemInt1(xANDz, z[i]);

		}

		MultRegression.PrintMultRegression(mrx,"\t","\n");

		MultRegression mrxANDz = new MultRegression(y,mat,xANDz);

		MultRegression.PrintMultRegression(mrxANDz,"\t","\n");

		double ssResidX = mrxANDz.syy - mrxANDz.resid;

		double ssRegXbarZ = (mrxANDz.syy-mrxANDz.resid) - (mrx.syy-mrx.resid);

		double ret = ssRegXbarZ/ssResidX;

		return ret;

		

	}

}

• Utils.InOutUtils

package Utils;
import java.io.BufferedReader;

import java.io.BufferedWriter;

import java.io.FileReader;

import java.io.FileWriter;

import java.io.IOException;

import java.io.*;

import DiscreteMathTools.MiscUtil;
public class InOutUtils {
	public static void out3File(BufferedWriter bw1,String st){

		

		//BufferedWriter bw1 = null;

		try{

			//bw1 = new BufferedWriter(new FileWriter(file));

			

			

			bw1.write(st);

			//bw1.close();

	

		}catch(Exception e){

			

			System.out.println(e);

			

		}

		

		

	}

/*

	public static double[][] readDatFileIOString(String filePath)  throws IOException {

		String stFile = readUsingIOString(filePath);

		//System.out.println("read File " + stFile);

		String[] data1 = stFile.split("\n");

		System.out.println("samples " + data1.length);

		String[] phenotype = {"1","2"};

		double[] phenovalue = {-1,1};

		

		String[] genotype = {"1/1","1/2","2/2","0/0"};

		double[] genovalue = {1,0,-1,ConditionOP.UNKNOWN};

		double[][] data = new double[data1.length-1][0];

		int counter =0;

		for(int i=1;i<data1.length;i++){

			

			

			String[] tmp = data1[i].split(",");

			tmp[tmp.length-1] = tmp[tmp.length-1].trim();

			//for(int k=0;k<tmp.length;k++){

				//System.out.println(">>|" + tmp[k] + "|<<");

			//}

			data[counter] = new double[tmp.length-1];

			for(int j=1;j<tmp.length;j++){

				if(j==1){

					for(int k=0;k<phenotype.length;k++){

						//String plusN = phenotype[k] + "\n";

						if(tmp[j].equals(phenotype[k])){

							//System.out.println("jpheno IN " + j);

							data[counter][j-1]=phenovalue[k];

						//}else if(tmp[j].equals(plusN)){

							//data[counter][j-1]=phenovalue[k];

						}

					}

				}else{

					for(int k=0;k<genotype.length;k++){

						//String plusN = genotype[k] + "\n";

						if(tmp[j].equals(genotype[k])){

							//System.out.println("jgeno IN " + j);

							data[counter][j-1]=genovalue[k];

							//System.out.println("data[counter][j-1] " + data[counter][j-1]);

						//}else if(tmp[j].equals(plusN)){

							//data[counter][j-1]=genovalue[k];

						}

					}

				}

				

				

			}

			counter++;

		}

		int numsnp = data[0].length-1;

		//System.out.println("num SNPs " + numsnp);

		return data;

	}

*/

	public static double[][] readDatFileIOStringTAB(String filePath)  throws IOException {

		String stFile = readUsingIOString(filePath);

		//System.out.println("read File " + stFile);

		String[] data1 = stFile.split("\n");

		double[][] data = new double[data1.length][0];

		int counter =0;

		for(int i=0;i<data1.length;i++){

			

			

			String[] tmp = data1[i].split("\t");

			tmp[tmp.length-1] = tmp[tmp.length-1].trim();

			

			data[counter] = new double[tmp.length];

			for(int j=0;j<tmp.length;j++){

				data[i][j]=Double.parseDouble(tmp[j]);

			}

			counter++;

		}

		

		return data;

	}

/*

	public static double[][] readDatFile(String in){

		String[] phenotype = {"1","2"};

		double[] phenovalue = {-1,1};

		

		String[] genotype = {"1/1","1/2","2/2","0/0"};

		double[] genovalue = {1,0,-1,ConditionOP.UNKNOWN};

		

		double[][] ret = new double[0][0];

		String[][] records = new String[0][0];

		try{

			FileReader fr = new FileReader(in);

			

			BufferedReader objData=new BufferedReader(fr);

			

			

			int linecounter=0;

			while(objData.ready()){

				String[] line = objData.readLine().split(",");

				records=MiscUtil.AddElemString2(records,line);

			}

			fr.close();

			objData.close();

			System.out.println("records line " + records.length);

			ret = new double[records.length-1][0];

			int counter =0;

			for(int i=1;i<records.length;i++){

				//if(counter>0){

					ret[counter] = new double[records[i].length-1];

					for(int j=1;j<records[i].length;j++){

						if(j==1){

							for(int k=0;k<phenotype.length;k++){

								if(records[i][j].equals(phenotype[k])){

									ret[counter][j-1]=phenovalue[k];

									break;

								}

							}

							

						}else{

							for(int k=0;k<genotype.length;k++){

								if(records[i][j].equals(genotype[k])){

									ret[counter][j-1]=genovalue[k];

									break;

								}

							}

						}

					}

				//}

				counter ++;

			}

			

		}catch(IOException e){

			System.out.println(e);

		}

		

		

		

		

		return ret;

	}

*/

/*

	public static double[][] readDatFileCut(String in){

		String[] phenotype = {"1","2"};

		double[] phenovalue = {-1,1};

		

		String[] genotype = {"7","8","9","6"};

		double[] genovalue = {1,0,-1,ConditionOP.UNKNOWN};

		

		double[][] ret = new double[0][0];

		String[][] records = new String[0][0];

		try{

			FileReader fr = new FileReader(in);

			

			BufferedReader objData=new BufferedReader(fr);

			

			

			int linecounter=0;

			while(objData.ready()){

				String[] line = objData.readLine().split(",");

				records=MiscUtil.AddElemString2(records,line);

			}

			fr.close();

			objData.close();

			System.out.println("records line " + records.length);

			ret = new double[records.length][0];

			int counter =0;

			for(int i=0;i<records.length;i++){

				//if(counter>0){

					ret[counter] = new double[records[i].length];

					for(int j=0;j<records[i].length;j++){

						if(j==0){

							for(int k=0;k<phenotype.length;k++){

								if(records[i][j].equals(phenotype[k])){

									ret[counter][j]=phenovalue[k];

									break;

								}

							}

							

						}else{

							for(int k=0;k<genotype.length;k++){

								if(records[i][j].equals(genotype[k])){

									ret[counter][j]=genovalue[k];

									break;

								}

							}

						}

					}

				//}

				counter ++;

			}

			

		}catch(IOException e){

			System.out.println(e);

		}

		

		

		

		

		return ret;

	}

*/

	public static byte[] readUsingIO(String filePath) throws IOException {

		BufferedInputStream in = new BufferedInputStream(new FileInputStream(filePath));

		byte[] data = new byte[in.available()];

		in.read(data);

		int size = in.read(data);

		System.out.println("size="+size);

		/*

		int i = bis.available();

	      byte[] buf = new byte[i];

	      bis.read(buf);

	      */

	      System.out.println("取得した文字列は : " + new String(data));

	      System.out.println("--- 読み込み完了 ---");

		return data;

		

		

	}
	public static String readUsingIOString(String filePath) throws IOException {

		BufferedInputStream in = new BufferedInputStream(new FileInputStream(filePath));

		byte[] data = new byte[in.available()];

		in.read(data);

		int size = in.read(data);

		//System.out.println("size="+size);

		/*

		int i = bis.available();

	      byte[] buf = new byte[i];

	      bis.read(buf);

	      */

		String ret = new String(data);

	      //System.out.println("取得した文字列は : " + new String(data));

	      //System.out.println("--- 読み込み完了 ---");

		return ret;

		

		

	}

	
}

• Utils.InvMat

package Utils;
import java.text.*;
public class InvMat {

	//こちらのページの公開ソースの行列表現部分のみを演算対象部分のみの表現に修正したものです。

	//http://www.geocities.jp/java_sample_program/index.htm

	

	

	//pivotは、消去演算を行う前に、対象となる行を基準とし、それ以降の

	//行の中から枢軸要素の絶対値が最大となる行を見つけ出し、対象の行と

	//その行とを入れ替えることを行う関数である。

	void pivot(int k,double a[][],int N){

		double max,copy;

		//ipは絶対値最大となるk列の要素の存在する行を示す変数で、

		//とりあえずk行とする

		int ip=k;

		//k列の要素のうち絶対値最大のものを示す変数maxの値をとりあえず

		//max=|a[k][k]|とする

		max=Math.abs(a[k][k]);

		//k+1行以降、最後の行まで、|a[i][k]|の最大値とそれが存在する行を

		//調べる

		for(int i=k+1; i<N; i++){	

		  if(max<Math.abs(a[i][k])){

			ip=i;

			max=Math.abs(a[i][k]);

		  }

		}

		if(ip!=k){

		  for(int j=0; j<2*N; j++){

			//入れ替え作業

		  	copy    =a[ip][j];

			a[ip][j]=a[k][j];

			a[k][j] =copy;

		  }

		}

	}
	//ガウス・ジョルダン法により、消去演算を行う

	void sweep(int k,double a[][],int N){

	  double piv,mmm;

	  //枢軸要素をpivとおく

	  piv=a[k][k];

	  //k行の要素をすべてpivで割る a[k][k]=1となる

	  for(int j=0; j<2*N; j++)

		a[k][j]=a[k][j]/piv;

	  //	

	  for(int i=0; i<N; i++){

		mmm=a[i][k];

		//a[k][k]=1で、それ以外のk列要素は0となる

		//k行以外

		if(i!=k){

		  //i行において、k列から2N-1列まで行う	

		  for(int j=k; j<2*N; j++)

			a[i][j]=a[i][j]-mmm*a[k][j];

		}

	  }

	}

	

	//逆行列を求める演算

	void gyaku(double a[][],int N){

	    for(int k=0; k<N; k++){

		      pivot(k,a,N);

		      sweep(k,a,N);

		}

	}

	
	public static void main(String[] args) {

		//DecimalFormatクラスにより小数点以下5桁まで表示

		DecimalFormat ft=new DecimalFormat("##0.00000");

		

		double b[][] ={{1,-1,-1,1,-1,1,1,-1},

				{1,1,-1,-1,-1,-1,1,1},

				{1,-1,1,-1,-1,1,-1,1},

				{1,-1,-1,1,1,-1,-1,1},

				{1,1,1,1,-1,-1,-1,-1},

				{1,1,-1,-1,1,1,-1,-1},

				{1,-1,1,-1,1,-1,1,-1},

				{1,1,1,1,1,1,1,1}};

		double[][] invb = inv(b);

		String out ="元行列\n";

		for(int i=0;i<b.length;i++){

			for(int j=0;j<b[i].length;j++){

				out += b[i][j] + "\t";

			}

			out += "\n";

		}

		System.out.println(out);

		out="";

		for(int i=0;i<invb.length;i++){

			for(int j=0;j<invb[i].length;j++){

				out += invb[i][j] + "\t";

			}

			out += "\n";

		}

		System.out.println(out);

		

	}

	

	public static double[][] matPlusEmat(double[][] m){

		double[][] ret = new double[m.length][m[0].length*2];

		for(int i=0;i<ret.length;i++){

			for(int j=0;j<m[i].length;j++){

				ret[i][j]=m[i][j];

				if(i==j){

					ret[i][m[i].length+j]=1;

				}else{

					ret[i][m[i].length+j]=0;

				}

				

			}

			

		}

		return ret;

	}

	public static double[][] inv(double[][] b){

		double[][] ret = new double[b.length][b[0].length];

		double[][] a = matPlusEmat(b);

		int N=a.length;
		

   		//Gyakugyouretuクラスのインスタンスの生成

		InvMat g=new InvMat();

		

		//gyakuの呼び出し

		g.gyaku(a,N);

		

		

		for(int i=0; i<N; i++){

		  for (int j=N; j<2*N; j++){

			ret[i][j-N]=a[i][j];

		  	

		  }

		}

		

		return ret;

	}

	

}

• Utils.MiscUtil

/*

 * 作成日: 2006/05/16

 *

 * TODO この生成されたファイルのテンプレートを変更するには次へジャンプ:

 * ウィンドウ - 設定 - Java - コード・スタイル - コード・テンプレート

 */

package Utils;



/**

 * @author Ryo Yamada

 *

 * TODO この生成された型コメントのテンプレートを変更するには次へジャンプ:

 * ウィンドウ - 設定 - Java - コード・スタイル - コード・テンプレート

 */

public class MiscUtil {

	public static double DEV= 0.00000001;

	public static void main(String[] args) {

		double[] d = {4,-1,2,3,6,2,6,2,2};

		int[][] a = bublSortRetInt2(d);

		

		double[] c = {1,2,3,4,5,6,7,8,9,10};

		

		/*

		 * 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,3.0,3.0,3.0,4.0,4.0,

0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,

		 */

		/*

		double[] b1 = {10.0,10.0,10.0,10.0,11.0,11.0,11.0,11.0,11.0,12.0,12.0,12.0,12.0,12.0,13.0,13.0,13.0,13.0,13.0,14.0,14.0,14.0,14.0,14.0,15.0,15.0,15.0,15.0,15.0,16.0,16.0,16.0,16.0};

		double[] b2 = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0,1.0,1.0};

		int[][] ret = PermutationOP.CompareGradVector(b1,b2);

		String out ="";

		for(int i=0;i<ret.length;i++){

			out += b1[i] + " ";

			for(int j=0;j<ret[i].length;j++){

				out += ret[i][j] + " ";

				//System.out.println("ret " + ret[i][j] + " ");

			}

			out += "\n";

		}

		System.out.println(out);

		

		int[] x = {0,5,0,3};

		int[] y = {0,5,0,3};

		double[] rect1 =  OmPm.PermutationOP.GradiationRect(x,y);

		for(int i=0;i<rect1.length;i++){

			System.out.println("rect1 " + rect1[i]);

		}

		*/

		int k=10;

		double[] moment = new double[k];

		for(int i=0;i<k;i++){

			moment[i] = momentX(c,0,i);

			System.out.println("i moment " + i + " " + moment[i]);

		}

		double[] centralMoment = new double[k];

		for(int i=0;i<k;i++){

			centralMoment[i] = momentCentralX(c,i);

			System.out.println("i central moment " + i + " " + centralMoment[i]);

		}

		double[] centralMomentSt = new double[k];

		for(int i=0;i<k;i++){

			centralMomentSt[i] = momentCentralXSt(c,i);

			System.out.println("i central moment standardized " + i + " " + centralMomentSt[i]);

		}

		

	}

	

	public static void bublSort(int[] p){

		boolean flag=true;

		boolean sameflag=false;

		for(int i=p.length-1;i>=0;i--){

			for(int j=1;j<=i;j++){

				if(p[j-1]==p[j]){

					sameflag=true;

					//ret = 0;

					//return ret;

				}

				if(p[j-1]>p[j]){

					int tmp = p[j-1];

					p[j-1]=p[j];

					p[j]=tmp;

					//ret+=1;

					if(flag){

						flag=false;

					}else{

						flag=true;

					}

				}

			}

		}

	}

	public static void bublSortDouble(double[] p){

		boolean flag=true;

		boolean sameflag=false;

		for(int i=p.length-1;i>=0;i--){

			for(int j=1;j<=i;j++){

				if(p[j-1]==p[j]){

					sameflag=true;

					//ret = 0;

					//return ret;

				}

				if(p[j-1]>p[j]){

					double tmp = p[j-1];

					p[j-1]=p[j];

					p[j]=tmp;

					//ret+=1;

					if(flag){

						flag=false;

					}else{

						flag=true;

					}

				}

			}

		}

	}

	/*

	 * Deep Copies

	 * @author Ryo Yamada

	 */

	public static int[] DeepCopyInt1(int[] in){

		int[] ret=new int[in.length];

		for(int i=0;i<in.length;i++){

			ret[i]=in[i];

		}

		return ret;

	}

	public static int[][] DeepCopyInt2(int[][] in){

		int[][] ret=new int[in.length][0];

		for(int i=0;i<in.length;i++){

			ret[i] = DeepCopyInt1(in[i]);

		}

		return ret;

	}

	public static double[] DeepCopyDouble1(double[] in){

		double[] ret=new double[in.length];

		for(int i=0;i<in.length;i++){

			ret[i]=in[i];

		}

		return ret;

	}

	public static double[][] DeepCopyDouble2(double[][] in){

		double[][] ret=new double[in.length][0];

		for(int i=0;i<in.length;i++){

			ret[i] = DeepCopyDouble1(in[i]);

		}

		return ret;

	}

	public static int[] AddElemInt1(int[] in,int e){

		int[] ret = new int[in.length+1];

		for(int i=0;i<in.length;i++){

			ret[i]=in[i];

		}

		ret[in.length]=e;

		return ret;

	}

	public static int[][] AddElemInt2(int[][] in, int[] e){

		int[][] ret=new int[in.length+1][0];

		for(int i=0;i<in.length;i++){

			ret[i]=DeepCopyInt1(in[i]);

			

		}

		ret[in.length] = DeepCopyInt1(e);

		return ret;

	}

	public static double[] AddElemDouble1(double[] in,double e){

		double[] ret = new double[in.length+1];

		for(int i=0;i<in.length;i++){

			ret[i]=in[i];

		}

		ret[in.length]=e;

		return ret;

	}

	public static double[][] AddElemDouule2(double[][] in, double[] e){

		double[][] ret=new double[in.length+1][0];

		for(int i=0;i<in.length;i++){

			ret[i]=DeepCopyDouble1(in[i]);

			

		}

		ret[in.length] = DeepCopyDouble1(e);

		return ret;

	}

	public static double[][] AddElemDouble3(double[][] in, double[][] e){

		double[][] ret=new double[in.length+e.length][0];

		for(int i=0;i<in.length;i++){

			ret[i]=DeepCopyDouble1(in[i]);

			

		}

		for(int i=0;i<e.length;i++){

			ret[in.length+i] = DeepCopyDouble1(e[i]);

		}

		

		return ret;

	}

	public static String[] DeepCopyString1(String[] in){

		String[] ret=new String[in.length];

		for(int i=0;i<in.length;i++){

			ret[i]=in[i];

		}

		return ret;

	}

	public static String[][] DeepCopyString2(String[][] in){

		String[][] ret=new String[in.length][0];

		for(int i=0;i<in.length;i++){

			ret[i] = DeepCopyString1(in[i]);

		}

		return ret;

	}

	public static String[] AddElemString1(String[] in,String e){

		String[] ret = new String[in.length+1];

		for(int i=0;i<in.length;i++){

			ret[i]=in[i];

		}

		ret[in.length]=e;

		return ret;

	}

	public static String[][] AddElemString2(String[][] in, String[] e){

		String[][] ret=new String[in.length+1][0];

		for(int i=0;i<in.length;i++){

			ret[i]=DeepCopyString1(in[i]);

			

		}

		ret[in.length] = DeepCopyString1(e);

		return ret;

	}

	/*

	 * Return order 

	 */

	public static int[] bublSortRetInt(double[] p){

		double[] copyP = DeepCopyDouble1(p);

		int[] ret = new int[p.length];

		for(int i=0;i<ret.length;i++){

			ret[i]=i;

		}

		boolean flag=true;

		boolean sameflag=false;

		for(int i=copyP.length-1;i>=0;i--){

			for(int j=1;j<=i;j++){

				if(copyP[j-1]==copyP[j]){

					sameflag=true;

					//ret = 0;

					//return ret;

				}

				if(copyP[j-1]>copyP[j]){

					double tmp = copyP[j-1];

					int tmpint = ret[j-1];

					copyP[j-1]=copyP[j];

					ret[j-1]=ret[j];

					copyP[j]=tmp;

					ret[j]=tmpint;

					//ret+=1;

					if(flag){

						flag=false;

					}else{

						flag=true;

					}

				}

			}

		}

		/*

		for(int i=0;i<copyP.length;i++){

			System.out.println(copyP[i]);

			System.out.println(ret[i]);

		}

		*/

		int[] ret2 = new int[ret.length];

		for(int i=0;i<ret.length;i++){

			ret2[ret[i]] = i;

		}

		return ret2;

	}

	 

	/*

	 * arg : list of double values, some of which are identical

	 * ret : (1) ascending order.  Different number for multiple values

	 *       (2) ascending order.  Same number for multiple values

	 *       (3) descending order. Same number for multiple values

	 *       (4) number of records of the same value with self.

	 *       (5) more deviated. Smaller value between (2) + (4) and (3) + (4)

	 */

	

	public static int[][] bublSortRetInt2(double[] p){

		double[] copyP = DeepCopyDouble1(p);

		int[] ret = new int[p.length];

		for(int i=0;i<ret.length;i++){

			ret[i]=i;

		}

		boolean flag=true;

		boolean sameflag=false;

		for(int i=copyP.length-1;i>=0;i--){

			for(int j=1;j<=i;j++){

				if(Math.abs(copyP[j-1]-copyP[j])<DEV){

				//if(copyP[j-1]==copyP[j]){

					sameflag=true;

				}

				if(copyP[j-1]>copyP[j]){

					double tmp = copyP[j-1];

					int tmpint = ret[j-1];

					copyP[j-1]=copyP[j];

					ret[j-1]=ret[j];

					copyP[j]=tmp;

					ret[j]=tmpint;

					if(flag){

						flag=false;

					}else{

						flag=true;

					}

				}

			}

		}

		

				

		int[] retorder = new int[ret.length];

		int[] sameVsameO = new int[ret.length];

		int[] revorder = new int[ret.length];

		int[] numsame = new int[ret.length];

		

		int[][] retset = new int[ret.length][5];

		

		

		

		for(int i=0;i<ret.length;i++){

			retset[ret[i]][0] = i;

			retorder[ret[i]] = i;

		}

		

		

		

		

		for(int i=0;i<retorder.length;i++){

			if(i>0){

				if(Math.abs(p[ret[i]]-p[ret[i-1]])<DEV){

				//if(p[ret[i]]==p[ret[i-1]]){

					numsame[ret[i]]=numsame[ret[i-1]];

					retset[ret[i]][3]=retset[ret[i-1]][3];

					sameVsameO[ret[i]]=sameVsameO[ret[i-1]];

					retset[ret[i]][1]=retset[ret[i-1]][1];

				}else{

					

					numsame[ret[i]]=1;

					retset[ret[i]][3]=1;

					sameVsameO[ret[i]]=retorder[ret[i]];

					retset[ret[i]][1]=retset[ret[i]][0];

					

					for(int j=i+1;j<ret.length;j++){

						if(Math.abs(p[ret[j]]-p[ret[j-1]])<DEV){

						//if(p[ret[j]]==p[ret[j-1]]){
							numsame[ret[i]]++;

							retset[ret[i]][3]++;

						}else{

							break;

						}

					}

				}

			}else{
				numsame[ret[i]]=1;

				retset[ret[i]][3]=1;

				sameVsameO[ret[i]]=retorder[ret[i]];

				retset[ret[i]][1]=retset[ret[i]][1];

				for(int j=i+1;j<ret.length;j++){

					if(Math.abs(p[ret[j]]-p[ret[j-1]])<DEV){

					//if(p[ret[j]]==p[ret[j-1]]){
						numsame[ret[i]]++;

						retset[ret[i]][3]++;

					}else{

						break;

					}

				}

				

			}

		}

		

		

		

		
		for(int i=0;i<ret.length;i++){

			revorder[i] = ret.length-(sameVsameO[i]+numsame[i]);

			retset[i][2] = ret.length-(retset[i][1]+retset[i][3]);

			if(retset[i][1]<retset[i][2]){

				retset[i][4]=retset[i][1]+retset[i][3];

			}else{

				retset[i][4]=retset[i][2]+retset[i][3];

			}

		}

		

		/*

		System.out.println("Input ");

		for(int i=0;i<p.length;i++){

			System.out.println("i=" + i  + "  "+ p[i]);

		}

		System.out.println("Ascending-ordered list ");

		for(int i=0;i<ret.length;i++){

			System.out.println("i=" + i + "  " + ret[i] + " " + p[ret[i]]);

		}

		System.out.println("(1) ascending order.  Different number for multiple values ");

		for(int i=0;i<retorder.length;i++){

			System.out.println("i=" + i + "  "+ retorder[i]);

		}

		System.out.println("(2) ascending order.  Same number for multiple values ");

		for(int i=0;i<sameVsameO.length;i++){

			System.out.println("i=" + i + "  "+ sameVsameO[i]);

		}

		System.out.println("(3) descending order. Same number for multiple values");

		for(int i=0;i<revorder.length;i++){

			System.out.println("i=" + i + "  "+ revorder[i]);

		}

		System.out.println("(4) number of records of the same value with self.");

		for(int i=0;i<numsame.length;i++){

			System.out.println("i=" + i + "  "+ numsame[i]);

		}

		*/

		//int[][] ret2 = new int[0][0];

		return retset;

	}
	/*

	 * arg : list of double values, some of which are identical

	 * ret : (1) ascending order.  Different number for multiple values

	 *       (2) ascending order.  Same number for multiple values

	 *       (3) descending order. Same number for multiple values

	 *       (4) number of records of the same value with self.

	 */

	

	public static int[][] bublSortSameVals(double[] p){

		double[] copyP = DeepCopyDouble1(p);

		int[] ret = new int[p.length];

		for(int i=0;i<ret.length;i++){

			ret[i]=i;

		}

		boolean flag=true;

		boolean sameflag=false;

		for(int i=copyP.length-1;i>=0;i--){

			for(int j=1;j<=i;j++){

				if(copyP[j-1]==copyP[j]){

					sameflag=true;

				}

				if(copyP[j-1]>copyP[j]){

					double tmp = copyP[j-1];

					int tmpint = ret[j-1];

					copyP[j-1]=copyP[j];

					ret[j-1]=ret[j];

					copyP[j]=tmp;

					ret[j]=tmpint;

					if(flag){

						flag=false;

					}else{

						flag=true;

					}

				}

			}

		}

		

		

		

		int[] retorder = new int[ret.length];

		for(int i=0;i<ret.length;i++){

			retorder[ret[i]] = i;

		}

		

		

		int[] numsame = new int[ret.length];

		int[] sameVsameO = new int[ret.length];

		

		for(int i=0;i<retorder.length;i++){

			if(i>0){

				if(p[ret[i]]==p[ret[i-1]]){

					numsame[ret[i]]=numsame[ret[i-1]];

					sameVsameO[ret[i]]=sameVsameO[ret[i-1]];

				}else{

					

					numsame[ret[i]]=1;

					sameVsameO[ret[i]]=retorder[ret[i]];

					for(int j=i+1;j<ret.length;j++){

						if(p[ret[j]]==p[ret[j-1]]){

	

							numsame[ret[i]]++;

						}else{

							break;

						}

					}

				}

			}else{

	

				numsame[ret[i]]=1;

				sameVsameO[ret[i]]=retorder[ret[i]];

				for(int j=i+1;j<ret.length;j++){

					if(p[ret[j]]==p[ret[j-1]]){

	

						numsame[ret[i]]++;

					}else{

						break;

					}

				}

				

			}

		}

		

		

		

		

		int[] revorder = new int[ret.length];

		for(int i=0;i<ret.length;i++){

			revorder[i] = ret.length-(sameVsameO[i]+numsame[i]);

		}

		

		System.out.println("Input ");

		for(int i=0;i<p.length;i++){

			System.out.println("i=" + i  + "  "+ p[i]);

		}

		System.out.println("Ascending-ordered list ");

		for(int i=0;i<ret.length;i++){

			System.out.println("i=" + i + "  " + ret[i] + " " + p[ret[i]]);

		}

		System.out.println("(1) ascending order.  Different number for multiple values ");

		for(int i=0;i<retorder.length;i++){

			System.out.println("i=" + i + "  "+ retorder[i]);

		}

		System.out.println("(2) ascending order.  Same number for multiple values ");

		for(int i=0;i<sameVsameO.length;i++){

			System.out.println("i=" + i + "  "+ sameVsameO[i]);

		}

		System.out.println("(3) descending order. Same number for multiple values");

		for(int i=0;i<revorder.length;i++){

			System.out.println("i=" + i + "  "+ revorder[i]);

		}

		System.out.println("(4) number of records of the same value with self.");

		for(int i=0;i<numsame.length;i++){

			System.out.println("i=" + i + "  "+ numsame[i]);

		}

		int[][] ret2 = new int[0][0];

		return ret2;

	}

	public static boolean IdenticalDouble(double a, double b, double DEV){

		boolean ret = false;

		if(Math.abs(a-b)<DEV){

			ret = true;

		}

		return ret;

	}

	public static double momentCentralXSt(double[] a,int dim){

		

		double mean = momentX(a,0,1);

		double sigma = momentX(a,mean,2);

		double ret = momentX(a,mean,dim);

		sigma = Math.sqrt(sigma);

		double tmp = Math.log(sigma)*dim;

		ret /=Math.exp(tmp);

		return ret;

	}

	public static double momentCentralX(double[] a,int dim){

		double mean = momentX(a,0,1);

		double sigma = momentX(a,mean,2);

		sigma = Math.sqrt(sigma);

		//System.out.println("mean " + mean);

		double ret = momentX(a,mean,dim);

		return ret;

	}

	public static double momentX(double[] a,double m,int dim){

		double ret = 0;

		for(int i=0;i<a.length;i++){

			double tmp = 0;

			double pm =1;

			if(a[i] != m){

				if(a[i]>m){

					//System.out.println("a[i] " + a[i]);

					double tmp2 = Math.log(a[i]-m);

					//System.out.println("tmp2 " + tmp2);

					tmp = tmp2 * dim;

					

				}else{

					//System.out.println("a[i] " + a[i]);

					double pm2 =-1;

					double tmp2 = Math.log(-(a[i]-m));

					//System.out.println("tmp2 " + tmp2);

					tmp = tmp2 * dim;

					for(int j=0;j<dim;j++){

						//tmp+=tmp2;

						pm*=pm2;

					}

					//System.out.println("pm " + pm2);

				}

//				System.out.println("tmp " + tmp);

				double ttt = Math.exp(tmp)*pm;

				//System.out.println("ttt " + ttt);

				ret += Math.exp(tmp)*pm;

			}

			

			

			

		}

		ret /= a.length;

		return ret;

	}
	public static double momentCentralXStUNK(double[] a,int dim,double UNK){

		

		double mean = momentXUNK(a,0,1,UNK);

		double sigma = momentXUNK(a,mean,2,UNK);

		double ret = momentXUNK(a,mean,dim,UNK);

		sigma = Math.sqrt(sigma);

		double tmp = Math.log(sigma)*dim;

		ret /=Math.exp(tmp);

		return ret;

	}
	public static double momentCentralXUNK(double[] a,int dim,double UNK){

		double mean = momentXUNK(a,0,1,UNK);

		double sigma = momentXUNK(a,mean,2,UNK);

		sigma = Math.sqrt(sigma);

		//System.out.println("mean " + mean);

		double ret = momentXUNK(a,mean,dim,UNK);

		return ret;

	}
	public static double momentXUNK(double[] a,double m,int dim,double UNK){

		double ret = 0;

		int counter = 0;

		for(int i=0;i<a.length;i++){

			if(a[i] != UNK){

				counter++;

				double tmp = 0;

				double pm =1;

				if(a[i] != m){

					if(a[i]>m){

						//System.out.println("a[i] " + a[i]);

						double tmp2 = Math.log(a[i]-m);

						//System.out.println("tmp2 " + tmp2);

						tmp = tmp2 * dim;

						

					}else{

						//System.out.println("a[i] " + a[i]);

						double pm2 =-1;

						double tmp2 = Math.log(-(a[i]-m));

						//System.out.println("tmp2 " + tmp2);

						tmp = tmp2 * dim;

						for(int j=0;j<dim;j++){

							//tmp+=tmp2;

							pm*=pm2;

						}

						//System.out.println("pm " + pm2);

					}

					ret += Math.exp(tmp)*pm;

				}

				//System.out.println("tmp " + tmp);

				//double ttt = Math.exp(tmp)*pm;

				//System.out.println("ttt " + ttt);

				

			}

			

			

			

		}

		

		//ret /= a.length;

		ret /= counter;

		return ret;

	}

	

	/*

	 * Add 1 column to mat

	 */

	public static double[][] addDummyCol(double[][] mat){

		double[][] mat2 = new double[mat.length][mat[0].length+1];

		for(int i=0;i<mat2.length;i++){

			mat2[i][0]=1;

			for(int j=0;j<mat[i].length;j++){

				mat2[i][j+1]=mat[i][j];

			}

		}

		return mat2;

	}
}

• Utils.MultRegression

package Utils;
import DiscreteMathTools.MatrixExec;
public class MultRegression {

	double meany;

	double vary;

	double[] meanx;

	double[] varx;

	double[][] covmat;

	double[][] cov;

	double[][] invcovmat;

	double[][] b;

	double b0;

	double resid;

	double reg;

	double[] y;

	double[] yhat;

	double syy;

	double[][] ssx;

	double  multcorrcoef;

	double[] meansquare;

	double anovaF;

	int[] anovaDF;

	int[] x;//基礎回帰対象の変数ID

	int[] z;//追加回帰対象の変数ID

	int[] xANDz;//基礎回帰もしくは追加回帰対象とする独立変数のID

	

	public MultRegression(double[] y_, double[][] mat){

		y = MiscUtil.DeepCopyDouble1(y_);

		syy = Calculator.ss(y);

		covmat = Calculator.covarmatrix(mat);

		//MatrixExec.PrintMat(covmat,"\t","\n");

		

		cov = Calculator.cov(y,mat);

		//MatrixExec.PrintMat(cov,"\t","\n");

		

		invcovmat = InvMat.inv(covmat);

		//MatrixExec.PrintMat(invcovmat, "\t","\n");

				

		

		

		b = MatrixExec.Prodmat(invcovmat,cov);

		//MatrixExec.PrintMat(b,"\t","\n");

		meany = Calculator.mean(y);

		double[] meanx = Calculator.meanV(mat);

		vary = Calculator.var(y);

		double[][] meanXmat = MatrixExec.Vec2Mat(meanx);

		double[][] bt=MatrixExec.Transversmat(b);

		double[][] tmp = MatrixExec.Prodmat(bt,meanXmat);

		b0 = meany - tmp[0][0];

		yhat = new double[y.length];

		for(int i=0;i<y.length;i++){

			yhat[i] = b0;

			for(int j=0;j<b.length;j++){

				yhat[i] += b[j][0] * mat[i][j];

			}

		}

		resid = 0;

		for(int i=0;i<y.length;i++){

			resid += (y[i]-yhat[i])*(y[i]-yhat[i]);

		}

		reg = syy-resid;

		multcorrcoef = 1- resid/syy;

		meansquare = new double[3];

		meansquare[0] = (syy-resid)/((double)(mat[0].length));

		meansquare[1] = resid/((double)(y.length-1-mat[0].length));

		meansquare[2] = syy/((double)(y.length-1));

		anovaF = meansquare[0]/meansquare[1];

		anovaDF = new int[2];

		anovaDF[0] = mat[0].length;

		anovaDF[1] = y.length-1-mat[0].length;

		x = new int[mat[0].length];

		xANDz = new int[mat[0].length];

		for(int i=0;i<x.length;i++){

			x[i]=i;

			xANDz[i]=i;

		}

		z = new int[0];

		

	}

	public MultRegression(double[] y_, double[][] mat,int[] x_){

		y = MiscUtil.DeepCopyDouble1(y_);

		double[][] mat2 = new double[mat.length][x_.length];

		for(int i=0;i<mat2.length;i++){

			for(int j=0;j<mat2[i].length;j++){

				mat2[i][j]=mat[i][x_[j]];

			}

		}

		syy = Calculator.ss(y);

		covmat = Calculator.covarmatrix(mat2);

		//MatrixExec.PrintMat(covmat,"\t","\n");

		

		cov = Calculator.cov(y,mat2);

		//MatrixExec.PrintMat(cov,"\t","\n");

		

		invcovmat = InvMat.inv(covmat);

		//MatrixExec.PrintMat(invcovmat, "\t","\n");

				

		

		

		b = MatrixExec.Prodmat(invcovmat,cov);

		//MatrixExec.PrintMat(b,"\t","\n");

		meany = Calculator.mean(y);

		double[] meanx = Calculator.meanV(mat2);

		vary = Calculator.var(y);

		double[][] meanXmat = MatrixExec.Vec2Mat(meanx);

		double[][] bt=MatrixExec.Transversmat(b);

		double[][] tmp = MatrixExec.Prodmat(bt,meanXmat);

		b0 = meany - tmp[0][0];

		yhat = new double[y.length];

		for(int i=0;i<y.length;i++){

			yhat[i] = b0;

			for(int j=0;j<b.length;j++){

				yhat[i] += b[j][0] * mat[i][j];

			}

		}

		resid = 0;

		for(int i=0;i<y.length;i++){

			resid += (y[i]-yhat[i])*(y[i]-yhat[i]);

		}

		reg = syy-resid;

		multcorrcoef = 1- resid/syy;

		meansquare = new double[3];

		meansquare[0] = (syy-resid)/((double)(mat2[0].length));

		meansquare[1] = resid/((double)(y.length-1-mat2[0].length));

		meansquare[2] = syy/((double)(y.length-1));

		anovaF = meansquare[0]/meansquare[1];

		anovaDF = new int[2];

		anovaDF[0] = mat[0].length;

		anovaDF[1] = y.length-1-mat[0].length;

		x = MiscUtil.DeepCopyInt1(x_);

		xANDz = MiscUtil.DeepCopyInt1(x_);

		

		z = new int[0];

		

	}

	public MultRegression(double[] y_, double[][] mat,int[] x_,int[] z_){

		y = MiscUtil.DeepCopyDouble1(y_);

		syy = Calculator.ss(y);

		covmat = Calculator.covarmatrix(mat);

		//MatrixExec.PrintMat(covmat,"\t","\n");

		

		cov = Calculator.cov(y,mat);

		//MatrixExec.PrintMat(cov,"\t","\n");

		

		invcovmat = InvMat.inv(covmat);

		//MatrixExec.PrintMat(invcovmat, "\t","\n");

				

		

		

		b = MatrixExec.Prodmat(invcovmat,cov);

		//MatrixExec.PrintMat(b,"\t","\n");

		meany = Calculator.mean(y);

		double[] meanx = Calculator.meanV(mat);

		vary = Calculator.var(y);

		double[][] meanXmat = MatrixExec.Vec2Mat(meanx);

		double[][] bt=MatrixExec.Transversmat(b);

		double[][] tmp = MatrixExec.Prodmat(bt,meanXmat);

		b0 = meany - tmp[0][0];

		yhat = new double[y.length];

		for(int i=0;i<y.length;i++){

			yhat[i] = b0;

			for(int j=0;j<b.length;j++){

				yhat[i] += b[j][0] * mat[i][j];

			}

		}

		resid = 0;

		for(int i=0;i<y.length;i++){

			resid += (y[i]-yhat[i])*(y[i]-yhat[i]);

		}

		reg = syy-resid;

		multcorrcoef = 1- resid/syy;

		meansquare = new double[3];

		meansquare[0] = (syy-resid)/((double)(mat[0].length));

		meansquare[1] = resid/((double)(y.length-1-mat[0].length));

		meansquare[2] = syy/((double)(y.length-1));

		anovaF = meansquare[0]/meansquare[1];

		anovaDF = new int[2];

		anovaDF[0] = mat[0].length;

		anovaDF[1] = y.length-1-mat[0].length;

		x = MiscUtil.DeepCopyInt1(x_);

		z = MiscUtil.DeepCopyInt1(z_);

		xANDz = MiscUtil.DeepCopyInt1(x_);

		for(int i=0;i<z.length;i++){

			xANDz = MiscUtil.AddElemInt1(xANDz, z[i]);

		}

		

	}

	

	

	

	public static void PrintMultRegression(MultRegression mr,String sep1,String sep2){

		String ret = StringMultRegression(mr,sep1,sep2);

		System.out.println(ret);

	}

	public static String StringMultRegression(MultRegression mr, String sep1, String sep2){

		String ret="";

		ret += "covmat\n";

		ret += MatrixExec.StringMat(mr.covmat, "\t", "\n");

		ret += "Eauation\n";

		ret += "Y_i = " + mr.b0 +" + ";

		for(int i=0;i<mr.b.length-1;i++){

			int xx = i+1;

			ret += mr.b[i][0] + " * X_i"+xx+" + " ;

		}

		ret += mr.b[mr.b.length-1][0] + " * X_i" + mr.b.length + "\n";

		

		ret += "y and yhat dif\n";

		for(int i=0;i<mr.y.length;i++){

			double dif = mr.y[i]-mr.yhat[i];

			ret += mr.y[i] + " " + mr.yhat[i] + " " + dif + "\n"; 

		}

		ret += "meany\n";

		ret += mr.meany +"\n";

		ret += "vary\n";

		ret += mr.vary + "\n";

		ret += "syy\n";

		ret += mr.syy + "\n";

		ret += "reg\n";

		ret += mr.reg + "\n";

		ret += "resid\n";

		ret += mr.resid + "\n";

		

		ret += "multiple correlation coefficiency\n";

		ret += mr.multcorrcoef + "\n";

		

		ret += "meansquare\n";

		for(int i=0;i<mr.meansquare.length;i++){

			ret += mr.meansquare[i] + " ";

		}

		ret += "\n";

		ret += "F of ANOVA\n";

		ret += mr.anovaF + "\n";

		ret += "df of F\n";

		ret += mr.anovaDF[0] + "," + mr.anovaDF[1] + "\n";

		ret += "x+z\n{";

		for(int i=0;i<mr.xANDz.length;i++){

			ret += mr.xANDz[i] + ",";

		}

		ret += "}\n";

		ret += "x\n{";

		for(int i=0;i<mr.x.length;i++){

			ret += mr.x[i] + ",";

		}

		ret += "}\n";

		ret += "z\n{";

		for(int i=0;i<mr.z.length;i++){

			ret += mr.z[i] + ",";

		}

		ret += "}\n";

		return ret;

		

	}

	
}

• Utils.TestRun

package Utils;

import DiscreteMathTools.*;
public class TestRun {
	/**

	 * @param args

	 */

	public static void main(String[] args) {

		// TODO 自動生成されたメソッド・スタブ

		double[] v= {1,3,2,4.5};

		double mean = Calculator.mean(v);

		System.out.println("mean " + mean);

		double var = Calculator.var(v);

		System.out.println("var " + var);

		double min = Calculator.min(v);

		System.out.println("min " + min);

		double max = Calculator.max(v);

		System.out.println("max " + max);

		double[] sortup = Calculator.sortUp(v);

		System.out.println("sortup");

		for(int i=0;i<sortup.length;i++){

			System.out.println(sortup[i]);

		}

		double[] sortdown = Calculator.sortDown(v);

		System.out.println("sortdown");

		for(int i=0;i<sortdown.length;i++){

			System.out.println(sortdown[i]);

		}

		

		double[] mode = Calculator.mode(v);

		System.out.println("mode");

		for(int i=0;i<mode.length;i++){

			System.out.println(mode[i]);

		}

		double median = Calculator.median(v);

		System.out.println("median " + median);

		

		double[] y = {0.9,1.3,2.0,1.8,2.2,3.5,1.9,2.7,2.1,3.6};

		double[][] mat = {{1.2,1.9},{1.6,2.7},{3.5,3.7},

				{4.0,3.1},{5.6,3.5},{5.7,7.5},{6.7,1.2},{7.5,3.7},{8.5,0.6},{9.7,5.1}};

		

		/*

		 * dummy column を加える

		 */

		double[][] mat2 = MiscUtil.addDummyCol(mat);

		/*

		double[][] mat = {{1,-1,-1,1,-1,1,1,-1},

				{1,1,-1,-1,-1,-1,1,1},

				{1,-1,1,-1,-1,1,-1,1},

				{1,-1,-1,1,1,-1,-1,1},

				{1,1,1,1,-1,-1,-1,-1},

				{1,1,-1,-1,1,1,-1,-1},

				{1,-1,1,-1,1,-1,1,-1},

				{1,1,1,1,1,1,1,1}};

		*/

		/*

		X1 1.2 1.6 3.5 4.0 5.6 5.7 6.7 7.5 8.5 9.7 

		X2 1.9 2.7 3.7 3.1 3.5 7.5 1.2 3.7 0.6 5.1 

		Y 0.9 1.3 2.0 1.8 2.2 3.5 1.9 2.7 2.1 3.6 

		*/

		

		double[][] cov = Calculator.covarmatrix(mat);

		MatrixExec.PrintMat(cov,"\t","\n");

		

		double[][] c = Calculator.cov(y,mat);

		MatrixExec.PrintMat(c,"\t","\n");

		

		double[][] invcov = InvMat.inv(cov);

		MatrixExec.PrintMat(invcov, "\t","\n");

				

		double[][] emat = MatrixExec.Prodmat(cov,invcov);

		MatrixExec.PrintMat(emat,"\t","\n");

		

		double[][] b = MatrixExec.Prodmat(invcov,c);

		MatrixExec.PrintMat(b,"\t","\n");

		

		double meany=Calculator.mean(y);

		double[] meanx = Calculator.meanV(mat);

		double[][] meanXmat = MatrixExec.Vec2Mat(meanx);

		double[][] bt=MatrixExec.Transversmat(b);

		double[][] tmp = MatrixExec.Prodmat(bt,meanXmat);

		double b0 = meany - tmp[0][0];

		System.out.println("b0 " + b0);

		

		System.out.println("

    




    

  

ryamada22 2006-09-17 16:25 読者になる

広告を非表示にする

関連記事

• 2008-09-18

  アニメーションGIF2

  使用したフリーアプリについては、こちらをどうぞ。 ３次元グラ…
• 2008-08-31

  配置する

  だいたいN個の点を半径１の球面上に「ほぼ均一に配置する」 N=1…
• 2008-05-30

  ２次元のヒストグラム

  今、ある値のペアのデータがあるとする。{{x1,y1},{x2,y2},...,…
• 

  2007-07-15

  へのフィット(続きの続きの続き)

  これに先立つ記事は7月10日(こちら)、14日(こちら)。 そこで登…
• 2007-07-08

  関連SNPの統計量空間の狭さ：ラベルパーミュテーション用エクセル２

  昨日、ラベルパーミュテーションをエクセルで簡易に行ってみた…

• もっと読む

コメントを書く