元理系院生の新入社員がPythonとJavaで色々頑張るブログ

プログラミングや機械学習について調べた事を書いていきます

【Java】javaで多次元正規分布

旋回作成した行列計算用のクラスを用いて多次元正規分布を求めるクラスを作成します。emoson.hateblo.jp

多次元正規分布については以下の記事で触れました。emoson.hateblo.jp


コードは次のようになります

import emoson.math.Matrix2D;

public class NormalDistribution {
	private Matrix2D VCM;		//分散共分散行列
	private Matrix2D invVCM;	//分散共分散行列の逆行列
	private Matrix2D Mean;		//平均ベクトル
	
	public NormalDistribution(Matrix2D vcm, Matrix2D me){
		this.VCM = vcm;
		this.invVCM = Matrix2D.inv(vcm);
		this.Mean = me;
	}
	
	public NormalDistribution(double[][] vcm, double[] me){
		this.VCM = new Matrix2D(vcm);
		this.invVCM = Matrix2D.inv(VCM);
		this.Mean = new Matrix2D(me);
	}
	
	//入力に対する確率密度関数の値を算出
	public double getProb(Matrix2D x){
		double a = Math.exp(Matrix2D.dot(Matrix2D.dot(Matrix2D.prod(-0.5, Matrix2D.sub(x, Mean)).T(), invVCM), Matrix2D.sub(x, Mean)).getArrays()[0][0]);
		double b = Matrix2D.prod(Math.pow(Math.PI*2, (double)VCM.getRow()/2), Matrix2D.sqrt(VCM)).getArrays()[0][0];
		return a / b;
	}
	
	//配列でもOK
	public double getProb(double[] x){
		return getProb(new Matrix2D(x));
	}

}

このクラスはコンストラクタで分散共分散行列と平均ベクトルを与えることで、任意の確率密度関数を作成する事が出来ます。

ある値に対する確率密度関数の値を求める際にはgetProb関数を使用します。

分散共分散行列 \Sigma=\left|
    \begin{array}{ccc}
      a_{1.5} & a_{0.3} \\
      a_{0.3} & a_{2.0} \\
    \end{array}
  \right|、平均ベクトルを \mu=\left(
    \begin{array}{c}
      0.0 \\
      0.0 \\
    \end{array}
  \right)
とした時の正規分布を求めるプログラムは次のようになります

Matrix2D vcm = new Matrix2D(new double[][]{{1.5, 0.3}, {0.3, 2.0}});
Matrix2D mu = new Matrix2D(new double[]{0.0, 0.0});
NormalDistribution g = new NormalDistribution(vcm, mu);
double y=-1.0;
while(y<=1.0){
	double x = -1.0;
	while(x <= 1.0){
	        System.out.print(g.getProb(new double[]{x, y})+"\t");
		x+=0.1;
	}
	System.out.println();
	y+=0.1;
}

このプログラムでは2つの変数を-1から+1まで0.1刻みで変化させ、その際の確率密度関数の値を求めています。
この結果をエクセルにペタッと貼って、値の大小に応じてグラデーションをかけると次のようになります。
f:id:emoson:20151030162440p:plain

いい感じに正規分布になっていますね。