확률질량함수, 확률밀도함수, 누적분포함수¶

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

# 왼쪽에 PDF, 오른쪽에 CDF, 가운데에는 둘의 관계를 나타내는 화살표를 놓는다.
fig, (ax_pdf, ax_arrow, ax_cdf) = plt.subplots(1, 3, figsize=(12, 3))

x = np.linspace(-3, 3, 100)

# 왼쪽: 확률밀도함수. 각 점에서 확률이 얼마나 빽빽한지를 나타낸다.
ax_pdf.set_title("PDF", fontsize=16)
ax_pdf.plot(x, stats.norm().pdf(x))

# 가운데: 두 함수를 잇는 연산.
#   PDF -> CDF 는 적분 (왼쪽에서 여기까지의 넓이를 쌓는다)
#   CDF -> PDF 는 미분 (누적이 늘어나는 속도가 곧 밀도다)
ax_arrow.arrow(0.1, 0.6, 0.8, 0, width=0.05, length_includes_head=True)
ax_arrow.arrow(0.9, 0.4, -0.8, 0, width=0.05, length_includes_head=True)
ax_arrow.annotate("Integrate", (0.38, 0.75), fontsize=14)
ax_arrow.annotate("Differentiate", (0.30, 0.2), fontsize=14)
for spine in ax_arrow.spines.values():
    spine.set_visible(False)
ax_arrow.set_xticks([])
ax_arrow.set_yticks([])

# CDF
ax_cdf.set_title("CDF", fontsize=16)
ax_cdf.plot(x, stats.norm().cdf(x))

for ax in (ax_pdf, ax_cdf):
    ax.set_ylim(0, 1)
plt.tight_layout()
plt.show()

# 그림에서 읽을 수치. F'(x) = f(x) 이므로 CDF 의 기울기가 곧 PDF 의 높이다.
nd = stats.norm()
print(f"{'x':>5}{'f(x) = F 의 기울기':>18}{'F(x)':>10}")
for v in (0.0, 1.0, 2.0, 3.0):
    print(f"{v:>5.1f}{nd.pdf(v):>18.6f}{nd.cdf(v):>10.6f}")

# 차분으로 기울기를 재어 밀도와 맞춰 본다.
h = 1e-5
for v in (0.0, 1.0, 3.0):
    slope = (nd.cdf(v + h) - nd.cdf(v - h)) / (2 * h)
    print(f"  x={v:.1f}:  (F(x+h)-F(x-h))/2h = {slope:.8f},  f(x) = {nd.pdf(v):.8f}")

print(f"\n[-1, 1] 에서 밀도 아래 넓이 = F(1) - F(-1) = {nd.cdf(1) - nd.cdf(-1):.6f}")
print(f"  같은 구간에서 CDF 가 오른 높이도 {nd.cdf(1) - nd.cdf(-1):.6f}")

# 세로축을 [0, 1] 로 잘라도 되는가. 분산이 작아지면 밀도가 1 을 넘는다.
print(f"\n밀도의 최댓값 1/(sigma*sqrt(2pi))")
for s in (1.0, 0.5, 0.2):
    print(f"  sigma = {s}:  {1 / (s * np.sqrt(2 * np.pi)):.6f}")

    x    f(x) = F 의 기울기      F(x)
  0.0          0.398942  0.500000
  1.0          0.241971  0.841345
  2.0          0.053991  0.977250
  3.0          0.004432  0.998650
  x=0.0:  (F(x+h)-F(x-h))/2h = 0.39894228,  f(x) = 0.39894228
  x=1.0:  (F(x+h)-F(x-h))/2h = 0.24197072,  f(x) = 0.24197072
  x=3.0:  (F(x+h)-F(x-h))/2h = 0.00443185,  f(x) = 0.00443185

[-1, 1] 에서 밀도 아래 넓이 = F(1) - F(-1) = 0.682689
  같은 구간에서 CDF 가 오른 높이도 0.682689

밀도의 최댓값 1/(sigma*sqrt(2pi))
  sigma = 1.0:  0.398942
  sigma = 0.5:  0.797885
  sigma = 0.2:  1.994711

from scipy import stats

mean, std_dev = 50, 10

# P(40 ≤ X ≤ 60) for X ~ N(50, 10²)
prob = stats.norm(mean, std_dev).cdf(60) - stats.norm(mean, std_dev).cdf(40)
print(f"P(40 ≤ X ≤ 60) = {prob * 100:.2f}%")

