分布の計算
カイ2乗分布の計算プログラムをPythonで作ってみた。
カイ2乗分布の確率密度関数は次式で与えられる。
カイ2乗の値
に対する累積確率
を求めるプログラムを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Calc_p_from_chi2.py＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
from adap_gl import
AdapGL as agl
#from calcroot import
Newton
#from
Calc_z_NormalDistri import calc_z
def sqr( x ):
return x * x
class
C_chi2_tan_kernel:
def __init__(self, df):
self.df = df
def v(self, theta):
x =
math.tan(theta)
v =
pow(x, (self.df/2.0) - 1.0) * math.exp(-0.5 * x)
return v / sqr(math.cos(theta))
class
C_chi2_tan_density:
def __init__(self, df):
self.df = df
chi2_kernel = C_chi2_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( 0.0, math.pi/2.0, chi2_kernel)
self.c = 1.0 / v
def v(self, theta):
x =
math.tan(theta)
v =
pow(x, (self.df/2.0) - 1.0) * math.exp(-0.5 * x)
return self.c * v / sqr(math.cos(theta))
class C_cum_chi2:
def __init__(self, df):
self.df = df
self.chi2_tan = C_chi2_tan_density(self.df)
self.gl = agl()
def v( self, x ):
chi2_cum = self.gl.integral(0.0, math.atan(x), self.chi2_tan)
return chi2_cum
while True:
try:
s =
input("df = ")
df =
float(s)
cum_chi2 = C_cum_chi2(df)
s =
input("chi^2 = ")
chi2 = float(s)
p =
cum_chi2.v(chi2)
print("p = %f
chi^2 = %f"%(p, chi2))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
実行例を図１に示す。
図１
自由度dfとカイ２乗chi^2の値を設定すると累積確率ｐが算出される。自由度dfに文字を設定するとプログラムの実行終了となる。図１では、文字ｐを設定して終了している。
累積確率ｐに対するカイ2乗の値を求めるプログラムを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Calc_chi2_from_p.py＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
from adap_gl import
AdapGL as agl
from calcroot import
Newton
#from
Calc_z_NormalDistri import calc_z
def sqr( x ):
return x * x
class
C_chi2_tan_kernel:
def __init__(self, df):
self.df = df
def v(self, theta):
x =
math.tan(theta)
v =
pow(x, (self.df/2.0) - 1.0) * math.exp(-0.5 * x)
return v / sqr(math.cos(theta))
class
C_chi2_tan_density:
def __init__(self, df):
self.df = df
chi2_kernel = C_chi2_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( 0.0, math.pi/2.0, chi2_kernel)
self.c = 1.0 / v
def v(self, theta):
x =
math.tan(theta)
v =
pow(x, (self.df/2.0) - 1.0) * math.exp(-0.5 * x)
return self.c * v / sqr(math.cos(theta))
class C_cum_chi2:
def __init__(self, df):
self.df = df
self.chi2_tan = C_chi2_tan_density(self.df)
self.gl = agl()
def v( self, x ):
chi2_cum = self.gl.integral(0.0, math.atan(x), self.chi2_tan)
return chi2_cum
class C_dcum_chi2:
def __init__(self, df):
self.df = df
chi2_kernel = C_chi2_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( 0.0, math.pi/2.0, chi2_kernel)
self.c = 1.0 / v
def v(self, x):
v =
self.c * pow(x, (self.df/2.0) - 1.0) * math.exp(-0.5 * x)
return v
def calc_chi2_from_p(
p, df ):
cum_chi2 = C_cum_chi2( df )
Lb = df
while (cum_chi2.v(Lb) >=
p):
Lb
*= 0.5
Ub = Lb + 1.0
while (cum_chi2.v(Ub) <=
p):
Ub
+= 1.0
dcum_chi2 = C_dcum_chi2( df
)
chi2 = Newton( cum_chi2,
dcum_chi2, p, Lb, Ub )
return chi2
while True:
try:
s =
input("df = ")
df =
float(s)
s =
input("p = ")
p =
float(s)
chi2
= calc_chi2_from_p( p, df )
print("chi^2 = %f p = %f df = %f"%(chi2, p, df))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
実行例を図２に示す。
図２
自由度dfと累積確率ｐを設定すると、対応するカイ2乗の値chi^2が計算される。自由度に文字を設定するとプログラムの実行終了となる。図２では、文字aを設定して終了している。
プログラムは圧縮ファイルchi2_distriPrgs.zipにまとめた。ファイル名chi2_distriPrgs.zipをクリックしてダウンロードすることができる。