ベータ分布と台形積分による数値解析

はじめに

ベータ分布は、確率分布の一種である。一般に二項分布はコンビネーションを用いて表すが、それだと連続性を議論できない。そこで、連続性を議論することが出来るように、拡張したものがベータ関数である。また、ベータ関数は、定積分によって定義される。そこで今回は、まず台形積分によってベータ関数の値を数値解析する。次に、数値解析した値を基に、ベータ分布を再現することを目的とする。

導入

ベータ分布について

一般的に、二項分布による確率は以下のように定義される。

P(m)=_nC_m p^m(1-p)^{n-m}

ここで、

_nC_m=\frac{n!}{m!(n-m)!} 

これをベータ関数で表したい。
一方で、ベータ関数の定義は以下のように与えられる。

beta(k,l)=\int_0^1x^{k-1}(1-x)^{l-1}dx

これにより、$k,l$が自然数の場合は、以下のことが成立する。

beta(k,l)=\frac{(k-1)!(l-1)!}{(k+l-1)!}=(_{k+l-2}C_{k-1})^{-1}\times\frac{1}{k+l-1}

したがって、

_nC_m=\frac{1}{beta(n-m+1,m+1)(n+1)}

となる。

プログラム

そこで、以下のようなプログラムを書いた。
ただし、$p=n,q=m$とした。

python beta_distribution.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math

S_ary=[]
#区分求積法による数値積分
##積分範囲
a=0
b=1
##分割数
n=10000
# 短冊の幅Δx
dx = (b-a)/n
# 積分する関数の定義
def f(x,p,q):
    return x**(p-1)*(1-x)**(q-1)
def beta_s(p,q):

  s = 0

  for i in range(1,n-1):
    x1 = a + dx*i
    x2 = a + dx*(i+1)
    f1 = f(x1,p,q)    # 上底
    f2 = f(x2,p,q)    # 下底
    # 面積
    s += dx * (f1+f2)/2
  return s
num=100
p1=0.01
sum_p_q=100
p_ary=np.linspace(0,sum_p_q,num)
for i in range(num):
  p=p_ary[i]
  q=sum_p_q-p
  P=(1/(beta_s(p,q)*(sum_p_q+1)))*p1**(p-1)*(1-p1)**(q-1)
  S_ary.append(P)
plt.plot(p_ary,S_ary)
title="ベータ分布p="+str(p1)
plt.title(title)
plt.xlabel("m回事象が起こる確率")
plt.ylabel("確率")
plt.savefig(title+".png")
plt.show()

結果

このプログラムを実行すると以下のようなグラフが出力される。

このことから、二項分布や正規分布、ポアソン分布と分布が似ていることが分かる。

まとめ

今回は、確率分布の一種であるベータ分布を二項分布の連続型への拡張とみて再現した。その際に、ベータ関数の値を求めるために台形積分を用いた。結果、二項分布や正規分布、成功する試行がまれな場合はポアソン分布に酷似するということが分かった。このように、確率統計の分野は線形代数や微分積分の計算と繋がっているので、それらを合わせて考えると面白い。

参考文献

ベータ分布について

台形積分の数値解析について