# P(X ≤ 55)
prob_55 = stats.norm(mean, std_dev).cdf(55)
print(f"P(X ≤ 55) = {prob_55 * 100:.2f}%")

P(40 ≤ X ≤ 60) = 68.27%
P(X ≤ 55) = 69.15%

import scipy.stats as stats

z_95 = stats.norm(0, 1).ppf(0.95)
print(f"95th percentile of N(0,1): {z_95:.4f}")

z_975 = stats.norm(0, 1).ppf(0.975)
print(f"97.5th percentile of N(0,1): {z_975:.4f}")

# 두 수가 각각 어떤 꼬리를 자르는지 확인한다.
nd = stats.norm()
print(f"\n{'z':>10}{'P(Z > z)':>12}{'P(|Z| <= z)':>14}")
for z in (z_95, z_975, 1.96):
    print(f"{z:>10.6f}{nd.sf(z):>12.6f}{2 * nd.cdf(z) - 1:>14.6f}")

# ppf 에 닫힌 꼴이 없으므로 수치로 푼다. 이분법으로 Phi(z) = 0.975 를 푼다.
lo, hi = 0.0, 4.0
print(f"\n이분법으로 Phi(z) = 0.975 를 푼다")
for step in range(1, 41):
    mid = (lo + hi) / 2
    if nd.cdf(mid) < 0.975:
        lo = mid
    else:
        hi = mid
    if step in (5, 10, 20, 40):
        print(f"  {step:>2}단계:  z = {mid:.10f}")
print(f"  scipy 의 ppf = {z_975:.10f}")

95th percentile of N(0,1): 1.6449
97.5th percentile of N(0,1): 1.9600

         z    P(Z > z)   P(|Z| <= z)
  1.644854    0.050000      0.900000
  1.959964    0.025000      0.950000
  1.960000    0.024998      0.950004

이분법으로 Phi(z) = 0.975 를 푼다
   5단계:  z = 1.8750000000
  10단계:  z = 1.9570312500
  20단계:  z = 1.9599647522
  40단계:  z = 1.9599639845
  scipy 의 ppf = 1.9599639845

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

fig, ax = plt.subplots(figsize=(12, 3))
ax.set_xlim(-3, 3)
ax.set_ylim(-0.2, 1.1)

# CDF 곡선
x = np.linspace(-3, 3, 100)
ax.plot(x, stats.norm().cdf(x), label='CDF')

# CDF와 PPF는 서로 역함수다.
#   CDF: 값 z 를 넣으면 누적확률 u 가 나온다      (가로 -> 세로)
#   PPF: 누적확률 u 를 넣으면 값 z 가 나온다      (세로 -> 가로)
# 아래 두 점이 같은 (z, u) 쌍을 축마다 표시한 것이다.
u = 0.975
z = stats.norm().ppf(u)      # 표준정규분포의 97.5 백분위수. 그 유명한 1.96이다.

ax.plot(0, u, 'or', markersize=8)
ax.plot(z, 0, 'or', markersize=8)
ax.annotate(f"U = {u}", (-1.2, u + 0.02), fontsize=14)
ax.annotate(f"Z = {z:.3f}", (z - 0.3, -0.12), fontsize=14)
ax.annotate("PPF →", (0.3, u + 0.03), fontsize=14)
ax.annotate("↓ CDF", (z + 0.1, 0.5), fontsize=14)

ax.spines[['right', 'top']].set_visible(False)
ax.spines['left'].set_position('zero')
ax.spines['bottom'].set_position('zero')
ax.legend(fontsize=14)
plt.show()

# 왕복이 제자리로 돌아오는지 확인한다.
nd = stats.norm()
print(f"F(F^-1(u)) = {nd.cdf(nd.ppf(u)):.10f}   (u = {u})")
print(f"F^-1(F(z)) = {nd.ppf(nd.cdf(z)):.10f}   (z = {z:.6f})")

# 분위수함수의 기울기는 밀도의 역수다.  (F^-1)'(p) = 1 / f(F^-1(p))
print(f"\n{'p':>8}{'z = F^-1(p)':>14}{'f(z)':>12}{'1/f(z)':>12}")
for p in (0.5, 0.975, 0.99, 0.999):
    zp = nd.ppf(p)
    print(f"{p:>8.3f}{zp:>14.6f}{nd.pdf(zp):>12.6f}{1 / nd.pdf(zp):>12.4f}")

