指数分布

############################
# 指数分布の計算
#      coded by Y.Suganuma
############################

##########################################
# 指数分布の計算（P(X = x), P(X < x)）
#      x : データ
#      ram : 母数
#      pr : P(X = x)
#      return : P(X < x)
##########################################

def exponential(x, ram, pr)

	f = 0.0

	if x < 0
		pr[0] = 0.0
	else
		y     = Math.exp(-ram * x)
		pr[0] = ram * y
		f     = 1.0 - y
	end

	return f
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

#################################
# 関数値（f(x)）とその微分の計算
#################################
exp_snx = Proc.new { |sw, x|
	if sw == 0
		$p - Math.exp(-$ram * x)
	else
		$ram * Math.exp(-$ram * x)
	end
}

			# 密度関数と分布関数の値
print("母数は？ ")
$ram = 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())
		f = exponential(x, $ram, pr)
		print("P(X = " + String(x) + ") = " + String(pr[0]) + ",  P( X < " + String(x) + ") = " + String(f) + " (母数 = " + String($ram) + ")\n")
			# グラフ出力
	else
		print("   密度関数のファイル名は？ ")
		file1 = gets().strip()
		print("   分布関数のファイル名は？ ")
		file2 = gets().strip()
		print("   データの上限は？ ")
		up = Integer(gets())
		print("   刻み幅は？ ")
		h    = Float(gets())
		out1 = open(file1, "w")
		out2 = open(file2, "w")
		x    = 0
		while x < up+0.5*h
			f = exponential(x, $ram, 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) + "％値 = ∞ (母数 = " + String($ram) + ")\n")
	else
		ok = Array.new(1)
		f  = newton(0.0, 1.0e-6, 1.0e-10, 100, ok, &exp_snx)
		print(String(x) + "％値 = " + String(f) + " (母数 = " + String($ram) + ")\n")
	end
end