ベータ関数をプロットしてみる
a,bを自然数に絞った手抜き版。
> import matplotlib > import pylab > x = [ 0.01*i for i in range(100)] > def minigamma(a): b = a-1 if b > 1: return b * minigamma(b) else: return 1 > def dbeta(x, a, b): return minigamma(a+b)/(minigamma(a)*minigamma(b)) * x **(a-1) * (1-x)**(b-1) >matplotlib.pyplot.plot([ dbeta(xi , 3, 2) for xi in x])
こんなふうにできた。
ちなみにRだと
x <- seq(0,1, length = 100)
plot(x, dbeta(x, 3, 2))
だけでいけます。pythonでは、Rのrbeta相当はrandom.betavariateでイケそうですが、dbetaがなさそうです。
Dirichlet
def ddirchlet(x, a): b = map(lambda n: x**(n-1), a) a0 = sum(a) return minigamma(a0)/(reduce(lambda x,y: x *minigamma(y),a)) * reduce( lambda x,y: x * y, b)
追記
ガンマ関数間違ってるので直した。グラフの形も若干変わるかも。
