指数分布

# -*- coding: UTF-8 -*-
from math import *
import numpy as np

##########################################
# 指数分布の計算（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     = exp(-ram * x)
		pr[0] = ram * y
		f     = 1.0 - y

	return f

############################################
# 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 exponential, newton

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

############################
# P(X > x) - 1 + p（関数値）
############################
def exponential_f(x) :
	return p - exp(-ram * x)

############################
# P(X = x)（関数の微分）
############################
def exponential_df(x) :
	return ram * exp(-ram * x)

			# 密度関数と分布関数の値
s   = input("母数は？ ")
ram = float(s)
print("目的とする結果は？ ")
print("     =0 : 確率の計算（ P(X = x) 及び P(X < x) の値）")
s  = input("     =1 : ｐ％値（ P(X > u) = 0.01p となるuの値） ")
sw = int(s)
pr = 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)
		f = exponential(x, ram, pr)
		print("P(X = " + str(x) + ") = " + str(pr[0]) + ",  P( X < " + str(x) + ") = " + str(f) + " (母数 = " + str(ram) + ")")
			# グラフ出力
	else :
		file1 = input("   密度関数のファイル名は？ ")
		file2 = input("   分布関数のファイル名は？ ")
		s     = input("   データの上限は？ ")
		up    = float(s)
		s     = input("   刻み幅は？ ")
		h     = float(s)
		out1  = open(file1, "w")
		out2  = open(file2, "w")
		x     = 0
		while x < up+0.5*h :
			f = exponential(x, ram, pr)
			out1.write(str(x) + " " + str(pr[0]) + "\n")
			out2.write(str(x) + " " + str(f) + "\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(ram) + ")")
	else :
		ok = np.empty(1, np.int)
		f  = newton(exponential_f, exponential_df, 0.0, 1.0e-6, 1.0e-10, 100, ok)
		print(str(x) + "％値 = " + str(f) + " (母数 = " + str(ram) + ")")