############################
# ｔ分布の計算
# coded by Y.Suganuma
############################
############################################
# Γ（ｘ）の計算（ガンマ関数，近似式）
# ier : =0 : normal
# =-1 : x=-n (n=0,1,2,･･･)
# return : 結果
# coded by Y.Suganuma
############################################
def gamma(x, ier)
ier[0] = 0
if x > 5.0
v = 1.0 / x
s = ((((((-0.000592166437354 * v + 0.0000697281375837) * v + 0.00078403922172) * v - 0.000229472093621) * v - 0.00268132716049) * v + 0.00347222222222) * v + 0.0833333333333) * v + 1.0
g = 2.506628274631001 * Math.exp(-x) * (x ** (x-0.5)) * s
else
err = 1.0e-20
w = x
t = 1.0
if x < 1.5
if x < err
k = Integer(x)
y = Float(k) - x
if y.abs() < err or (1.0-y).abs() < err
ier[0] = -1
end
end
if ier[0] == 0
while w < 1.5
t /= w
w += 1.0
end
end
else
if w > 2.5
while w > 2.5
w -= 1.0
t *= w
end
end
end
w -= 2.0
g = (((((((0.0021385778 * w - 0.0034961289) * w + 0.0122995771) * w - 0.00012513767) * w + 0.0740648982) * w + 0.0815652323) * w + 0.411849671) * w + 0.422784604) * w + 0.999999926
g *= t
end
return g
end
################################################
# 標準正規分布N(0,1)の計算（P(X = x), P(X < x)）
# w : P(X = x)
# return : P(X < x)
################################################
def normal(x, w)
# 確率密度関数（定義式）
w[0] = Math.exp(-0.5 * x * x) / Math.sqrt(2.0*Math::PI)
# 確率分布関数（近似式を使用）
y = 0.70710678118654 * x.abs()
z = 1.0 + y * (0.0705230784 + y * (0.0422820123 + y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 + y * 0.0000430638)))))
pp = 1.0 - z ** (-16.0)
if x < 0.0
pp = 0.5 - 0.5 * pp
else
pp = 0.5 + 0.5 * pp
end
return pp
end
################################################
# ｔ分布の計算（P(X = tt), P(X < tt)）
# （自由度が∞の時の値はN(0,1)を利用して下さい）
# dd : P(X = tt)
# df : 自由度
# return : P(X < tt)
################################################
def t(tt, dd, df)
if tt < 0.0
sign = -1.0
else
sign = 1.0
end
if tt.abs() < 1.0e-10
tt = sign * 1.0e-10
end
t2 = tt * tt
x = t2 / (t2 + df)
if df%2 != 0
u = Math.sqrt(x*(1.0-x)) / Math::PI
p = 1.0 - 2.0 * Math.atan2(Math.sqrt(1.0-x), Math.sqrt(x)) / Math::PI
ia = 1
else
u = Math.sqrt(x) * (1.0 - x) / 2.0
p = Math.sqrt(x)
ia = 2
end
if ia != df
i1 = ia
while i1 < df-1
p += 2.0 * u / i1
u *= ((1.0 + i1) / i1 * (1.0 - x))
i1 += 2
end
end
dd[0] = u / tt.abs()
pp = 0.5 + 0.5 * sign * p
return pp
end
############################################
# 二分法による非線形方程式(f(x)=0)の解
# x1,x2 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜f(x(k))｜＜eps2）
# max : 最大試行回数
# ind : 実際の試行回数
# （負の時は解を得ることができなかった）
# fn : f(x)を計算する関数名
# return : 解
# coded by Y.Suganuma
############################################
def bisection(x1, x2, eps1, eps2, max, ind, &fn)
x0 = 0.0
f1 = fn.call(x1)
f2 = fn.call(x2)
if f1*f2 > 0.0
ind[0] = -1
else
ind[0] = 0
if f1*f2 == 0.0
if f1 == 0.0
x0 = x1
else
x0 = x2
end
else
sw = 0
while sw == 0 && ind[0] >= 0
sw = 1
ind[0] += 1
x0 = 0.5 * (x1 + x2)
f0 = fn.call(x0)
if f0.abs() > eps2
if ind[0] <= max
if (x1-x2).abs() > eps1 && (x1-x2).abs() > eps1*x2.abs()
sw = 0
if f0*f1 < 0.0
x2 = x0
f2 = f0
else
x1 = x0
f1 = f0
end
end
else
ind[0] = -1
end
end
end
end
end
return x0
end
############################################
# Newton法による非線形方程式(f(x)=0)の解
# x0 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜f(x(k))｜＜eps2）
# max : 最大試行回数
# ind : 実際の試行回数
# （負の時は解を得ることができなかった）
# fn : f(x)とその微分を計算する関数名
# return : 解
# coded by Y.Suganuma
############################################
def newton(x0, eps1, eps2, max, ind, &fn)
x1 = x0
x = x1
ind[0] = 0
sw = 0
while sw == 0 and ind[0] >= 0
sw = 1
ind[0] += 1
g = fn.call(0, x1)
if g.abs() > eps2
if ind[0] <= max
dg = fn.call(1, x1)
if dg.abs() > eps2
x = x1 - g / dg
if (x-x1).abs() > eps1 && (x-x1).abs() > eps1*x.abs()
x1 = x
sw = 0
end
else
ind[0] = -1
end
else
ind[0] = -1
end
end
end
return x
end
################################################################
# 標準正規分布N(0,1)のｐ％値（P(X > u) = 0.01p）（二分法を使用）
# ind : >= 0 : normal（収束回数）
# = -1 : 収束しなかった
################################################################
def p_normal(ind)
normal_snx = Proc.new { |x|
y = Array.new(1)
1.0 - $p - normal(x, y)
}
u = bisection(-7.0, 7.0, 1.0e-6, 1.0e-10, 100, ind, &normal_snx)
return u
end
#################################################
# ｔ分布のｐ％値（P(X > u) = 0.01p）
# （自由度が∞の時の値はN(0,1)を利用して下さい）
# ind : >= 0 : normal（収束回数）
# = -1 : 収束しなかった
#################################################
def p_t(ind)
t_snx = Proc.new { |sw, x|
y = Array.new(1)
z = t(x, y, $dof)
if sw == 0
z - 1.0 + $p
else
z = y[0]
end
}
tt = 0.0
pis = Math.sqrt(Math::PI)
df = Float($dof)
df2 = 0.5 * df
# 自由度＝１
if $dof == 1
tt = Math.tan(Math::PI*(0.5-$p))
else
# 自由度＝２
if $dof == 2
if $p > 0.5
c = -1.0
else
c = 1.0
end
p2 = (1.0 - 2.0 * $p)
p2 *= p2
tt = c * Math.sqrt(2.0 * p2 / (1.0 - p2))
# 自由度＞２
else
yq = p_normal(ind) # 初期値計算のため
if ind[0] >= 0
x = 1.0 - 1.0 / (4.0 * df)
e = x * x - yq * yq / (2.0 * df)
if e > 0.5
t0 = yq / Math.sqrt(e)
else
x = Math.sqrt(df) / (pis * $p * df * gamma(df2, ind) / gamma(df2+0.5, ind))
t0 = x ** (1.0/df)
end
# ニュートン法
tt = newton(t0, 1.0e-6, 1.0e-10, 100, ind, &t_snx)
end
end
end
return tt
end
# 密度関数と分布関数の値
print("自由度は？ ")
$dof = Integer(gets())
print("目的とする結果は？ \n")
print(" =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n")
print(" =1 : ｐ％値（ P(X > u) = 0.01p となるuの値） ")
sw = Integer(gets())
z = Array.new(1)
if sw == 0
print("グラフ出力？（=1: yes, =0: no） ")
sw = Integer(gets())
if sw == 0
# 密度関数と分布関数の値
print(" データは？ ")
x = Float(gets())
y = t(x, z, $dof)
print("P(X = " + String(x) + ") = " + String(z[0]) + ", P( X < " + String(x) + ") = " + String(y) + "(自由度 = " + String($dof) + ")\n")
# グラフ出力
else
print(" 密度関数のファイル名は？ ")
file1 = gets().strip()
print(" 分布関数のファイル名は？ ")
file2 = gets().strip()
print(" データの下限は？ ")
x1 = Float(gets())
print(" データの上限は？ ")
x2 = Float(gets())
print(" 刻み幅は？ ")
h = Float(gets())
out1 = open(file1, "w")
out2 = open(file2, "w")
x = x1
while x < x2+0.5*h
y = t(x, z, $dof)
out1.print(String(x) + " " + String(z[0]) + "\n")
out2.print(String(x) + " " + String(y) + "\n")
x += h
end
out1.close()
out2.close()
end
# ％値
else
print("％の値は？ ")
x = Float(gets())
$p = 0.01 * x
if $p < 1.0e-7
print(String(x) + "％値 = ∞ (自由度 = " + String($dof) + ")\n")
elsif (1.0-$p)< 1.0e-7
print(String(x) + "％値 = -∞ (自由度 = " + String($dof) + ")\n")
else
ind = Array.new(1)
y = p_t(ind)
print(String(x) + "％値 = " + String(y) + " 収束 " + String(ind[0]) + " (自由度 = " + String($dof) + ")\n")
end
end