정규분포¶

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats

mu, sigma = 0, 1                  # 표준정규분포
# 평균에서 좌우 3 표준편차. 확률의 99.7%가 이 안에 있다.
x = np.linspace(mu - 3*sigma, mu + 3*sigma, 200)

fig, ax = plt.subplots(figsize=(12, 3))
# 두 함수를 같은 축에 겹쳐 관계를 본다.
#   PDF는 평균에서 가장 높고 좌우로 떨어진다.
#   CDF는 0에서 1로 단조 증가하며, PDF가 가장 높은 곳에서 가장 가파르다.
# CDF의 기울기가 곧 PDF이기 때문이다.
ax.plot(x, stats.norm(mu, sigma).pdf(x), label='PDF')
ax.plot(x, stats.norm(mu, sigma).cdf(x), label='CDF')
ax.spines[['top', 'right']].set_visible(False)
ax.legend()
plt.show()

import numpy as np
from scipy import integrate, optimize, stats

phi, Phi = stats.norm.pdf, stats.norm.cdf

# 봉우리 높이가 1/sqrt(2pi) 인가.
print(f"phi(0)        = {phi(0):.9f}")
print(f"1/sqrt(2*pi)  = {1 / np.sqrt(2 * np.pi):.9f}")

# phi'' = (x^2 - 1) phi 를 중심차분으로 맞춰 본다. 부호가 |x|=1 에서 바뀐다.
h = 1e-5
print(f"{'x':>6}{'phi_xx':>14}{'(x^2-1)phi':>14}")
for x in (-2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0):
    num = (phi(x + h) - 2 * phi(x) + phi(x - h)) / h**2
    print(f"{x:>6.1f}{num:>14.7f}{(x*x - 1) * phi(x):>14.7f}")

# 변곡점의 높이는 봉우리의 exp(-1/2) 배다.
print(f"phi(1)/phi(0) = {phi(1) / phi(0):.6f},  exp(-1/2) = {np.exp(-0.5):.6f}")

# 두 곡선이 만나는 곳: Phi(x) = phi(x).
x_star = optimize.brentq(lambda x: Phi(x) - phi(x), -3, 0)
print(f"교차점 x* = {x_star:.6f},  공통값 = {phi(x_star):.6f}")

# 교차가 한 번뿐인지 격자에서 부호변화를 세어 확인한다.
g = np.linspace(-3, 3, 600_001)
d = Phi(g) - phi(g)
print(f"[-3,3] 에서 부호변화 횟수 = {int(np.sum(np.sign(d[:-1]) != np.sign(d[1:])))}")

# 그림이 쓰는 세로 범위. CDF 는 0~1 을 다 쓰고 PDF 는 0.4 아래에 갇힌다.
print(f"PDF 의 최댓값 = {phi(0):.4f},  CDF 의 범위 = {Phi(-3):.4f} ~ {Phi(3):.4f}")
print(f"phi(3)/phi(0) = {phi(3) / phi(0):.4f}   <- 그림 양끝에서 밀도는 봉우리의 1% 남짓")
print(f"quad 로 잰 전체 질량 = {integrate.quad(phi, -np.inf, np.inf)[0]:.12f}")

phi(0)        = 0.398942280
1/sqrt(2*pi)  = 0.398942280
     x        phi_xx    (x^2-1)phi
  -2.0     0.1619729     0.1619729
  -1.0    -0.0000003     0.0000000
  -0.5    -0.2640488    -0.2640490
   0.0    -0.3989420    -0.3989423
   0.5    -0.2640488    -0.2640490
   1.0    -0.0000003     0.0000000
   2.0     0.1619729     0.1619729
phi(1)/phi(0) = 0.606531,  exp(-1/2) = 0.606531
교차점 x* = -0.302631,  공통값 = 0.381086
[-3,3] 에서 부호변화 횟수 = 1
PDF 의 최댓값 = 0.3989,  CDF 의 범위 = 0.0013 ~ 0.9987
phi(3)/phi(0) = 0.0111   <- 그림 양끝에서 밀도는 봉우리의 1% 남짓
quad 로 잰 전체 질량 = 1.000000000000

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

mu, sigma = 1, 2
x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 400)

dist = stats.norm(loc=mu, scale=sigma)
y_cdf = dist.cdf(x)
y_pdf = dist.pdf(x)

fig, ax_cdf = plt.subplots(figsize=(12, 3))

# CDF는 왼쪽 축(0~1). PDF는 오른쪽 축(밀도).
# 두 함수의 눈금 규모가 달라 한 축에 그리면 한쪽이 납작해지므로 축을 나눈다.
# 이 절에서는 두 축의 관계가 고정되어 있어(CDF는 PDF의 적분) 안전한 사용이다.
ax_cdf.plot(x, y_cdf, lw=2, label="CDF P(X ≤ x)")
ax_cdf.set_xlabel("x")
ax_cdf.set_ylabel("P(X ≤ x)")
ax_cdf.set_ylim(-0.02, 1.02)

# 기준점 세 개를 표시한다: 평균에서 -1, 0, +1 표준편차.
# CDF 값이 각각 약 0.159, 0.500, 0.841 이 나온다.
# 0.841 - 0.159 = 0.682 가 곧 "68% 규칙"이다.
for xv in [mu - sigma, mu, mu + sigma]:
    yv = dist.cdf(xv)
    ax_cdf.axvline(xv, linestyle='--', color='gray', alpha=0.7)
    ax_cdf.text(xv, yv + 0.05, f"P(X≤{xv:.0f})={yv:.3f}",
                ha='center', fontsize=9)

