正規分布

############################
# 正規分布の計算
#      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

############################################
# 二分法による非線形方程式(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

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

			# 密度関数と分布関数の値
print("平均値は？ ")
mean = Float(gets())
print("標準偏差は？ ")
sd = Float(gets())
print("目的とする結果は？\n")
print("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）\n")
print("     =1 : ｐ％値（ P(X > u) = 0.01p となるuの値） ")
sw = Integer(gets())
pr = Array.new(1)

if sw == 0
	print("グラフ出力？（=1: yes,  =0: no） ")
	sw = Integer(gets())
	if sw == 0
			# 密度関数と分布関数の値
		print("   データは？ ")
		x  = Float(gets())
		xx = (x - mean) / sd
		f  = normal(xx, pr)
		print("P(X = " + String(x) + ") = " + String(pr[0]) + ",  P( X < " + String(x) + ") = " + String(f) + "\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
			xx = (x - mean) / sd
			f  = normal(xx, pr)
			out1.print(String(x) + " " + String(pr[0]) + "\n")
			out2.print(String(x) + " " + String(f) + "\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) + "％値 = ∞\n")
	elsif (1.0-$p)< 1.0e-7
		print(String(x) + "％値 = -∞\n")
	else
		ind = Array.new(1)
		y   = sd * p_normal(ind) + mean
		print(String(x) + "％値 = " + String(y) + "  収束 " + String(ind[0]) + "\n")
	end
end