二項分布

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

############
# 組合せ nCr
############

def comb(n, r) :

	c = 1.0
	x = 1.0
	y = 1.0

	if r > 0 :
		if r > n-r :
			r = n - r
		for i1 in range(n, n-r, -1) :
			x *= i1
		for i1 in range(2, r+1) :
			y *= i1
		c = x / y

	return c

######################################
# 二項分布の計算（P(X = x), P(X < x)）
#      x : データ
#      n,p : パラメータ
#      pr : P(X = x)
#      return : P(X < x)
#      coded by Y.Suganuma
######################################

def binomial(x, n, p, pr) :

	q     = 1.0 - p
	pr[0] = comb(n, x) * pow(p, x) * pow(q, n-x)
	f     = pr[0]
	for i1 in range(0, x) :
		f += comb(n, i1) * pow(p, i1) * pow(q, n-i1)

	return f

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

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

############################
# 二項分布の計算
#      coded by Y.Suganuma
###########################

s  = input("n は？ ")
n  = int(s)
s  = input("p は？ ")
p  = float(s)
s  = input("グラフ出力？（=1: yes,  =0: no） ")
sw = int(s)
pr = np.empty(1, np.float)
			# 密度関数と分布関数の値
if sw == 0 :
	s = input("   データは？ ")
	x = int(s)
	f = binomial(x, n, p, pr)
	print("P(X = " + str(x) + ") = " + str(pr[0]) + ",  P( X < " + str(x) + ") = " + str(f) + " (n = " + str(n) + ", p = " + str(p) + ")")
			# グラフ出力
else :
	file1 = input("   密度関数のファイル名は？ ")
	file2 = input("   分布関数のファイル名は？ ")
	out1  = open(file1, "w")
	out2  = open(file2, "w")
	for x in range(0, n+1) :
		f = binomial(x, n, p, pr)
		out1.write(str(x) + " " + str(pr[0]) + "\n")
		out2.write(str(x) + " " + str(f) + "\n")
	out1.close()
	out2.close()