# 오른쪽 축에 밀도함수
ax_pdf = ax_cdf.twinx()
ax_pdf.plot(x, y_pdf, lw=2, color='tab:red', label="PDF (density)")
ax_pdf.set_ylabel("Density", color='tab:red')

ax_cdf.set_title(f"Normal({mu}, {sigma}) — CDF with PDF Overlay")
plt.tight_layout()
plt.show()

import numpy as np
from scipy import optimize, stats

mu, sigma = 1, 2
dist = stats.norm(loc=mu, scale=sigma)

# (1) F 가 가장 가파른 자리를 격자로 찾아 평균과 맞춰 본다.
g = np.linspace(mu - 3*sigma, mu + 3*sigma, 600_001)
slope = np.gradient(dist.cdf(g), g)
print(f"기울기가 최대인 자리 = {g[slope.argmax()]:.5f}   (mu = {mu})")
print(f"그때의 기울기        = {slope.max():.6f}")
print(f"f(mu) = 1/(sigma*sqrt(2*pi)) = {1 / (sigma * np.sqrt(2 * np.pi)):.6f}")

# F'' = 0 인 자리, 곧 밀도의 봉우리. 중심차분으로 F'' 를 재어 영점을 찾는다.
h = 1e-4
F2 = lambda x: (dist.cdf(x + h) - 2 * dist.cdf(x) + dist.cdf(x - h)) / h**2
print(f"F'' 의 영점 = {optimize.brentq(F2, mu - sigma, mu + sigma):.6f}   (변곡점은 평균 하나)")

# (2) 그림이 표시한 세 자리의 CDF 값.
print()
for xv in (mu - sigma, mu, mu + sigma):
    print(f"F({xv:+.0f}) = {dist.cdf(xv):.6f}   Phi({(xv - mu)/sigma:+.0f}) = "
          f"{stats.norm.cdf((xv - mu)/sigma):.6f}")
print(f"F(mu+sigma) - F(mu-sigma) = {dist.cdf(mu+sigma) - dist.cdf(mu-sigma):.6f}"
      f"   <- 68% 규칙의 참값")

# 모수와 무관하다는 것을 네 쌍에서 확인한다.
print()
print(f"{'mu':>6}{'sigma':>8}{'F(mu-sigma)':>14}{'F(mu+sigma)':>14}")
for m, s in ((0, 1), (1, 2), (-5, 0.3), (100, 15)):
    d = stats.norm(loc=m, scale=s)
    print(f"{m:>6}{s:>8}{d.cdf(m - s):>14.9f}{d.cdf(m + s):>14.9f}")

기울기가 최대인 자리 = 1.00000   (mu = 1)
그때의 기울기        = 0.199471
f(mu) = 1/(sigma*sqrt(2*pi)) = 0.199471
F'' 의 영점 = 1.000000   (변곡점은 평균 하나)

F(-1) = 0.158655   Phi(-1) = 0.158655
F(+1) = 0.500000   Phi(+0) = 0.500000
F(+3) = 0.841345   Phi(+1) = 0.841345
F(mu+sigma) - F(mu-sigma) = 0.682689   <- 68% 규칙의 참값

    mu   sigma   F(mu-sigma)   F(mu+sigma)
     0       1   0.158655254   0.841344746
     1       2   0.158655254   0.841344746
    -5     0.3   0.158655254   0.841344746
   100      15   0.158655254   0.841344746

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

mu, sigma = 0, 1
prob = 0.975      # 95% 신뢰구간의 한쪽 끝. 양쪽 꼬리에 2.5%씩 남긴다.

dist = stats.norm(loc=mu, scale=sigma)
x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 1000)
pdf = dist.pdf(x)

# ppf는 CDF의 역함수다. "누적확률이 이만큼 되는 지점은 어디인가"에 답한다.
#   cdf: 값 -> 확률
#   ppf: 확률 -> 값
# ppf(0.975)가 그 유명한 1.96 이며, 신뢰구간 공식의 z값이 여기서 나온다.
z = dist.ppf(prob)

fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(x, pdf, color='b', lw=2, label='PDF')
ax.plot([z, z], [0, dist.pdf(z)], color='k', lw=3)   # 경계선
# 왼쪽 97.5%를 칠한다. 칠해진 넓이가 곧 확률이라는 점이 요점이다.
ax.fill_between(x[x <= z], pdf[x <= z], 0,
                interpolate=True, color='r', alpha=0.25,
                label=f"P(X ≤ {z:.2f}) = {prob}")
ax.text(z + 0.05, dist.pdf(z) / 2,
        f"ppf({prob}) = {z:.4f}", fontsize=11, va='center')
ax.set_title(f"Normal({mu}, {sigma}) — PPF (Quantile Function)")
ax.legend(loc='upper left', frameon=False)
plt.tight_layout()
plt.show()

import numpy as np
from scipy import integrate, stats

ppf, cdf, pdf = stats.norm.ppf, stats.norm.cdf, stats.norm.pdf

# (2) z_{0.975} 와 어림수 1.96 의 신뢰수준.
z975 = ppf(0.975)
print(f"z_0.975              = {z975:.9f}")
print(f"Phi(1.96)            = {cdf(1.96):.9f}")
print(f"1.96 의 양측 꼬리확률 = {2 * (1 - cdf(1.96)):.9f}")
print(f"1.96 의 신뢰수준      = {100 * (2 * cdf(1.96) - 1):.6f}%")
print(f"quad 로 잰 Phi(z975) = {integrate.quad(pdf, -np.inf, z975)[0]:.12f}")

