Fonction génératrice et preuve du TCL

import numpy as np

# Y ~ N(0, 1) and X = exp(Y) is log-normal. Its "sisters" reweight the density
# of ln X by (1 + a sin(2 pi y)), |a| <= 1.
y = np.linspace(-12, 16, 400_001)
dy = y[1] - y[0]
phi = np.exp(-y**2 / 2) / np.sqrt(2 * np.pi)

for a in (0.0, 0.5, 1.0):
    w = phi * (1 + a * np.sin(2 * np.pi * y))       # density of ln X for this sister
    assert w.min() >= 0                             # still a density
    moments = [np.sum(np.exp(k * y) * w) * dy for k in range(5)]   # E[X^k]
    cdf2 = np.sum(w[y <= np.log(2)]) * dy           # P(X <= 2)
    m = "  ".join(f"{v:.3f}" for v in moments)
    print(f"a = {a}:  E[X^k], k = 0..4: {m}   P(X <= 2) = {cdf2:.4f}")

m = "  ".join(f"{np.exp(k * k / 2):.3f}" for k in range(5))
print(f"theory e^(k^2/2):          {m}")

a = 0.0:  E[X^k], k = 0..4: 1.000  1.649  7.389  90.017  2980.958   P(X <= 2) = 0.7559
a = 0.5:  E[X^k], k = 0..4: 1.000  1.649  7.389  90.017  2980.958   P(X <= 2) = 0.7675
a = 1.0:  E[X^k], k = 0..4: 1.000  1.649  7.389  90.017  2980.958   P(X <= 2) = 0.7791
theory e^(k^2/2):          1.000  1.649  7.389  90.017  2980.958

Mêmes moments jusqu'à l'ordre 4, et en fait à tous les ordres, mais trois fonctions de répartition différentes : P(X ≤ 2) vaut 0,756, 0,768 ou 0,779.

import numpy as np

lam, h = 3.0, 1e-3
M = lambda t: np.exp(lam * (np.exp(t) - 1))         # Poisson MGF

d1 = (M(h) - M(-h)) / (2 * h)                       # slope at 0
d2 = (M(h) - 2 * M(0) + M(-h)) / h**2               # curvature at 0
d3 = (M(2 * h) - 2 * M(h) + 2 * M(-h) - M(-2 * h)) / (2 * h**3)
print(f"M'(0)   = {d1:6.3f}   theory lam                 = {lam}")
print(f"M''(0)  = {d2:6.3f}   theory lam^2 + lam         = {lam**2 + lam}")
print(f"M'''(0) = {d3:6.3f}   theory lam^3 + 3lam^2 + lam = {lam**3 + 3 * lam**2 + lam}")
print(f"Var = M''(0) - M'(0)^2 = {d2 - d1**2:.3f}")

x = np.random.default_rng(0).poisson(lam, 1_000_000)
print("sample E[X], E[X^2], E[X^3]:", [round(float(np.mean(x**k)), 2) for k in (1, 2, 3)])
mu, sd = x.mean(), x.std()
print(f"skewness {np.mean(((x - mu) / sd)**3):.3f}   theory 1/sqrt(lam) = {lam**-0.5:.3f}")

