ベータ分布 - gerumaniumの日記

前回ガンマ関数を書いたので、これを使ってベータ分布を書きます。

ベータ分布は定義としては、以下のように与えられます。
$Beta(\mu|a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\mu^{a-1}(1-\mu)^{b-1}$

#-*- coding:utf-8 -*-
#ベータ分布

from numpy import *
from pylab import *

PI=3.14159265358979324
LOG_2PI=1.837877066400934548
N=8
B0=1
B1=(-1.0/2.0)
B2=(1.0/6.0)
B4=(1.0/30.0)
B6=(1.0/42.0)
B8=(-1.0/30.0)
B10=(5.0/66.0)
B12=(-691.0/2730.0)
B14=(7.0/6.0)
B16=(-3617.0/510.0)
B18=(43867.0/798.0)
B20=(-174611.0/330.0)

def beta(x,a,b):
	return (gamma(a+b)/( gamma(a)*gamma(b) )) * (x**(a-1.0)) * ((1.0-x)**(b-1.0))

def log_gamma(x):
	v=1
	while x < N:
		v*=x
		x+=1
	w=1/(x*x)
	return 0.5*LOG_2PI - log(v) -x + (x-0.5) *log(x) +((((((((  (B16/(16*15))*w + (B14/(14*13))) *w + (B12/(12*11))) *w + (B10*9)))*w +(B8/(8*7)))*w + (B6/(6*5)))*w +(B4/(4*3))) *w + (B2/(2*1)))/x


def gamma(x):
	if x < 0:
		return PI/(sin(PI*x)*exp(log_gamma(1-x)))
	return exp(log_gamma(x))

if __name__=="__main__":
	print "超パラメータa,bを決めます。"
	print "aを入力してください"
	a=input()
	print "bを入力してください"
	b=input()

	xlist=linspace(0,1,10000)
	ylist=[]
	for i in xlist:
		ylist.append(beta(i,a,b))
	
	plot(xlist,ylist,'r-')
	ylim(0,3)
	show()

PRMLの図表に従って、出力した。

以下引用です。
"一分でつくれるベータ分布"

#ガンマ関数で書く。
import matplotlib
import pylab
import scipy

x=[0.01*i for i in range(100)]
scipy.pkglod("special")
gamma=scipy.special.gamma
def dbeta(x,a,b):
 return gamma(a+b)/(gamma(a)*gamma(b))*x**(a-1)*(1-x)**(b-1)
matplotlib.pyplot.plot([dbeta(xi,3,2) for xi in x])

#ベータ確率密度関数を用いる。
import matplotlib
import pylab
import scipy

x=[0.01*i for i in range(100)]
scipy.pkgload("stats")
a=[scipy.stats.beta.pdf(xi,3,2) for xi in x]
matplotlib.pyplot.plot(a)
matplotlib.pyplot.show()