# (1) 분위수함수가 역함수임을 왕복으로 확인한다. 오차가 기계오차 수준이어야 한다.
print()
print("ppf(cdf(x)) - x :", [f"{ppf(cdf(x)) - x:.1e}" for x in (-3, -1, 0, 1, 2.5)])
print("ppf(1-p) + ppf(p):", [f"{ppf(1 - p) + ppf(p):.1e}" for p in (0.9, 0.975, 0.999)])

# (1) 도함수 1/phi(z) 가 꼬리에서 터지는 모습. 중심차분과 닫힌 꼴을 나란히 둔다.
print()
print(f"{'p':>9}{'z=ppf(p)':>12}{'1/phi(z)':>14}{'수치 도함수':>16}")
h = 1e-7
for p in (0.5, 0.75, 0.9, 0.975, 0.999, 0.99999):
    num = (ppf(p + h) - ppf(p - h)) / (2 * h)
    print(f"{p:>9}{ppf(p):>12.6f}{1 / pdf(ppf(p)):>14.4f}{num:>16.4f}")

# 쪽의 분위수 표를 다시 계산한다.
print()
for q in (0.500, 0.900, 0.950, 0.975, 0.995):
    print(f"ppf({q:.3f}) = {ppf(q):.6f}")

z_0.975              = 1.959963985
Phi(1.96)            = 0.975002105
1.96 의 양측 꼬리확률 = 0.049995790
1.96 의 신뢰수준      = 95.000421%
quad 로 잰 Phi(z975) = 0.975000000000

ppf(cdf(x)) - x : ['-4.4e-16', '1.1e-16', '0.0e+00', '-1.1e-16', '-1.3e-15']
ppf(1-p) + ppf(p): ['0.0e+00', '0.0e+00', '0.0e+00']

        p    z=ppf(p)      1/phi(z)          수치 도함수
      0.5    0.000000        2.5066          2.5066
     0.75    0.674490        3.1469          3.1469
      0.9    1.281552        5.6981          5.6981
    0.975    1.959964       17.1101         17.1101
    0.999    3.090232      296.9924        296.9924
  0.99999    4.264891    22327.7432      22328.4367

ppf(0.500) = 0.000000
ppf(0.900) = 1.281552
ppf(0.950) = 1.644854
ppf(0.975) = 1.959964
ppf(0.995) = 2.575829

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

mu, sigma = 0, 1
dist = stats.norm(loc=mu, scale=sigma)

x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 400)
# 생존함수 SF(x) = P(X > x) = 1 - CDF(x).
# 수학적으로는 CDF의 여집합일 뿐이지만, 계산 방식이 다르다.
# scipy는 sf를 1-cdf로 계산하지 않고 꼬리를 직접 적분하므로
# 아주 작은 확률에서도 정밀도를 잃지 않는다(아래 절 참고).
cdf = dist.cdf(x)
sf = dist.sf(x)

fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(x, cdf, lw=2, label='CDF  P(X ≤ x)')
ax.plot(x, sf, lw=2, label='SF   P(X > x)')
# 두 곡선이 만나는 지점을 표시한다.
# 대칭분포에서는 평균에서 CDF = SF = 0.5 로 교차한다.
ax.axvline(0, ls=':', color='gray', alpha=0.6)
ax.axhline(0.5, ls=':', color='gray', alpha=0.6)
ax.annotate("CDF + SF = 1", xy=(1.2, 0.5), fontsize=12,
            bbox=dict(boxstyle='round,pad=0.3', fc='lightyellow', ec='gray'))
ax.set_xlabel('x')
ax.set_ylabel('Probability')
ax.set_ylim(-0.03, 1.03)
ax.legend(loc='center left', frameon=False)
ax.set_title(f"Normal({mu}, {sigma}) — CDF vs Survival Function")
ax.grid(True, linestyle=':', alpha=0.5)
plt.tight_layout()
plt.show()

import numpy as np
from scipy import optimize, stats

sf, cdf, pdf = stats.norm.sf, stats.norm.cdf, stats.norm.pdf

# (1) S(x) = Phi(-x) 를 격자 전체에서 확인한다. 비트까지 같은지도 본다.
g = np.linspace(-6, 6, 2001)
print(f"max |S(x) - Phi(-x)| = {max(abs(sf(x) - cdf(-x)) for x in g):.3e}")
print(f"비트까지 같은가       = {all(sf(x) == cdf(-x) for x in g)}")

# (1) 교차점 Phi(x) = S(x).
print(f"교차점 = {optimize.brentq(lambda x: cdf(x) - sf(x), -2, 2):.6f},"
      f"  공통값 = {cdf(0):.6f}")

# (2) 비 S(x+1)/S(x) 와 점근식. 지수분포 Exp(1) 의 비는 상수다.
print()
print(f"{'x':>4}{'S(x)':>13}{'S(x+1)/S(x)':>14}{'점근식':>12}{'Exp(1) 의 비':>14}")
for x in (1, 2, 3, 4, 5):
    ratio = sf(x + 1) / sf(x)
    approx = x / (x + 1) * np.exp(-x - 0.5)
    print(f"{x:>4}{sf(x):>13.4e}{ratio:>14.6f}{approx:>12.6f}{np.exp(-1):>14.6f}")

# 그림이 쓰는 세로축에서 꼬리가 차지하는 몫.
print()
print(f"S(3) = {sf(3):.9f}  ->  세로축 눈금(1.06)의 {sf(3) / 1.06 * 100:.4f}%")
print(f"S(6) = {sf(6):.4e}  ->  세로축 눈금의 {sf(6) / 1.06 * 100:.2e}%")

