Python Probability

import random
def random_face():
    return random.choice(["back","front"])

both_front = 0
first_front = 0
either_front = 0

random.seed()
for _ in range(10000):
    f = random_face()
    s = random_face()
    if f == "front":
        first_front += 1
    if f == "front" and s == "front":
        both_front += 1
    if f == "front" or s == "front":
        either_front += 1

print("P(both|first_front)", both_front / first_front)
print("P(both|first_front)", both_front / either_front)
>>P(both|first_front) 0.48842229199372056(~0.5)
>>P(both|first_front) 0.3275862068965517(~0.333..)

#다음은 균등 분포의 학률밀도함수이다.
def uniform_pdf(x):
    return 1 if x >= 0 and x < 1 else 0

"""
확률변수의 값이 특정 값보다 작거나 클 확률을 나타내는 함수는
누적분포함수(Cumulative distribution function, cdf)라고 한다.
"""

def uniform_cdf(x):
    if x < 0: return 0       # 균등분포의 확률은 0을 초과해야 함
    elif x < 1: return x     # 예시 P(X <= 0.4) = 0.4
    else: return 1           # 균등분포의 확률 < 1

from matplotlib import pyplot as plt
from math import sqrt, pi, exp, erf

def normal_pdf(x, mu=0, sigma=1):
    sqrt_two_pi = sqrt(2*pi)
    return (exp(-(x-mu) ** 2 / 2 / sigma ** 2) / (sqrt_two_pi * sigma))

xs = [x / 10.0 for x in range(-50, 50)]
plt.plot(xs, [normal_pdf(x,sigma=1) for x in xs],'-',label='mu=0,sigma=1 - standard')
plt.plot(xs, [normal_pdf(x,sigma=2) for x in xs],'--',label='mu=0,sigma=2')
plt.plot(xs, [normal_pdf(x,sigma=0.5) for x in xs],':',label='mu=0,sigma=0.5')
plt.plot(xs, [normal_pdf(x,mu=-1) for x in xs],'-.',label='mu=0,sigma=1')
plt.legend(prop={'size': 8})
plt.title("Various Normal PDFs")
plt.show()

def normal_cdf(x, mu=0, sigma=1):
    return (1 + erf((x - mu) / sqrt(2) / sigma)) / 2

xs = [x / 10.0 for x in range(-50, 50)]
plt.plot(xs,[normal_cdf(x,sigma=1) for x in xs],'-',label='mu=0, sigma=1')
plt.plot(xs,[normal_cdf(x,sigma=2) for x in xs],'--',label='mu=0, sigma=2')
plt.plot(xs,[normal_cdf(x,sigma=0.5) for x in xs],':',label='mu=0, sigma=0.5')
plt.plot(xs,[normal_cdf(x,mu=-1) for x in xs],'-.',label='mu=-1, sigma=1')
plt.legend(loc=2,prop={'size':8}) # location 1 2 3 4 택1
plt.title("Various Normal CDFs")
plt.show()

def inverse_normal_cdf(p, mu=0, sigma=1, tolerance=0.00001):
    """이진 검색을 통해 역함수 근사"""

    # 표준정규분포가 아니라면 표준정규분포로 변환
    if mu != 0 or sigma != 1:
        return mu + sigma * inverse_normal_cdf(p, tolerance=tolerance)

    low_z, low_p = -10.0, 0     # normal_cdf(-10)는 0에 근접
    hi_z, hi_p = 10.0, 1        # normal_cdf(10)는 1에 근접
    while hi_z - low_z > tolerance:
        mid_z = (low_z + hi_z) / 2
        mid_p = normal_cdf(mid_z)
        if mid_p < p:
            # 중간 값이 작으면 더 큰 값 검색
            low_z, low_p = mid_z, mid_p
        elif mid_p > p:
            hi_z, hi_p = mid_z, mid_p
        else:
            break

    return mid_z

def bernoulli_trial(p):
    return 1 if random.random() < p else 0

def binomial(n, p):
    return sum(bernoulli_trial(p) for _ in range(n))

# n이 적당히 클 때 mu는 np
binomial(10000, 0.5)
>>4996(~5000)

from collections import Counter
def make_hist(p, n, num_points):
    data = [binomial(n, p) for _ in range(num_points)]
    histogram = Counter(data)
    plt.bar([x - 0.4 for x in histogram.keys()],
            [v / num_points for v in histogram.values()],
            0.8,
            color='0.75')
    mu = p * n
    sigma = sqrt(n * p * (1 - p))

    # 근사된 정규분포를 라인 차트로 표현
    xs = range(min(data), max(data) + 1)
    ys = [normal_cdf(i + 0.5, mu, sigma) -  normal_cdf(i - 0.5, mu, sigma)
         for i in xs]
    plt.plot(xs,ys)
    plt.title("Binomial Distribution & Normal Approximation")
    plt.show()

make_hist(0.5, 100, 10000)