# 같은 폭 0.015 를 가운데와 꼬리에서 각각 움직여 본다.
print(f"\n같은 Delta u = 0.015 가 z 를 얼마나 움직이는가")
for p in (0.500, 0.975):
    dz = nd.ppf(p + 0.015) - nd.ppf(p)
    print(f"  p: {p:.3f} -> {p + 0.015:.3f}   Delta z = {dz:.6f}")

F(F^-1(u)) = 0.9750000000   (u = 0.975)
F^-1(F(z)) = 1.9599639845   (z = 1.959964)

       p   z = F^-1(p)        f(z)      1/f(z)
   0.500      0.000000    0.398942      2.5066
   0.975      1.959964    0.058445     17.1101
   0.990      2.326348    0.026652     37.5204
   0.999      3.090232    0.003367    296.9924

같은 Delta u = 0.015 가 z 를 얼마나 움직이는가
  p: 0.500 -> 0.515   Delta z = 0.037608
  p: 0.975 -> 0.990   Delta z = 0.366384

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

np.random.seed(0)

# 역변환 표집: 균등난수만 있으면 어떤 분포든 만들어 낼 수 있다.
# 1단계 — 0과 1 사이 균등난수를 뽑는다. 이것이 "누적확률"에 해당한다.
u = stats.uniform().rvs(10_000)

# 2단계 — 그 누적확률에 대응하는 값을 PPF로 되찾는다.
# U ~ Uniform(0,1) 이면 F^{-1}(U) 는 정확히 F를 분포함수로 갖는다.
z = stats.norm().ppf(u)

plt.figure(figsize=(12, 3))
plt.hist(z, bins=100, density=True, alpha=0.7, label='Inverse Transform Samples')
x = np.linspace(-4, 4, 200)
# 만들어 낸 표본의 히스토그램이 참 정규 밀도와 겹치는지 확인한다
plt.plot(x, stats.norm().pdf(x), 'r--', lw=2, label='N(0,1) PDF')
plt.legend()
plt.show()

# 만들어 낸 표본이 정말 N(0,1) 인지 수로 확인한다.
n = len(z)
nd = stats.norm()
print(f"표본 {n}개")
print(f"  평균     {z.mean():+.4f}  (이론 0,  SE {1 / np.sqrt(n):.4f},"
      f"  z = {z.mean() * np.sqrt(n):+.2f})")
print(f"  표준편차 {z.std(ddof=1):.4f}  (이론 1,  SE {1 / np.sqrt(2 * n):.4f},"
      f"  z = {(z.std(ddof=1) - 1) * np.sqrt(2 * n):+.2f})")

print(f"\n{'p':>7}{'표본 분위수':>12}{'이론':>10}{'SE':>9}{'z':>8}")
for p in (0.05, 0.25, 0.5, 0.75, 0.95):
    q_hat, q = np.quantile(z, p), nd.ppf(p)
    se = np.sqrt(p * (1 - p) / n) / nd.pdf(q)
    print(f"{p:>7.2f}{q_hat:>12.4f}{q:>10.4f}{se:>9.4f}{(q_hat - q) / se:>8.2f}")

ks = stats.kstest(z, "norm")
print(f"\n콜모고로프-스미르노프  D = {ks.statistic:.5f},  p = {ks.pvalue:.4f}")

# F 에 계단(점질량)이 있어도 같은 방법이 통한다.
# P(X=0)=0.5, P(X=1)=0.3, P(X=2)=0.2 인 분포를 균등난수로 만들어 본다.
probs = np.array([0.5, 0.3, 0.2])
cuts = np.cumsum(probs)                  # F 의 계단 높이 0.5, 0.8, 1.0
x_disc = np.searchsorted(cuts, u)        # F^{-1}(u) = inf{x : F(x) >= u}
print(f"\n점질량이 있는 F 에도 통한다")
for k in range(3):
    print(f"  P(X={k}): 모의 {np.mean(x_disc == k):.4f}   참값 {probs[k]:.4f}"
          f"   SE {np.sqrt(probs[k] * (1 - probs[k]) / n):.4f}")

표본 10000개
  평균     -0.0145  (이론 0,  SE 0.0100,  z = -1.45)
  표준편차 1.0078  (이론 1,  SE 0.0071,  z = +1.11)

      p      표본 분위수        이론       SE       z
   0.05     -1.6715   -1.6449   0.0211   -1.26
   0.25     -0.6893   -0.6745   0.0136   -1.09
   0.50     -0.0163    0.0000   0.0125   -1.30
   0.75      0.6691    0.6745   0.0136   -0.39
   0.95      1.6373    1.6449   0.0211   -0.36

