正規分布の計算
累積確率密度関数
と確率密度関数
が次式で与えられる標準正規分布の計算プログラムをPythonで作成した。
まず、確率変数の値がz以下である確率を求めるプログラムは以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝ CumNormalDistri.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
from adap_gl import
AdapGL as agl
import math
class NormalDistri:
def __init__(self):
self.c = 1.0 / math.sqrt(2.0 * math.pi)
def v(self, z):
return self.c * math.exp(-0.5 * z * z)
n_func =
NormalDistri()
gl = agl()
def cum_normal( z ):
if z > 0.0:
return 0.5 + gl.integral(0.0, z, n_func)
elif z < 0.0:
return 0.5 - gl.integral(z, 0.0, n_func)
else:
return 0.5
while True:
s = input("z = ")
try:
z =
float(s)
p = cum_normal(z)
print("P(x < %f) = %f"%(z, p))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
確率密度関数は、積分を行うクラス型AdapGL (agl) のオブジェクトglの引数とするため、クラス型NormalDistriのオブジェクトn_funcとして用意している。積分を行うクラス型AdapGLについては、＜このウェブサイト＞を参照されたい。
プログラムCumNormalDistri.pyの実行例を図１に示す。
図１
プログラムを実行すると
「z = 」
と、zの値を聞いてくる。計算したい値を設定して、Enterキーを押すと、計算が始まり、
「P(x < …」
の形式で結果が表示される。
ｚの入力部は、以下のコードで行われている。
s = input("z = ")
try:
z = float(s)
p = cum_normal(z)
print("P(x < %f) = %f"%(z, p))
except:
break
実数値に変換できない文字列が入力されると、例外が生成されて、breakが実行され、プログラムの実行終了となる。
上のプログラムとは逆に、累積確率pに対応するｚの値を求めるプログラムが次のCalc_z_NormalDistri.pyである。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Calc_z_NormalDistri.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
from adap_gl import
AdapGL as agl
import math
from calcroot import
Newton
class NormalDistri:
def __init__(self):
self.c = 1.0 / math.sqrt(2.0 * math.pi)
def v(self, z):
return self.c *
math.exp(-0.5 * z * z)
n_func =
NormalDistri()
gl = agl()
def cum_normal( z ):
if z > 0.0:
return 0.5 + gl.integral(0.0, z, n_func)
elif z < 0.0:
return 0.5 - gl.integral(z, 0.0, n_func)
else:
return
0.5
class Ccum_normal:
def __init__(self):
pass
def v(self, z):
return cum_normal(z)
f_cum_n =
Ccum_normal()
class DFunc_normal:
def __init__(self):
pass
def v(self, z):
return n_func.v(z)
df_cum_n =
DFunc_normal()
def calc_z( p ):
Lz = 0.0
while (f_cum_n.v(Lz) >=
p):
Lz
-= 1.0
Uz = Lz + 1.0
while (f_cum_n.v(Uz) <=
p):
Uz
+= 1.0
z = Newton( f_cum_n,
df_cum_n, p, Lz, Uz )
return z
while True:
try:
s =
input("p = ")
p =
float(s)
if
((p <= 0.0) or (p >= 1.0)):
break
z =
calc_z(p)
print("P(z <= %f) = %f"%(z, p))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
累積確率を求める関数cum_normal( z )の方程式に対して根を求める関数Newtonを適用して、累積確率ｐに対するｚを求める関数calc_z(p)が構成されている。方程式の根を求める関数Newtonについては＜このウェブサイト＞を参照されたい。
累積確率に対するｚの値を求めるプログラムCalc_z_NormalDistri.pyの実行例を図２に示す。
図２
累積確率ｐの設定を
「p = 」
と聞いてくるので、計算したいｐの値を設定する。Enterキーを押すと、計算が始まり
「P(z <= …」
の形式で計算結果が表示される。
プログラムの実行を終了したいときは、数値に変換できない文字列か確率の範囲外の値を設定すればよい。図２では、文字qが入力されている。
プログラムを圧縮ファイルnormal_distri_prgs.zipにまとめた。ファイル名normal_distri_prgs.zipをクリックしてダウンロードすることができる。