# -*- coding: UTF-8 -*-
from math import *
import numpy as np
################################################
# 標準正規分布N(0,1)の計算（P(X = x), P(X < x)）
# w : P(X = x)
# return : P(X < x)
################################################
def normal(x, w) :
# 確率密度関数（定義式）
w[0] = exp(-0.5 * x * x) / sqrt(2.0*pi)
# 確率分布関数（近似式を使用）
y = 0.70710678118654 * abs(x)
z = 1.0 + y * (0.0705230784 + y * (0.0422820123 + y * (0.0092705272 + y * (0.0001520143 + y * (0.0002765672 + y * 0.0000430638)))))
P = 1.0 - z ** (-16.0)
if x < 0.0 :
P = 0.5 - 0.5 * P
else :
P = 0.5 + 0.5 * P
return P
##########################################
# χ２乗分布の計算（P(X = ch), P(X < ch)）
# dd : P(X = ch)
# df : 自由度
# return : P(X < ch)
##########################################
def chi(ch, dd, df) :
if ch < 1.0e-10 :
ch = 1.0e-10
pis = sqrt(pi)
chs = sqrt(ch)
x = 0.5 * ch
y = np.empty(1, np.float)
if df%2 != 0 :
u = sqrt(x) * exp(-x) / pis
pp = 2.0 * (normal(chs, y) - 0.5)
ia = 1
else :
u = x * exp(-x)
pp = 1.0 - exp(-x)
ia = 2
if ia != df :
for i1 in range(ia, df-1, 2) :
pp -= 2.0 * u / i1
u *= (ch / i1)
dd[0] = u / ch
return pp
############################################
# 二分法による非線形方程式(f(x)=0)の解
# fn : f(x)を計算する関数名
# x1,x2 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜f(x(k))｜＜eps2）
# max : 最大試行回数
# ind : 実際の試行回数
# （負の時は解を得ることができなかった）
# return : 解
# coded by Y.Suganuma
############################################
def bisection(fn, x1, x2, eps1, eps2, max, ind) :
x0 = 0.0
f1 = fn(x1)
f2 = fn(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
else :
sw = 0
while sw == 0 and ind[0] >= 0 :
sw = 1
ind[0] += 1
x0 = 0.5 * (x1 + x2)
f0 = fn(x0)
if abs(f0) > eps2 :
if ind[0] <= max :
if abs(x1-x2) > eps1 and abs(x1-x2) > eps1*abs(x2) :
sw = 0
if f0*f1 < 0.0 :
x2 = x0
f2 = f0
else :
x1 = x0
f1 = f0
else :
ind[0] = -1
return x0
############################################
# Newton法による非線形方程式(f(x)=0)の解
# fn : f(x)を計算する関数名
# dfn : f(x)の微分を計算する関数名
# x0 : 初期値
# eps1 : 終了条件１（｜x(k+1)-x(k)｜＜eps1）
# eps2 : 終了条件２（｜f(x(k))｜＜eps2）
# max : 最大試行回数
# ind : 実際の試行回数
# （負の時は解を得ることができなかった）
# return : 解
# coded by Y.Suganuma
############################################
def newton(fn, dfn, x0, eps1, eps2, max, ind) :
x1 = x0
x = x1
ind[0] = 0
sw = 0
while sw == 0 and ind[0] >= 0 :
sw = 1
ind[0] += 1
g = fn(x1)
if abs(g) > eps2 :
if ind[0] <= max :
dg = dfn(x1)
if abs(dg) > eps2 :
x = x1 - g / dg
if abs(x-x1) > eps1 and abs(x-x1) > eps1*abs(x) :
x1 = x
sw = 0
else :
ind[0] = -1
else :
ind[0] = -1
return x
----------------------------------
# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import chi, normal, bisection, newton
############################
# χ2分布の計算
# coded by Y.Suganuma
############################
##########################################
# 1.0 - p - P(X>x)（関数値, 標準正規分布）
##########################################
def normal_f(x) :
y = np.empty(1, np.float)
return 1.0 - p - normal(x, y)
################################################################
# 標準正規分布N(0,1)のｐ％値（P(X > u) = 0.01p）（二分法を使用）
# ind : >= 0 : normal（収束回数）
# = -1 : 収束しなかった
################################################################
def p_normal(ind) :
u = bisection(normal_f, -7.0, 7.0, 1.0e-6, 1.0e-10, 100, ind)
return u
##########################################
# 1.0 - p - P(X > x)（関数値, χ２乗分布）
##########################################
def chi_f(x) :
y = np.empty(1, np.float)
return (chi(x, y, dof) - 1.0 + p)
####################################
# P(X = x)（関数の微分, χ２乗分布）
####################################
def chi_df(x) :
y = np.empty(1, np.float)
z = chi(x, y, dof)
return y[0]
##########################################
# χ２乗分布のｐ％値（P(X > u) = 0.01p）
# ind : >=0 : normal（収束回数）
# =-1 : 収束しなかった
##########################################
def p_chi(ind) :
xx = 0.0
# 自由度＝１（正規分布を使用）
if dof == 1 :
global p # この宣言がないとエラー
po = p
p *= 0.5
x = p_normal(ind)
xx = x * x
p = po
else :
# 自由度＝２
if dof == 2 :
xx = -2.0 * log(p)
# 自由度＞２
else :
x = p_normal(ind) # 初期値計算のため
if ind[0] >= 0 :
w = 2.0 / (9.0 * dof)
x = 1.0 - w + x * sqrt(w)
x0 = (x ** 3.0) * dof # ニュートン法の初期値
xx = newton(chi_f, chi_df, x0, 1.0e-6, 1.0e-10, 100, ind)
return xx
# 密度関数と分布関数の値
s = input("自由度は？ ")
dof = int(s)
print("目的とする結果は？ ")
print(" =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）")
s = input(" =1 : ｐ％値（ P(X > u) = 0.01p となるuの値） ")
sw = int(s)
z = np.empty(1, np.float)
if sw == 0 :
s = input("グラフ出力？（=1: yes, =0: no） ")
sw = int(s)
if sw == 0 :
# 密度関数と分布関数の値
s = input(" データは？ ")
x = float(s)
y = chi(x, z, dof)
print("P(X = " + str(x) + ") = " + str(z[0]) + ", P( X < " + str(x) + ") = " + str(y) + "(自由度 = " + str(dof) + ")")
# グラフ出力
else :
file1 = input(" 密度関数のファイル名は？ ")
file2 = input(" 分布関数のファイル名は？ ")
s = input(" データの下限は？ ")
x1 = float(s)
s = input(" データの上限は？ ")
x2 = float(s)
s = input(" 刻み幅は？ ")
h = float(s)
out1 = open(file1, "w")
out2 = open(file2, "w")
x = x1
while x < x2+0.5*h :
y = chi(x, z, dof)
out1.write(str(x) + " " + str(z[0]) + "\n")
out2.write(str(x) + " " + str(y) + "\n")
x += h
out1.close()
out2.close()
# ％値
else :
s = input("％の値は？ ")
x = float(s)
p = 0.01 * x
if p < 1.0e-7 :
print(str(x) + "％値 = ∞ (自由度 = " + str(dof) + ")")
elif (1.0-p)< 1.0e-7 :
print(str(x) + "％値 = 0 (自由度 = " + str(dof) + ")")
else :
ind = np.empty(1, np.int)
y = p_chi(ind)
print(str(x) + "％値 = " + str(y) + " 収束 " + str(ind[0]) + " (自由度 = " + str(dof) + ")")