콜모고로프-스미르노프  D = 0.00814,  p = 0.5195

점질량이 있는 F 에도 통한다
  P(X=0): 모의 0.5064   참값 0.5000   SE 0.0050
  P(X=1): 모의 0.2964   참값 0.3000   SE 0.0046
  P(X=2): 모의 0.1972   참값 0.2000   SE 0.0040

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

np.random.seed(42)
data = stats.norm.rvs(size=200)      # 표준정규에서 200개

fig, ax = plt.subplots(figsize=(12, 4))

# 경험적 PDF: 히스토그램이 밀도함수의 표본 버전이다
counts, bin_edges, _ = ax.hist(data, bins=20, density=True, alpha=0.6, label="Empirical PDF")

# 경험적 CDF: 히스토그램의 도수를 왼쪽부터 누적하면 된다.
# "PDF를 적분하면 CDF"라는 관계를 이산 버전으로 실행한 것이다.
empirical_cdf = np.cumsum(counts) / np.sum(counts)
# 계단으로 그린다. 경험적 분포함수는 본래 계단함수이기 때문이다.
ax.step(bin_edges[1:], empirical_cdf, where='mid', label="Empirical CDF", lw=2)

# 이론적 CDF를 겹쳐 표본이 모집단을 얼마나 잘 따라가는지 본다
ax.plot(bin_edges, stats.norm.cdf(bin_edges), 'r', lw=2, label="Theoretical CDF")

ax.legend()
ax.spines[['right', 'top']].set_visible(False)
plt.show()

width = bin_edges[1] - bin_edges[0]
print(f"표본 {len(data)}개,  구간 폭 = {width:.4f},  범위 [{data.min():.3f}, {data.max():.3f}]")
print(f"히스토그램 최고 높이 = {counts.max():.4f}  (이론 최댓값 {1 / np.sqrt(2 * np.pi):.4f})")

# 계단의 값이 각 구간 오른쪽 끝에서의 경험적 분포함수와 같은지 확인한다.
print(f"\n{'오른쪽 끝':>10}{'누적':>9}{'참 ECDF':>10}{'이론 F':>10}")
for i in (4, 9, 14, 19):
    edge = bin_edges[i + 1]
    print(f"{edge:>10.4f}{empirical_cdf[i]:>9.4f}"
          f"{np.mean(data <= edge):>10.4f}{stats.norm.cdf(edge):>10.4f}")

# where='mid' 는 계단의 모서리를 오른쪽 끝이 아니라 그 반 칸 왼쪽에 놓는다.
print(f"\nwhere='mid' 가 계단을 옮기는 거리 = 반 칸 = {width / 2:.4f}")

# 경험적 분포함수와 이론 분포함수의 최대 거리(콜모고로프-스미르노프 통계량).
ks = stats.kstest(data, "norm")
print(f"\n경험적 F 와 이론 F 의 최대 거리  D = {ks.statistic:.5f},  p = {ks.pvalue:.4f}")
print(f"  n = {len(data)} 에서 5% 임계값 ~ 1.36/sqrt(n) = {1.36 / np.sqrt(len(data)):.5f}")

표본 200개,  구간 폭 = 0.2670,  범위 [-2.620, 2.720]
히스토그램 최고 높이 = 0.5618  (이론 최댓값 0.3989)

     오른쪽 끝       누적    참 ECDF      이론 F
   -1.2848   0.0900    0.0900    0.0994
    0.0502   0.5100    0.5100    0.5200
    1.3852   0.9200    0.9200    0.9170
    2.7202   1.0000    1.0000    0.9967

where='mid' 가 계단을 옮기는 거리 = 반 칸 = 0.1335

경험적 F 와 이론 F 의 최대 거리  D = 0.06500,  p = 0.3515
  n = 200 에서 5% 임계값 ~ 1.36/sqrt(n) = 0.09617

import numpy as np

rng = np.random.default_rng(0)
pk = np.array([0.1, 0.3, 0.4, 0.2])
F = np.cumsum(pk)

u = rng.random(500_000)
x = np.searchsorted(F, u, side="left")        # 일반화 역함수

print("일반화 역함수로 이산분포 표집")
for v in range(4):
    print(f"  P(X={v}): 모의 {np.mean(x == v):.5f}   참값 {pk[v]:.5f}")

