Ｆ分布の計算
Ｆ分布に従う確率変数Ｘに対して、Ｆの値に対する累積確率
、および累積確率ｐに対するＦの値を求める。Ｆ分布は自由度と呼ばれている２つのパラメータdf1とdf2をもち、次式（１）で与えられる。
プログラムでは、Ｆ分布の自由度と確率密度関数（１）の係数部分をクラス型のアトリビュートに設定して、（１）式を１変数の関数としてガウスの積分公式を適応的に適用している。具体的には以下のプログラムとした。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝ calc_p_for_F.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
from adap_gl import
AdapGL as agl
from calcroot import
Newton
class C_F_tan_kernel:
def __init__(self, df1,
df2):
self.df1 = df1
self.df2 = df2
def v(self, theta):
x =
math.tan(theta)
v =
pow(x, (self.df1 / 2.0) - 1.0) / \
pow(1.0 + (self.df1 * x / self.df2), (self.df1 + self.df2) / 2.0)
return v / (math.cos(theta) ** 2.0)
class C_F_tan_density:
def __init__(self, df1, df2):
self.df1 = df1
self.df2 = df2
F_kernel =C_F_tan_kernel(self.df1, self.df2)
gl =
agl()
v =
gl.integral(0.0, math.pi/2.0, F_kernel)
self.c = 1.0 / v
def v(self, theta):
x =
math.tan(theta)
v = pow(x, (self.df1 / 2.0) -
1.0) / \
pow(1.0 + (self.df1 * x / self.df2), (self.df1 + self.df2) / 2.0)
return self.c * v / (math.cos(theta) ** 2.0)
class C_cum_F:
def __init__(self, df1,
df2):
self.df1 = df1
self.df2 = df2
self.F_tan = C_F_tan_density(self.df1, self.df2)
self.gl = agl()
def v(self, F):
F_cum = self.gl.integral(0.0, math.atan(F), self.F_tan)
return F_cum
while True:
try:
s =
input("df1 = ")
df1
= float(s)
s =
input("df2 = ")
df2
= float(s)
s =
input("F = ")
F =
float(s)
cum_F = C_cum_F(df1, df2)
p =
cum_F.v(F)
print("p = %f F =
%f"%(p, F))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
無限区間の場合にもガウスの積分公式を適用するため、変数変換
を行っている。
上のプログラムcalc_p_for_F.pyの実行例を図１に示す。
図１
２つの自由度とＦの値を設定すると、対応する累積確率が計算されて表示される。自由度などの設定において、数値以外の文字を入力するとプログラムの実行終了となる。図１では文字 q が入力されて、プログラムの実行終了となっている。
累積確率
に対して対応するＦの値を求めるプログラムcalc_F_for_p.pyを以下のように作成した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝ calc_F_for_p.py＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import math
from adap_gl import
AdapGL as agl
from calcroot import
Newton
class C_F_tan_kernel:
def __init__(self, df1,
df2):
self.df1 = df1
self.df2 = df2
def v(self, theta):
x =
math.tan(theta)
v =
pow(x, (self.df1 / 2.0) - 1.0) / \
pow(1.0 + (self.df1 * x / self.df2), (self.df1 + self.df2) / 2.0)
return v / (math.cos(theta) ** 2.0)
class C_F_tan_density:
def __init__(self, df1,
df2):
self.df1 = df1
self.df2 =
df2
F_kernel =C_F_tan_kernel(self.df1, self.df2)
gl =
agl()
v =
gl.integral(0.0, math.pi/2.0, F_kernel)
self.c = 1.0 / v
def v(self, theta):
x =
math.tan(theta)
v =
pow(x, (self.df1 / 2.0) - 1.0) / \
pow(1.0 + (self.df1 * x / self.df2), (self.df1 + self.df2) / 2.0)
return self.c * v / (math.cos(theta) ** 2.0)
class C_cum_F:
def __init__(self, df1,
df2):
self.df1 = df1
self.df2 = df2
self.F_tan = C_F_tan_density(self.df1, self.df2)
self.gl = agl()
def v(self, F):
F_cum = self.gl.integral(0.0, math.atan(F), self.F_tan)
return F_cum
class C_dcum_F:
def __init__(self, df1,
df2):
self.df1 = df1
self.df2 = df2
F_kernel =C_F_tan_kernel(self.df1, self.df2)
gl =
agl()
v =
gl.integral(0.0, math.pi/2.0, F_kernel)
self.c = 1.0 / v
def v(self, F):
v =
pow(F, (self.df1 / 2.0) - 1.0) / \
pow(1.0 +
(self.df1 * F / self.df2), (self.df1 + self.df2) / 2.0)
return self.c * v
def calc_F_for_p( p,
df1, df2 ):
cum_F = C_cum_F(df1, df2)
Lb = 1.0
while (cum_F.v(Lb) >= p):
Lb
*= 0.5
Ub = Lb * 2.0
while (cum_F.v(Ub) <= p):
Ub
+= 1.0
dcum_F = C_dcum_F(df1, df2)
F = Newton( cum_F, dcum_F,
p, Lb, Ub )
return F
while True:
try:
s =
input("df1 = ")
df1
= float(s)
s =
input("df2 = ")
df2
= float(s)
s =
input("p = ")
p =
float(s)
F =
calc_F_for_p( p, df1, df2 )
print("F = %f p =
%f df1 = %f df2 = %f"%(F, p, df1, df2))
except:
break
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
プログラムの実行例を図２に示す。
図２
自由度df1とdf2、および累積確率pを設定すると対応するＦの値が算出される。自由度に数値に変換されない文字を入力するとプログラムの実行終了となる。図２では文字 a が入力されてプログラムの実行終了となっている。
プログラムファイルは圧縮ファイルF_DistriPrgs.zipにまとめた。ファイル名F_DistriPrgs.zipをクリックしてダウンロードできる。