max |S(x) - Phi(-x)| = 0.000e+00
비트까지 같은가       = True
교차점 = 0.000000,  공통값 = 0.500000

   x         S(x)   S(x+1)/S(x)         점근식    Exp(1) 의 비
   1   1.5866e-01      0.143393    0.111565      0.367879
   2   2.2750e-02      0.059336    0.054723      0.367879
   3   1.3499e-03      0.023462    0.022648      0.367879
   4   3.1671e-05      0.009051    0.008887      0.367879
   5   2.8665e-07      0.003442    0.003406      0.367879

S(3) = 0.001349898  ->  세로축 눈금(1.06)의 0.1273%
S(6) = 9.8659e-10  ->  세로축 눈금의 9.31e-08%

from scipy import stats

# 꼬리 확률을 두 가지 방법으로 구해 비교한다.
#   나쁜 방법: 1 - CDF.  CDF가 1에 아주 가까우면 뺄셈에서 유효숫자가 날아간다(상쇄).
#   좋은 방법: SF.       꼬리를 직접 계산하므로 상쇄가 일어나지 않는다.
for x in (6, 8, 10, 12):
    bad = 1 - stats.norm.cdf(x)
    good = stats.norm.sf(x)
    print(f"x={x:>3}:  1-cdf = {bad:.6e}   sf = {good:.6e}")

x=  6:  1-cdf = 9.865877e-10   sf = 9.865876e-10
x=  8:  1-cdf = 6.661338e-16   sf = 6.220961e-16
x= 10:  1-cdf = 0.000000e+00   sf = 7.619853e-24
x= 12:  1-cdf = 0.000000e+00   sf = 1.776482e-33

import numpy as np
from scipy import optimize, stats

sf, cdf = stats.norm.sf, stats.norm.cdf
eps = np.finfo(float).eps        # = 2^-52

print(f"eps          = {eps!r}")
print(f"1 의 아래 이웃 = {np.nextafter(1.0, 0.0)!r}")
print(f"그 간격        = {1 - np.nextafter(1.0, 0.0)!r}   (= 2^-53)")
print(f"접히는 문턱    = {2.0**-54!r}   (= 2^-54)")

# (1) 의 예측: S(x) = 2^-54 인 자리.
x_pred = stats.norm.isf(2.0**-54)

# 실제 경계: (1 - cdf(x)) == 0 이 되는 가장 작은 double 을 이분법으로 찾는다.
lo, hi = 8.0, 9.0
while hi != np.nextafter(lo, hi):
    mid = (lo + hi) / 2
    if (1 - cdf(mid)) == 0.0:
        hi = mid
    else:
        lo = mid

print()
print(f"예측 경계 isf(2^-54) = {x_pred!r}")
print(f"실제 경계            = {hi!r}")
print(f"차이                 = {hi - x_pred!r}")

# sf 가 언더플로하는 자리도 같은 방법으로 찾는다.
lo2, hi2 = 30.0, 45.0
while hi2 != np.nextafter(lo2, hi2):
    mid = (lo2 + hi2) / 2
    if sf(mid) == 0.0:
        hi2 = mid
    else:
        lo2 = mid
print(f"sf 가 0 이 되는 경계 = {hi2:.6f}   (그 직전 sf = {sf(lo2):.3e})")
print(f"배정밀도 최소 양수   = {5e-324!r}")
print(f"S(x) = 5e-324 가 되는 자리 = "
      f"{optimize.brentq(lambda x: stats.norm.logsf(x) - np.log(5e-324), 30, 50):.6f}")
print(f"sf(38) = {sf(38)!r},  logsf(38) = {stats.norm.logsf(38):.6f},"
      f"  exp(logsf(38)) = {np.exp(stats.norm.logsf(38)):.3e}")

# 상대오차가 2^-54/S(x) 를 따라 커지는지 표로 본다.
print()
print(f"{'x':>6}{'sf(x)':>13}{'1-cdf 의 상대오차':>20}{'2^-54/S(x)':>14}")
for x in (2, 4, 6, 7, 8, 8.2, 8.3, 10):
    good = sf(x)
    rel = abs((1 - cdf(x)) - good) / good
    print(f"{x:>6}{good:>13.3e}{rel:>20.2e}{2.0**-54 / good:>14.2e}")

eps          = 2.220446049250313e-16
1 의 아래 이웃 = 0.9999999999999999
그 간격        = 1.1102230246251565e-16   (= 2^-53)
접히는 문턱    = 5.551115123125783e-17   (= 2^-54)

예측 경계 isf(2^-54) = 8.292361075813597
실제 경계            = 8.292361075813597
차이                 = 0.0
sf 가 0 이 되는 경계 = 37.677121   (그 직전 sf = 5.886e-311)
배정밀도 최소 양수   = 5e-324
S(x) = 5e-324 가 되는 자리 = 38.467406
sf(38) = 0.0,  logsf(38) = -726.557216,  exp(logsf(38)) = 2.885e-316

     x        sf(x)        1-cdf 의 상대오차    2^-54/S(x)
     2    2.275e-02            4.58e-16      2.44e-15
     4    3.167e-05            3.64e-15      1.75e-12
     6    9.866e-10            5.61e-08      5.63e-08
     7    1.280e-12            4.11e-05      4.34e-05
     8    6.221e-16            7.08e-02      8.92e-02
   8.2    1.202e-16            7.63e-02      4.62e-01
   8.3    5.206e-17            1.00e+00      1.07e+00
    10    7.620e-24            1.00e+00      7.29e+06

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)
data = stats.norm(loc=0, scale=1).rvs(10_000)     # 참 모수는 (0, 1)