일반화 역함수로 이산분포 표집
  P(X=0): 모의 0.10026   참값 0.10000
  P(X=1): 모의 0.29910   참값 0.30000
  P(X=2): 모의 0.40091   참값 0.40000
  P(X=3): 모의 0.19974   참값 0.20000

import numpy as np

lam = 0.7                                      # 지수분포
print(f"{'t':>6}{'S(t)':>10}{'f(t)':>10}{'h(t)':>10}{'H(t)':>10}{'exp(-H)':>10}")
for t in np.linspace(0.2, 3.0, 5):
    S = np.exp(-lam * t)
    f = lam * np.exp(-lam * t)
    print(f"{t:>6.2f}{S:>10.5f}{f:>10.5f}{f / S:>10.5f}{lam * t:>10.5f}"
          f"{np.exp(-lam * t):>10.5f}")

     t      S(t)      f(t)      h(t)      H(t)   exp(-H)
  0.20   0.86936   0.60855   0.70000   0.14000   0.86936
  0.90   0.53259   0.37281   0.70000   0.63000   0.53259
  1.60   0.32628   0.22840   0.70000   1.12000   0.32628
  2.30   0.19989   0.13992   0.70000   1.61000   0.19989
  3.00   0.12246   0.08572   0.70000   2.10000   0.12246

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n = 400_000
z = rng.multivariate_normal([0, 0], [[1, 0.7], [0.7, 1]], n)

u1, u2 = stats.norm.cdf(z[:, 0]), stats.norm.cdf(z[:, 1])   # 코퓰라로 이동
x1 = stats.expon.ppf(u1)                                     # 주변분포만 갈아 끼운다
x2 = stats.lognorm.ppf(u2, 1)

print("같은 코퓰라, 다른 주변분포")
print(f"{'':>22}{'피어슨':>10}{'스피어만':>11}")
print(f"{'(정규, 정규)':>22}{np.corrcoef(z.T)[0,1]:>10.4f}"
      f"{stats.spearmanr(z)[0]:>11.4f}")
print(f"{'(지수, 로그정규)':>22}{np.corrcoef(x1, x2)[0,1]:>10.4f}"
      f"{stats.spearmanr(x1, x2)[0]:>11.4f}")

같은 코퓰라, 다른 주변분포
                             피어슨       스피어만
              (정규, 정규)    0.7000     0.6828
            (지수, 로그정규)    0.6024     0.6828

import warnings
import numpy as np
from scipy import integrate

warnings.simplefilter("ignore")               # 발산 경고를 잠시 끈다

candidates = [("1/x        (x>=1)", lambda x: 1 / x, 1),
              ("1/x^2      (x>=1)", lambda x: 1 / x ** 2, 1),
              ("1/(x ln^2 x) (x>=2)", lambda x: 1 / (x * np.log(x) ** 2), 2)]

print("f 의 부분적분 — 분포가 되는가")
print(f"{'f(x)':>20}{'[lo, 10]':>12}{'[lo, 10^3]':>13}{'[lo, 10^6]':>13}")
for name, g, lo in candidates:
    vals = [integrate.quad(g, lo, T, limit=400)[0] for T in (10, 1e3, 1e6)]
    print(f"{name:>20}" + "".join(f"{v:>13.5f}" for v in vals))

print("\nx f(x) 의 부분적분 — 평균이 존재하는가")
print(f"{'f(x)':>20}{'[lo, 10]':>12}{'[lo, 10^3]':>13}{'[lo, 10^6]':>13}")
for name, g, lo in [candidates[1], candidates[2], ("3/x^4      (x>=1)", lambda x: 3 / x ** 4, 1)]:
    vals = [integrate.quad(lambda x: x * g(x), lo, T, limit=400)[0] for T in (10, 1e3, 1e6)]
    print(f"{name:>20}" + "".join(f"{v:>13.5f}" for v in vals))

f 의 부분적분 — 분포가 되는가
                f(x)    [lo, 10]   [lo, 10^3]   [lo, 10^6]
   1/x        (x>=1)      2.30259      6.90776     13.81551
   1/x^2      (x>=1)      0.90000      0.99900      1.00000
 1/(x ln^2 x) (x>=2)      1.00840      1.29793      1.37031

x f(x) 의 부분적분 — 평균이 존재하는가
                f(x)    [lo, 10]   [lo, 10^3]   [lo, 10^6]
   1/x^2      (x>=1)      2.30259      6.90776     13.81551
 1/(x ln^2 x) (x>=2)      3.66288     34.68506   6246.97574
   3/x^4      (x>=1)      1.48500      1.50000     -0.00000