Up

正規分布の計算

 

累積確率密度関数と確率密度関数が次式で与えられる標準正規分布の計算プログラムを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の実行例を図1に示す。

図1

 

プログラムを実行すると

z =

と、zの値を聞いてくる。計算したい値を設定して、Enterキーを押すと、計算が始まり、

P(x < …」

の形式で結果が表示される。

zの入力部は、以下のコードで行われている。

    s = input("z = ")

    try:

        z = float(s)

        p = cum_normal(z)

        print("P(x < %f) = %f"%(z, p))

    except:

        break

実数値に変換できない文字列が入力されると、例外が生成されて、breakが実行され、プログラムの実行終了となる。

 

 

上のプログラムとは逆に、累積確率pに対応するzの値を求めるプログラムが次の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を適用して、累積確率pに対するzを求める関数calc_z(p)が構成されている。方程式の根を求める関数Newtonについては<このウェブサイト>を参照されたい。

累積確率に対するzの値を求めるプログラムCalc_z_NormalDistri.pyの実行例を図2に示す。

図2

 

累積確率pの設定を

p =

と聞いてくるので、計算したいpの値を設定する。Enterキーを押すと、計算が始まり

P(z <= …」

の形式で計算結果が表示される。

プログラムの実行を終了したいときは、数値に変換できない文字列か確率の範囲外の値を設定すればよい。図2では、文字qが入力されている。

 

プログラムを圧縮ファイルnormal_distri_prgs.zipにまとめた。ファイル名normal_distri_prgs.zipをクリックしてダウンロードすることができる。

 

 

Up