χ²分布

############################
# χ2分布の計算
#      coded by Y.Suganuma
############################

################################################
# 標準正規分布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 = 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
	end

	pis = Math.sqrt(Math::PI)
	chs = Math.sqrt(ch)
	x   = 0.5 * ch
	y   = Array.new(1)

	if df%2 != 0
		u  = Math.sqrt(x) * Math.exp(-x) / pis
		pp = 2.0 * (normal(chs, y) - 0.5)
		ia = 1

	else
		u  = x * Math.exp(-x)
		pp = 1.0 - Math.exp(-x)
		ia = 2
	end

	if ia != df
		i1 = ia
		while i1 < df-1
			pp -= 2.0 * u / i1
			u  *= (ch / i1)
			i1 += 2
		end
	end

	dd[0] = u / ch

	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）
#      ind : >=0 : normal（収束回数）
#            =-1 : 収束しなかった
##########################################
def p_chi(ind)

	chi_f = Proc.new { |sw, x|
		y = Array.new(1)
		z = chi(x, y, $dof)
		if sw == 0
			z - 1.0 + $p
		else
			z = y[0]
		end
	}

	xx = 0.0
			# 自由度＝１（正規分布を使用）
	if $dof == 1
		po  = $p
		$p *= 0.5
		x   = p_normal(ind)
		xx  = x * x
		$p  = po

	else
			# 自由度＝２
		if $dof == 2
			xx = -2.0 * Math.log($p)
			# 自由度＞２
		else
			x = p_normal(ind)   # 初期値計算のため

			if ind[0] >= 0
				w  = 2.0 / (9.0 * $dof)
				x  = 1.0 - w + x * Math.sqrt(w)
				x0 = (x ** 3.0) * $dof   # ニュートン法の初期値
				xx = newton(x0, 1.0e-6, 1.0e-10, 100, ind, &chi_f)
			end
		end
	end

	return xx
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 = chi(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 = chi(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) + "％値 = 0 (自由度 = " + String($dof) + ")\n")
	else
		ind = Array.new(1)
		y   = p_chi(ind)
		print(String(x) + "％値 = " + String(y) + "  収束 " + String(ind[0]) + " (自由度 = " + String($dof) + ")\n")
	end
end