fig, ax = plt.subplots(figsize=(12, 3))
# density=True 로 넓이를 1로 맞춰야 밀도곡선과 같은 눈금에 놓인다
_, bins, _ = ax.hist(data, bins=100, density=True, color='blue', alpha=0.7, label="Samples")
# 참 모수 (0, 1)이 아니라 **표본에서 추정한** 평균과 표준편차로 곡선을 그린다.
# 실제 분석에서는 참값을 모르기 때문이다.
ax.plot(bins, stats.norm(data.mean(), data.std()).pdf(bins),
        '--r', lw=3, label="Estimated Normal PDF")
ax.legend()
plt.show()

import numpy as np
from scipy import special, stats

np.random.seed(0)
data = stats.norm(loc=0, scale=1).rvs(10_000)
n = len(data)

# (1) 두 추정값과 그 이론 기댓값·표준오차.
#     numpy 의 .std() 는 ddof=0 이 기본이므로 최대가능도추정값이다.
c4 = np.exp(0.5 * np.log(2 / (n - 1))
            + special.gammaln(n / 2) - special.gammaln((n - 1) / 2))
E_sd = c4 * np.sqrt((n - 1) / n)        # E[std(ddof=0)]

print(f"n = {n}")
print(f"mu-hat  = {data.mean():+.6f}   이론 E = {0:+.6f}   이론 SE = {1/np.sqrt(n):.6f}"
      f"   z = {data.mean()*np.sqrt(n):+.3f}")
print(f"sig-hat = {data.std():.6f}   이론 E = {E_sd:.6f}   이론 SE = {1/np.sqrt(2*n):.6f}"
      f"   z = {(data.std()-E_sd)/(1/np.sqrt(2*n)):+.3f}")
print(f"sig-hat^2 = {data.var():.6f}  이론 E = (n-1)/n = {(n-1)/n:.6f}"
      f"   이론 SE = {np.sqrt(2*(n-1))/n:.6f}")

# 적합곡선의 봉우리가 참 봉우리보다 얼마나 높은가.
print(f"적합 봉우리 = {1/(data.std()*np.sqrt(2*np.pi)):.6f},"
      f"  참 봉우리 = {1/np.sqrt(2*np.pi):.6f},"
      f"  차이 = {100*(1/data.std() - 1):+.2f}%")

# (2) 밀도 히스토그램 막대의 이론 기댓값과 표준오차.
counts, edges = np.histogram(data, bins=100)
D = edges[1] - edges[0]
height = counts / (n * D)
p = stats.norm.cdf(edges[1:]) - stats.norm.cdf(edges[:-1])   # 칸 확률
exp_h = p / D                                                # 칸 평균 밀도
se_h = np.sqrt(p * (1 - p) / n) / D                          # 이항 표준오차

print()
print(f"칸 폭 Delta = {D:.6f},  범위 = [{edges[0]:.4f}, {edges[-1]:.4f}]")
print(f"{'칸 중앙':>10}{'관측 높이':>12}{'이론 E':>11}{'이론 SE':>10}{'z':>8}")
for i in (10, 30, 45, 50, 55, 70, 90):
    mid = (edges[i] + edges[i + 1]) / 2
    z = (height[i] - exp_h[i]) / se_h[i]
    print(f"{mid:>10.4f}{height[i]:>12.5f}{exp_h[i]:>11.5f}{se_h[i]:>10.5f}{z:>8.3f}")

z = (height - exp_h) / se_h
print()
print(f"100 개 막대의 표준화잔차:  평균 {z.mean():+.4f}   표준편차 {z.std():.4f}")
print(f"|z| <= 2 인 비율 = {np.mean(np.abs(z) <= 2):.3f}   max|z| = {np.abs(z).max():.3f}")

n = 10000
mu-hat  = -0.018434   이론 E = +0.000000   이론 SE = 0.010000   z = -1.843
sig-hat = 0.987557   이론 E = 0.999925   이론 SE = 0.007071   z = -1.749
sig-hat^2 = 0.975268  이론 E = (n-1)/n = 0.999900   이론 SE = 0.014141
적합 봉우리 = 0.403969,  참 봉우리 = 0.398942,  차이 = +1.26%

칸 폭 Delta = 0.075418,  범위 = [-3.7401, 3.8017]
      칸 중앙       관측 높이       이론 E     이론 SE       z
   -2.9482     0.00663    0.00518   0.00262   0.554
   -1.4399     0.15911    0.14152   0.01363   1.291
   -0.3086     0.38320    0.38031   0.02213   0.131
    0.0685     0.39116    0.39791   0.02262  -0.299
    0.4456     0.37524    0.36117   0.02158   0.652
    1.5768     0.10873    0.11512   0.01230  -0.519
    3.0852     0.00133    0.00343   0.00213  -0.986

100 개 막대의 표준화잔차:  평균 -0.0878   표준편차 0.9426
|z| <= 2 인 비율 = 0.960   max|z| = 2.476

import pandas as pd

url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
df = pd.read_csv(url)
mean, std, n = df.x.mean(), df.x.std(), len(df.x)

n1 = len(df.x[(mean - std < df.x) & (df.x < mean + std)])
n2 = len(df.x[(mean - 2*std < df.x) & (df.x < mean + 2*std)])
n3 = len(df.x[(mean - 3*std < df.x) & (df.x < mean + 3*std)])

print(f"Within 1σ: {n1/n*100:.2f}%")   # ≈ 68%
print(f"Within 2σ: {n2/n*100:.2f}%")   # ≈ 95%
print(f"Within 3σ: {n3/n*100:.2f}%")   # ≈ 99.7%

