確率計算

  • ディリクレ分布については、Wikipediaの記事がこちら
    • このWikipedia記事のは、3つの独立事象について、Dir(a=\{a1,a2,a3\})a1=a2=a3の場合について、0.3<=a1<=2の範囲で0<=x1,x2,x3<=1,x1=x2=x3の範囲について計算し、その結果を、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)
}
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;
	}