正規分布

# -*- 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

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

----------------------------------

# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import normal, bisection

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

			# 密度関数と分布関数の値
s    = input("平均値は？ ")
mean = float(s)
s    = input("標準偏差は？ ")
sd   = 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)
		xx = (x - mean) / sd
		f  = normal(xx, pr)
		print("P(X = " + str(x) + ") = " + str(pr[0]) + ",  P( X < " + str(x) + ") = " + str(f))
			# グラフ出力
	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 :
			xx = (x - mean) / sd
			f  = normal(xx, 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) + "％値 = ∞")
	elif (1.0-p)< 1.0e-7 :
		print(str(x) + "％値 = -∞")
	else :
		ok = np.empty(1, np.int)
		y  = sd * p_normal(ok) + mean
		print(str(x) + "％値 = " + str(y) + "  収束 " + str(ok[0]))