Within 1σ: 72.69%
Within 2σ: 95.00%
Within 3σ: 98.66%

import numpy as np
import pandas as pd
from scipy import integrate, stats

# (1) 정규분포의 참값. cdf 와 quad 두 가지로 낸다.
print(f"{'k':>4}{'2*Phi(k)-1':>14}{'quad':>14}")
for k in (1, 2, 3):
    exact = 2 * stats.norm.cdf(k) - 1
    num = integrate.quad(stats.norm.pdf, -k, k)[0]
    print(f"{k:>4}{exact:>14.9f}{num:>14.9f}")
print(f"95% 를 정확히 주는 k = {stats.norm.ppf(0.975):.6f}")
print(f"k = 2 가 주는 확률    = {2 * stats.norm.cdf(2) - 1:.6f}")

# (2) 자료는 무엇이 다른가.
url = ('https://raw.githubusercontent.com/gedeck/'
       'practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv')
x = pd.read_csv(url).x
m, s = x.mean(), x.std()

print()
print(f"n = {len(x)},  평균 = {m:,.0f},  표준편차 = {s:,.0f},  중앙값 = {x.median():,.0f}")
print(f"왜도 = {stats.skew(x):.4f},  초과첨도 = {stats.kurtosis(x):.4f}")
print(f"자료 범위 = [{x.min():,}, {x.max():,}]")

print()
print(f"{'k':>4}{'자료':>10}{'정규':>10}{'차':>9}{'아래 밖':>10}{'위 밖':>9}")
for k in (1, 2, 3):
    emp = np.mean((m - k*s < x) & (x < m + k*s))
    thy = 2 * stats.norm.cdf(k) - 1
    print(f"{k:>4}{100*emp:>9.2f}%{100*thy:>9.2f}%{100*(emp-thy):>+8.2f}%"
          f"{100*np.mean(x <= m - k*s):>9.3f}%{100*np.mean(x >= m + k*s):>8.3f}%")

print()
print(f"평균 - 3*표준편차 = {m - 3*s:,.0f}   <- 소득이 음수일 수 없으므로 아래쪽은 비어 있다")
print(f"정규분포라면 위쪽 3 시그마 밖이 {100*stats.norm.sf(3):.3f}% 인데"
      f" 자료는 {100*np.mean(x >= m + 3*s):.3f}%")

   k    2*Phi(k)-1          quad
   1   0.682689492   0.682689492
   2   0.954499736   0.954499736
   3   0.997300204   0.997300204
95% 를 정확히 주는 k = 1.959964
k = 2 가 주는 확률    = 0.954500

n = 50000,  평균 = 68,761,  표준편차 = 32,872,  중앙값 = 62,000
왜도 = 1.0488,  초과첨도 = 1.0808
자료 범위 = [4,000, 199,000]

   k        자료        정규        차      아래 밖      위 밖
   1    72.69%    68.27%   +4.42%   12.746%  14.564%
   2    95.00%    95.45%   -0.45%    0.000%   5.002%
   3    98.66%    99.73%   -1.07%    0.000%   1.338%

평균 - 3*표준편차 = -29,856   <- 소득이 음수일 수 없으므로 아래쪽은 비어 있다
정규분포라면 위쪽 3 시그마 밖이 0.135% 인데 자료는 1.338%

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def shade_area(z_bounds, side='left', ax=None):
    """표준정규곡선 아래의 한 구간을 칠한다."""
    x = np.linspace(-4, 4, 200)
    ax.plot(x, stats.norm().pdf(x), color='k', alpha=0.9)

    if side == 'left':
        x_shade = np.linspace(-4, z_bounds, 200)
    elif side == 'right':
        x_shade = np.linspace(z_bounds, 4, 200)
    else:  # center
        x_shade = np.linspace(z_bounds[0], z_bounds[1], 200)

    ax.fill_between(x_shade, stats.norm().pdf(x_shade), alpha=0.2, color='k')
    ax.spines[['left', 'right', 'top']].set_visible(False)
    ax.spines['bottom'].set_position('zero')
    ax.set_yticks([])

# 왼쪽 넓이
z = -1.2
print(f"P(Z ≤ {z}) = {stats.norm().cdf(z):.4f}")

# 오른쪽 넓이
z = 1.2
print(f"P(Z ≥ {z}) = {stats.norm().sf(z):.4f}")

# 가운데 넓이
z1, z2 = -2.1, 1.2
print(f"P({z1} ≤ Z ≤ {z2}) = {stats.norm().cdf(z2) - stats.norm().cdf(z1):.4f}")

P(Z ≤ -1.2) = 0.1151
P(Z ≥ 1.2) = 0.1151
P(-2.1 ≤ Z ≤ 1.2) = 0.8671

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 제목에 한글을 쓴다
plt.rcParams["axes.unicode_minus"] = False

fig, axes = plt.subplots(1, 3, figsize=(12, 2.6))

shade_area(-1.2, side='left', ax=axes[0])
axes[0].set_title("왼쪽 넓이  $P(Z \\leq -1.2) = 0.1151$", fontsize=10)

shade_area(1.2, side='right', ax=axes[1])
axes[1].set_title("오른쪽 넓이  $P(Z \\geq 1.2) = 0.1151$", fontsize=10)

shade_area((-2.1, 1.2), side='center', ax=axes[2])
axes[2].set_title("가운데 넓이  $P(-2.1 \\leq Z \\leq 1.2) = 0.8671$", fontsize=10)

for ax in axes:
    ax.set_xticks([-4, -2, 0, 2, 4])
    ax.set_xlabel("$z$")

