確率計算

• ディリクレ分布については、Wikipediaの記事がこちら。
  • このWikipedia記事の図は、３つの独立事象について、 をの場合について、の範囲での範囲について計算し、その結果を、x1,x2の軸についてプロットしたものである。
  • これをエクセルで計算するとこんな感じ
  • Rで、関数にするとこんな感じ(確率の対数を返す。xには、総和が1となるようなベクトルを与えること)

pdirichletLN<-function(x,a){
b=0
for(i in 1:length(a)) b=b+log(gamma(a[i]))
c=0
for(i in 1:length(a)) c=c+a[i]
d=log(gamma(c))
g=0
for(i in 1:length(a)) g=g+log(x[i])*(a[i]-1)

f=g-b+d
return(f)
}

• • Javaだと

public static double pDirichletLN(double[] x,double[] a){
		double ret=0;
		double b=0;
		for(int i=0;i<a.length;i++){
			b+=org.apache.commons.math.special.Gamma.logGamma(a[i]);
		}
		double c=0;
		for(int i=0;i<a.length;i++){
			c+=a[i];
		}
		double d=org.apache.commons.math.special.Gamma.logGamma(c);
		double g=0;
		for(int i=0;i<a.length;i++){
			g+=Math.log(x[i])*(a[i]-1);
		}
		ret=g-b+d;
		return ret;
	}