ｔ分布の計算
確率分布としてｔ分布と呼ばれているものについて、t 以下の値をとる累積確率ｐを計算するプログラムと、累積確率 p に対応する t を求めるプログラムを作成した。ｔ分布の確率密度関数は次式で与えられる。
ここで、d は自由度と呼ばれているパラメータである。パラメータ d と、分母の値をクラス型の中にアトリビュート（attribute）として設定して、クラス型内の１変数関数を表すメソッドを積分するクラス型により累積確率
を算出する。
次のプログラムcalc_p_from_t.pyでは、クラス型C_cum_tのオブジェクトcum_tによって（１）式のｔに対する累積確率が計算されている。cum_tは、自由度df によって
cum_t = C_cum_t(df)
と生成されている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ calc_p_from_t.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
from adap_gl import
AdapGL as agl
def sqr( x ):
return x * x
class C_t_tan_kernel:
def __init__(self, df):
self.df = df
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return v / sqr(math.cos(theta))
class C_t_tan_density:
def __init__(self, df):
self.df = df
f_kernel = C_t_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( -math.pi/2.0, math.pi/2.0, f_kernel)
self.c = 1.0 / v
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return self.c * v / sqr(math.cos(theta))
class C_cum_t:
def __init__(self, df):
self.df = df
self.t_tan = C_t_tan_density(self.df)
self.gl = agl()
def v( self, t ):
t_cum = self.gl.integral(-math.pi/2.0, math.atan(t), self.t_tan)
return t_cum
while True:
try:
s =
input("df = ")
df =
float(s)
cum_t = C_cum_t(df)
s =
input("t = ")
t =
float(s)
p =
cum_t.v(t)
print("p = %f t =
%f"%(p, t))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のプログラムcalc_p_from_t.pyの実行例を図１に示す。
図１
数値の入力時に英字など数値に変換できない入力を行うとプログラムの実行終了となる。図１では、最後に英字 q が入力されている。
次のプログラムCalc_t_from_p.pyでは、累積確率ｐに対応するｔの値が関数calc_p_to_tによって求められている。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Calc_t_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_t_tan_kernel:
def __init__(self, df):
self.df = df
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return v / sqr(math.cos(theta))
class C_t_tan_density:
def __init__(self, df):
self.df = df
f_kernel = C_t_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( -math.pi/2.0, math.pi/2.0, f_kernel)
self.c = 1.0 / v
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return self.c * v / sqr(math.cos(theta))
class C_cum_t:
def __init__(self, df):
self.df = df
self.t_tan = C_t_tan_density(self.df)
self.gl = agl()
def v( self, t ):
t_cum = self.gl.integral(-math.pi/2.0, math.atan(t), self.t_tan)
return t_cum
class C_dcum_t:
def __init__(self, df):
self.df = df
f_kernel = C_t_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( -math.pi/2.0, math.pi/2.0, f_kernel)
self.c = 1.0 / v
def v(self, t):
tv =
self.c * pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return tv
def calc_p_to_t( p, df
):
cum_t = C_cum_t( df )
Lb = 0.0
while (cum_t.v(Lb) >= p):
Lb
-= 1.0
Ub = Lb + 1.0
while (cum_t.v(Ub) <= p):
Ub
+= 1.0
dcum_t = C_dcum_t( df )
t = Newton( cum_t, dcum_t,
p, Lb, Ub )
return t
while True:
try:
s =
input("df = ")
df =
float(s)
s =
input("p = ")
p =
float(s)
t =
calc_p_to_t( p, df )
print("t = %f p = %f df = %f"%(t, p, df))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
実行例を図２に示す。
図２
図２の場合も、数値に変換されない文字ｑの入力でプログラムの実行終了となっている。
ｔ分布
は、自由度 d が大きくなると正規分布に近づく。自由度d=1000とおいて正規分布と比較してみた。
まず、累積確率をｔ分布と標準正規分布とで比較するプログラムを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝ Check_p_from_t_ByNormalD.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
from adap_gl import
AdapGL as agl
def sqr( x ):
return x * x
class C_t_tan_kernel:
def __init__(self, df):
self.df = df
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return v / sqr(math.cos(theta))
class C_t_tan_density:
def __init__(self, df):
self.df = df
f_kernel =
C_t_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( -math.pi/2.0, math.pi/2.0, f_kernel)
print("v = %f"%v)
self.c = 1.0 / v
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return self.c * v / sqr(math.cos(theta))
class C_cum_t:
def __init__(self, df):
self.df = df
self.t_tan = C_t_tan_density(self.df)
self.gl = agl()
def v( self, t ):
t_cum = self.gl.integral(-math.pi/2.0, math.atan(t), self.t_tan)
return t_cum
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
cum_t = C_cum_t(1000)
t = -10.0
while True:
p = cum_t.v(t)
pn = cum_normal(t)
print("t = %f p = %f pn = %f p - pn = %f"%(t,
p, pn, p - pn))
t += 0.5
if t > 10.1:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
実行結果を図３に示す。
図３
次に、逆に累積確率ｐから対応するｔ値を求めて、標準正規分布の場合と比較するプログラムを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ Check_t_from_p_ByNormalD.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_t_tan_kernel:
def __init__(self, df):
self.df = df
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return v / sqr(math.cos(theta))
class C_t_tan_density:
def __init__(self, df):
self.df = df
f_kernel = C_t_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( -math.pi/2.0, math.pi/2.0, f_kernel)
self.c = 1.0 / v
def v(self, theta):
t =
math.tan(theta)
v =
pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return self.c * v / sqr(math.cos(theta))
class C_cum_t:
def __init__(self, df):
self.df = df
self.t_tan = C_t_tan_density(self.df)
self.gl = agl()
def v( self, t ):
t_cum = self.gl.integral(-math.pi/2.0, math.atan(t), self.t_tan)
return t_cum
class C_dcum_t:
def __init__(self, df):
self.df = df
f_kernel = C_t_tan_kernel(self.df)
gl =
agl()
v =
gl.integral( -math.pi/2.0, math.pi/2.0, f_kernel)
self.c = 1.0 / v
def v(self, t):
tv =
self.c * pow(1.0 / (1.0 + sqr(t) / self.df), (self.df + 1.0) / 2.0)
return tv
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
def calc_p_to_t( p, df
):
cum_t = C_cum_t( df )
Lb = 0.0
while (cum_t.v(Lb) >= p):
Lb
-= 1.0
Ub = Lb + 1.0
while (cum_t.v(Ub) <= p):
Ub
+= 1.0
dcum_t = C_dcum_t( df )
t = Newton( cum_t, dcum_t,
p, Lb, Ub )
return t
p = 0.05
while (True):
t = calc_p_to_t( p, 1000 )
z = calc_z( p )
print("p = %f t = %f z = %f t - z = %f"%(p,
t, z, t - z))
p += 0.05
if (p > 0.97):
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
実行結果は図４のようになる。
図４
プログラムは圧縮ファイルt_distri_prgs.zipにまとめた。ファイル名t_distri_prgs.zipをクリックしてダウンロードできる。