plt.tight_layout()
plt.show()

import numpy as np
from scipy import integrate, stats

cdf, sf, pdf = stats.norm.cdf, stats.norm.sf, stats.norm.pdf

# 세 넓이를 닫힌 꼴(Phi)과 수치적분(quad)으로 각각 구해 맞춰 본다.
rows = [
    ("P(Z <= -1.2)", cdf(-1.2), integrate.quad(pdf, -np.inf, -1.2)),
    ("P(Z >=  1.2)", sf(1.2), integrate.quad(pdf, 1.2, np.inf)),
    ("P(-2.1<=Z<=1.2)", cdf(1.2) - cdf(-2.1), integrate.quad(pdf, -2.1, 1.2)),
]
print(f"{'넓이':>17}{'Phi 로':>14}{'quad 로':>14}{'quad 오차한계':>16}")
for name, closed, (num, err) in rows:
    print(f"{name:>17}{closed:>14.9f}{num:>14.9f}{err:>16.2e}")

# 앞의 두 값이 같은 것은 우연이 아니다. 비트까지 같은지 본다.
print()
print(f"cdf(-1.2) == sf(1.2) ?  {cdf(-1.2) == sf(1.2)}")

# 가운데 넓이는 두 꼬리를 뺀 것이기도 하다.
print(f"1 - P(Z<-2.1) - P(Z>1.2) = {1 - cdf(-2.1) - sf(1.2):.9f}")

# 그림이 그리는 구간 [-4, 4] 밖의 확률. 칠한 넓이에 들어가지 않는다.
print(f"[-4, 4] 밖의 확률 = {2 * sf(4):.3e}")
print(f"왼쪽 칸이 실제로 칠한 넓이 = {cdf(-1.2) - cdf(-4):.9f}"
      f"  (참값 {cdf(-1.2):.9f} 보다 {cdf(-4):.2e} 작다)")

               넓이         Phi 로        quad 로       quad 오차한계
     P(Z <= -1.2)   0.115069670   0.115069670        1.31e-10
     P(Z >=  1.2)   0.115069670   0.115069670        1.31e-10
  P(-2.1<=Z<=1.2)   0.867065909   0.867065909        1.75e-14

cdf(-1.2) == sf(1.2) ?  True
1 - P(Z<-2.1) - P(Z>1.2) = 0.867065909
[-4, 4] 밖의 확률 = 6.334e-05
왼쪽 칸이 실제로 칠한 넓이 = 0.115037999  (참값 0.115069670 보다 3.17e-05 작다)

import numpy as np
import scipy.stats as stats

np.random.seed(42)
samples = stats.norm.rvs(size=10)
print(samples)  # 시드가 42면 언제나 같은 값이 나온다

[ 0.49671415 -0.1382643   0.64768854  1.52302986 -0.23415337 -0.23413696
  1.57921282  0.76743473 -0.46947439  0.54256004]

import numpy as np
import scipy.stats as stats

# (1) 같은 씨앗은 같은 배열을, 다른 씨앗은 다른 배열을 준다.
np.random.seed(42); a = stats.norm.rvs(size=10)
np.random.seed(42); b = stats.norm.rvs(size=10)
np.random.seed(43); c = stats.norm.rvs(size=10)
print(f"seed(42) 를 두 번:  비트까지 같은가 = {np.array_equal(a, b)}")
print(f"seed(42) 대 seed(43): 비트까지 같은가 = {np.array_equal(a, c)}")

# 그러나 표본오차는 그대로 남아 있다.
print()
print(f"seed(42) 표본의 평균 = {a.mean():+.6f}   (참 평균은 0)")
print(f"seed(43) 표본의 평균 = {c.mean():+.6f}")
print(f"이론 SE = 1/sqrt(10) = {1/np.sqrt(10):.6f}"
      f"   seed(42) 의 z = {a.mean()*np.sqrt(10):+.3f}")

# 씨앗을 1000 개 바꿔 보면 표본평균이 이론 SE 만큼 흩어진다.
means = []
for s in range(1000):
    np.random.seed(s)
    means.append(stats.norm.rvs(size=10).mean())
means = np.array(means)
print(f"씨앗 1000 개의 표본평균:  평균 {means.mean():+.4f}"
      f"   표준편차 {means.std(ddof=1):.4f}   (이론 {1/np.sqrt(10):.4f})")
print(f"|평균| <= 2*SE 인 비율 = {np.mean(np.abs(means) <= 2/np.sqrt(10)):.3f}")

# (2) 전역 시드와 새 생성기는 같은 씨앗이라도 다른 수열을 준다.
rng = np.random.default_rng(42)
d = stats.norm.rvs(size=10, random_state=rng)
print()
print(f"seed(42)        앞 세 개 = {np.round(a[:3], 8)}")
print(f"default_rng(42) 앞 세 개 = {np.round(d[:3], 8)}")
print(f"같은가 = {np.array_equal(a, d)}")

seed(42) 를 두 번:  비트까지 같은가 = True
seed(42) 대 seed(43): 비트까지 같은가 = False

seed(42) 표본의 평균 = +0.448061   (참 평균은 0)
seed(43) 표본의 평균 = +0.221260
이론 SE = 1/sqrt(10) = 0.316228   seed(42) 의 z = +1.417
씨앗 1000 개의 표본평균:  평균 -0.0071   표준편차 0.3105   (이론 0.3162)
|평균| <= 2*SE 인 비율 = 0.962

seed(42)        앞 세 개 = [ 0.49671415 -0.1382643   0.64768854]
default_rng(42) 앞 세 개 = [ 0.30471708 -1.03998411  0.7504512 ]
같은가 = False

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