M'(0)   =  3.000   theory lam                 = 3.0
M''(0)  = 12.000   theory lam^2 + lam         = 12.0
M'''(0) = 57.000   theory lam^3 + 3lam^2 + lam = 57.0
Var = M''(0) - M'(0)^2 = 3.000
sample E[X], E[X^2], E[X^3]: [3.0, 12.0, 56.98]
skewness 0.576   theory 1/sqrt(lam) = 0.577

import numpy as np

rng = np.random.default_rng(1)
n, t = 100_000, 0.3
cases = {"Poisson(2)": (rng.poisson(2, n), np.exp(2 * (np.exp(t) - 1))),
         "N(0, 1)":    (rng.standard_normal(n), np.exp(t**2 / 2)),
         "Exp(1)":     (rng.exponential(1, n), 1 / (1 - t))}
for name, (x, exact) in cases.items():
    est = np.mean(np.exp(t * x))
    print(f"{name:10s}  mean of e^(tX) = {est:.4f}   formula = {exact:.4f}")

# Exp(1) near the edge (t -> 1): Var(e^(tX)) = M(2t) - M(t)^2 is infinite once 2t >= 1
x = rng.exponential(1, (5, n))                      # five independent samples
for t in (0.3, 0.9):
    est = np.mean(np.exp(t * x), axis=1)
    print(f"t = {t}: five estimates {np.round(est, 2)}   formula {1 / (1 - t):.2f}")

Poisson(2)  mean of e^(tX) = 2.0104   formula = 2.0132
N(0, 1)     mean of e^(tX) = 1.0452   formula = 1.0460
Exp(1)      mean of e^(tX) = 1.4326   formula = 1.4286
t = 0.3: five estimates [1.42 1.43 1.43 1.43 1.43]   formula 1.43
t = 0.9: five estimates [6.22 7.64 8.46 7.22 7.83]   formula 10.00

À t = 0,9, M(0,9) = 10 existe, mais e^0,9X a une variance infinie, M(1,8) − M(0,9)² : les cinq estimations s'étalent de 6,2 à 8,5, toutes sous la vraie valeur : une bonne part de M(0,9) vient de tirages énormes et rares, que 100 000 tirages contiennent trop peu.

import numpy as np

rng = np.random.default_rng(6)
h = rng.normal(69, 3, 200_000)                      # true heights, inches
zone = (h >= 70.5) & (h < 72)                       # from 5'10.5" to just under 6'
bump = zone & (rng.random(h.size) < 0.5)            # half of them say "6 feet"
r = np.where(bump, 72.0, h)                         # reported heights
print(f"share reporting exactly 72 in: {np.mean(r == 72):.4f}")

for w in (2, 1, 0.5, 0.1):
    edges = 72 + w * (np.arange(-150, 151) - 0.5)   # bins centred on 72 (bin 150)
    dens, _ = np.histogram(r, edges, density=True)
    other = np.delete(dens, 150).max()
    print(f"bin width {w:>3}: bar at 72 = {dens[150]:.3f}   tallest other = {other:.3f}")

F = lambda x: np.mean(r <= x)                       # empirical CDF
print(f"F(71.999) = {F(71.999):.4f}   F(72) = {F(72):.4f}   jump = {F(72) - F(71.999):.4f}")

share reporting exactly 72 in: 0.0749
bin width   2: bar at 72 = 0.094   tallest other = 0.124
bin width   1: bar at 72 = 0.133   tallest other = 0.134
bin width 0.5: bar at 72 = 0.210   tallest other = 0.134
bin width 0.1: bar at 72 = 0.811   tallest other = 0.137
F(71.999) = 0.7668   F(72) = 0.8417   jump = 0.0749

La hauteur de la classe de 72 dépend de la largeur choisie, celle des autres non : le signe qu'il n'y a pas de densité en 72. Le saut de la fonction de répartition, lui, ne dépend d'aucun choix.

import numpy as np

rng = np.random.default_rng(2)
n = 1_000_000                                       # one-second windows
a = rng.poisson(2.0, n)                             # requests from service A
b = rng.poisson(3.0, n)                             # requests from service B, independent
s = a + b                                           # merged stream
print(f"merged:     mean {s.mean():.3f}  var {s.var():.3f}  P(0) {np.mean(s == 0):.5f}")
print(f"Poisson(5): mean 5      var 5      P(0) {np.exp(-5):.5f}")

t = 0.3
m = lambda x: np.mean(np.exp(t * x))                # empirical MGF at t
exact = np.exp(5 * (np.exp(t) - 1))
print(f"M_A * M_B = {m(a) * m(b):.4f}   M_(A+B) = {m(s):.4f}   exact = {exact:.4f}")
print(f"dependent, B = A:  M_(2A) = {m(2 * a):.4f}   but M_A^2 = {m(a)**2:.4f}")

merged:     mean 5.005  var 4.998  P(0) 0.00661
Poisson(5): mean 5      var 5      P(0) 0.00674
M_A * M_B = 5.7599   M_(A+B) = 5.7590   exact = 5.7505
dependent, B = A:  M_(2A) = 5.1836   but M_A^2 = 4.0575

import numpy as np

t = 1.0
coin = lambda u: np.cosh(u)                         # X = +1 or -1: mean 0, sigma 1
cexp = lambda u: np.exp(-u) / (1 - u)               # Exp(1) - 1: mean 0, sigma 1, skewed
print("     n   coin M_Zn(1)   exp M_Zn(1)   n * eps_n (exp)")
for n in (2, 10, 100, 1000, 10_000):
    u = t / np.sqrt(n)                              # step 4: M at t / (sigma sqrt n)
    eps = cexp(u) - 1 - t**2 / (2 * n)              # step 5: left after 1 + t^2/(2n)
    print(f"{n:6d}   {coin(u)**n:.5f}        {cexp(u)**n:.5f}       {n * eps:.5f}")
print(f"limit e^(t^2/2) = {np.exp(t**2 / 2):.5f}")

# Monte Carlo: standardized sums of 100 exponentials
rng = np.random.default_rng(3)
z = (rng.exponential(1, (200_000, 100)).sum(axis=1) - 100) / np.sqrt(100)
print(f"E[e^Z] = {np.mean(np.exp(z)):.4f}   (exact {cexp(0.1)**100:.4f})")
print(f"P(Z <= 1) = {np.mean(z <= 1):.4f}   (Phi(1) = 0.8413)")

     n   coin M_Zn(1)   exp M_Zn(1)   n * eps_n (exp)
     2   1.58909        2.83398       0.86688
    10   1.63539        1.89464       0.15989
   100   1.64735        1.70924       0.03749
  1000   1.64858        1.66662       0.01093
 10000   1.64871        1.65427       0.00337
limit e^(t^2/2) = 1.64872
E[e^Z] = 1.7160   (exact 1.7092)
P(Z <= 1) = 0.8410   (Phi(1) = 0.8413)

La pièce converge vite, l'exponentielle lentement, et n · ε_n décroît comme 1/√n : à n = 100, il vaut 0,037, près de 1/(3√100) = 0,033.