# 로그정규분포: log(X)가 N(mu, sigma^2)을 따르는 분포다.
# mu와 sigma는 **로그를 취한 뒤의** 평균과 표준편차이지 X 자체의 것이 아니다.
mu = 0
sigmas = [0.5, 1.0, 1.5, 2.0]
x = np.linspace(0.001, 8, 500)     # X > 0 이므로 0에서 시작한다

fig, ax = plt.subplots(figsize=(12, 4))
for sigma in sigmas:
    #   s     = 로그 척도의 표준편차 sigma
    #   scale = exp(mu)  <- loc가 아니라 scale에 넣는다
    rv = stats.lognorm(s=sigma, scale=np.exp(mu))
    # sigma가 커질수록 봉우리가 0쪽으로 밀리고 오른쪽 꼬리가 길어진다.
    # mu=0 이라 중앙값은 네 곡선 모두 exp(0)=1 로 같다는 점을 확인하라.
    ax.plot(x, rv.pdf(x), label=rf'$\sigma={sigma}$')
ax.set_xlabel('x')
ax.set_ylabel('f(x)')
ax.set_title(r'Log-Normal Distribution — PDF ($\mu=0$, varying $\sigma$)')
ax.legend()
ax.set_ylim(bottom=-0.02)
plt.tight_layout()
plt.show()

import numpy as np
from scipy import integrate, optimize, stats

mu = 0
sigmas = [0.5, 1.0, 1.5, 2.0]

# (1) 세 대푯값의 닫힌 꼴과 수치값.
print(f"{'sigma':>6}{'최빈값':>11}{'중앙값':>10}{'평균':>11}"
      f"{'격자 최빈':>12}{'quad 평균':>12}")
for s in sigmas:
    rv = stats.lognorm(s=s, scale=np.exp(mu))
    g = np.linspace(1e-6, 20, 2_000_001)
    mode_num = g[rv.pdf(g).argmax()]
    mean_num = integrate.quad(lambda y: y * rv.pdf(y), 0, np.inf, limit=200)[0]
    print(f"{s:>6}{np.exp(mu - s*s):>11.6f}{np.exp(mu):>10.6f}{np.exp(mu + s*s/2):>11.6f}"
          f"{mode_num:>12.6f}{mean_num:>12.6f}")

# 분산도 맞춰 본다.
print()
for s in sigmas:
    rv = stats.lognorm(s=s, scale=np.exp(mu))
    m = integrate.quad(lambda y: y * rv.pdf(y), 0, np.inf, limit=200)[0]
    m2 = integrate.quad(lambda y: y*y * rv.pdf(y), 0, np.inf, limit=200)[0]
    closed = np.exp(2*mu + s*s) * (np.exp(s*s) - 1)
    print(f"sigma={s}: Var quad = {m2 - m*m:>12.6f}   닫힌 꼴 = {closed:>12.6f}"
          f"   CV = {np.sqrt(np.exp(s*s) - 1):.4f}")

# (2) 봉우리 높이. 닫힌 꼴 exp(s^2/2 - mu)/(s*sqrt(2pi)) 와 실제 pdf 를 맞춘다.
print()
print(f"{'sigma':>6}{'봉우리 높이':>14}{'닫힌 꼴':>12}{'그림에 보이는 최대':>20}{'P(Y>8)':>10}")
x = np.linspace(0.001, 8, 500)          # 쪽의 그림이 쓰는 격자
for s in sigmas:
    rv = stats.lognorm(s=s, scale=np.exp(mu))
    peak = rv.pdf(np.exp(mu - s*s))
    closed = np.exp(s*s/2 - mu) / (s * np.sqrt(2*np.pi))
    print(f"{s:>6}{peak:>14.6f}{closed:>12.6f}{rv.pdf(x).max():>20.6f}{rv.sf(8):>10.4f}")

r = optimize.minimize_scalar(lambda s: np.exp(s*s/2) / (s*np.sqrt(2*np.pi)),
                             bounds=(0.1, 3), method='bounded')
print(f"봉우리가 가장 낮아지는 sigma = {r.x:.6f}   그 높이 = {r.fun:.6f}"
      f"   sqrt(e/(2*pi)) = {np.sqrt(np.e/(2*np.pi)):.6f}")

 sigma        최빈값       중앙값         평균       격자 최빈     quad 평균
   0.5   0.778801  1.000000   1.133148    0.778801    1.133148
   1.0   0.367879  1.000000   1.648721    0.367881    1.648721
   1.5   0.105399  1.000000   3.080217    0.105401    3.080217
   2.0   0.018316  1.000000   7.389056    0.018311    7.389056

sigma=0.5: Var quad =     0.364696   닫힌 꼴 =     0.364696   CV = 0.5329
sigma=1.0: Var quad =     4.670774   닫힌 꼴 =     4.670774   CV = 1.3108
sigma=1.5: Var quad =    80.529395   닫힌 꼴 =    80.529395   CV = 2.9134
sigma=2.0: Var quad =  2926.359837   닫힌 꼴 =  2926.359837   CV = 7.3211

 sigma        봉우리 높이        닫힌 꼴          그림에 보이는 최대    P(Y>8)
   0.5      0.904122    0.904122            0.903948    0.0000
   1.0      0.657745    0.657745            0.657737    0.0188
   1.5      0.819219    0.819219            0.818289    0.0828
   2.0      1.473903    1.473903            1.472928    0.1492
봉우리가 가장 낮아지는 sigma = 1.000000   그 높이 = 0.657745   sqrt(e/(2*pi)) = 0.657745