A jar holds $n=5$ distinct tokens; draw $k=3$. Count the outcomes when the draw is (a) ordered with replacement, (b) ordered without replacement, (c) unordered without replacement, (d) unordered with replacement. Show these are $n^k$, $P(n,k)=n!/(n-k)!$, $\binom{n}{k}$, and $\binom{n+k-1}{k}$, and explain the strict ordering $\binom{n}{k}\le P(n,k)\le n^k$ (removing order divides by $k!$; removing replacement removes the repeated draws). Then count the license plates of 2 letters followed by 3 digits (ordered, with replacement, mixed alphabets): $26^2\cdot10^3$. Check: ordered_with_replacement(5,3)$=125$, permutations(5,3)$=60$, combinations(5,3)$=10$, combinations_with_replacement(5,3)$=35$; multiplication_principle([26,26,10,10,10])$=676000$.

Ten sprinters race; gold, silver, and bronze medals are awarded. How many distinct podiums are possible? This is an ordered selection of $k=3$ from $n=10$ without replacement, $P(10,3)=10\cdot9\cdot8=10!/7!$. Verify the "choose then order" identity $P(n,k)=\binom{n}{k}\,k!$ by counting the same podium as: pick the 3 medalists ($\binom{10}{3}$ ways), then order them on the podium ($3!$ ways). Check: permutations(10,3)$=720$ and combinations(10,3)*factorial(3) $=120\cdot6=720$.

(a) How many 5-card hands can be dealt from a 52-card deck? Order is irrelevant, no replacement: $\binom{52}{5}$. (b) A committee of 3 is chosen from 12 people; show that choosing the 3 members ($\binom{12}{3}$) equals choosing the 9 non-members ($\binom{12}{9}$) — the symmetry $\binom{n}{k}=\binom{n}{n-k}$. Check: combinations(52,5)$=2598960$; combinations(12,3)$=$combinations(12,9)$=220$.

Ten identical candies are handed out to 4 children. (a) If a child may receive none, count the distributions: non-negative solutions of $x_1+x_2+x_3+x_4=10$, which is $\binom{10+4-1}{4-1}=\binom{13}{3}$. (b) If every child must get at least one, give one to each first and distribute the remaining $6$ freely: $\binom{6+4-1}{4-1}=\binom{9}{3}$. Note this count is exactly the ~SM-01 Einstein-solid multiplicity $\Omega=\binom{q+N-1}{q}$ with $q=10$ quanta, $N=4$ oscillators. Check: stars_and_bars(10,4)$=286$ and stars_and_bars(6,4)$=84$.

Count the distinguishable arrangements of the letters of MISSISSIPPI (1 M, 4 I, 4 S, 2 P, total 11): of the $11!$ orderings of distinct tiles, those permuting like letters coincide, an overcount of $1!\,4!\,4!\,2!$, giving $\binom{11}{1,4,4,2}= 11!/(1!\,4!\,4!\,2!)$. Do the same for STATISTICS (3 S, 3 T, 1 A, 2 I, 1 C, total 10). Recall this multinomial coefficient is the statistical multiplicity $W$ of ~SM-01. Check: multinomial(11,[1,4,4,2])$=34650$ and multinomial(10,[3,3,1,2,1])$=50400$.

(a) Build row 6 of Pascal's triangle by the rule $\binom{6}{k}=\binom{5}{k-1}+ \binom{5}{k}$ and confirm the row sums to $2^6$ (the number of subsets of a 6-set, from $(1+1)^6$); confirm the alternating sum is $0$. (b) Specialize the binomial theorem to $x=p,\,y=1-p$ to show $\sum_{k=0}^{n}\binom{n}{k}p^k(1-p)^{n-k}=1$ — this is the normalization of the binomial pmf of ~ST-07, with $\binom{n}{k}$ counting the orderings of $k$ successes in $n$ trials. Check: binomial_coefficient_sum(6)$=64$, binomial_coefficient_sum(6, signed=True)$=0$, binomial_coefficient_sum(10)$=1024$; binomial_theorem_check(0.3,0.7,6)[0]$=1.0$.

The discrete $\binom{n}{k}$ has a continuous counterpart through the Beta integral $$\int_0^1 x^{k}(1-x)^{\,n-k}\,dx=\frac{k!\,(n-k)!}{(n+1)!}=\frac{1}{(n+1)\binom{n}{k}}.$$ Take $n=5,\ k=2$: predict $1/[(6)\binom{5}{2}]=1/60$, then integrate numerically. Separately confirm the factorial is the Gamma integral $\Gamma(6)=\int_0^\infty t^5 e^{-t}dt=5!$. Check: beta_integral(2,5)$=1/60\approx0.0166667$ and _integrate(lambda x: x2(1-x)3, 0, 1)$\approx0.0166667$; _integrate(lambda t: t**5exp(-t), 0, 60)$\approx120=$factorial(5).

Two fair coins are tossed. Using $P(A\mid B)=P(A\cap B)/P(B)$, find (a) the probability the first toss is heads given the second is heads, and (b) the probability the two tosses agree ($C$) given the first is heads ($A$). For (a), $P(A\cap B)=P(HH)=\tfrac14$ and $P(B)=\tfrac12$, so $P(A\mid B)=\tfrac14/\tfrac12 =\tfrac12$. For (b), $C\cap A=\{HH\}$ so $P(C\mid A)=\tfrac{1/4}{1/2}=\tfrac12$. In both cases the conditional equals the marginal $\tfrac12$ — heralding independence (P3). Check: m,A,B,C = two_coin_space(); conditional_event(m,A,B) $=0.5$ and conditional_event(m,C,A) $=0.5$.

Draw two cards from a $52$-card deck without replacement. By the multiplication rule $P(\text{ace}_1\cap\text{ace}_2)=P(\text{ace}_1)\,P(\text{ace}_2\mid \text{ace}_1)$. The first factor is $4/52$; given an ace is gone, $3$ aces remain among $51$ cards, so the second factor is $3/51$. Hence $$P(\text{two aces})=\frac{4}{52}\cdot\frac{3}{51}=\frac{12}{2652}=\frac{1}{221} \approx0.004525 .$$ The second factor differs from $4/52$, so the draws are dependent. Check: deck,fa,sa = two_card_deck(); chain_rule_factors(deck,[fa,sa]) $=[0.07692\ldots, 0.05882\ldots]$ (i.e. $4/52,\,3/51$); chain_rule(_) $=0.0045249=$ prob_inter(deck,fa,sa); and independent_events(deck,fa,sa) is False.

Show the two fair-coin events $A$ (first $H$) and $B$ (second $H$) are independent in both senses: the product rule $P(A\cap B)=\tfrac14=\tfrac12\cdot\tfrac12 =P(A)P(B)$, and the conditional rule $P(A\mid B)=\tfrac12=P(A)$. Then contrast with $A$ and its complement $A^{c}$ (first toss tails): these are disjoint, so $P(A\cap A^{c})=0\ne P(A)P(A^{c})=\tfrac14$ — disjoint events of positive probability are the opposite of independent. Check: independent_events(m,A,B) and conditional_equals_marginal(m,A,B) are both True; with Ac the complement of A, independent_events(m,A,Ac) is False (and prob_inter(m,A,Ac) $=0$).

For the two coins let $C=$ "the tosses agree." Show $A,B,C$ are pairwise independent — $P(A\cap B)=P(A\cap C)=P(B\cap C)=\tfrac14=$ (product of the two marginals each $\tfrac12$) — yet not mutually independent, because $A$ and $B$ together force $C$: $$P(A\cap B\cap C)=P(HH)=\tfrac14\ \ne\ P(A)P(B)P(C)=\tfrac18 .$$ This is why mutual independence ($2^{n}-n-1$ equations) is strictly stronger than pairwise ($\binom n2$ equations). Check: pairwise_independent(m,[A,B,C]) $=$ True while mutually_independent(m,[A,B,C]) $=$ False; prob(m, set(A)&set(B)&set(C)) $=0.25$ but prob(m,A)prob(m,B)prob(m,C) $=0.125$.

Fix $B=$ "second toss $H$" and define $Q(\cdot)=P(\cdot\mid B)$. Verify $Q$ obeys the Kolmogorov axioms (~ST-01): it is nonnegative, $Q(\Omega)=Q(B)=1$, and the outcomes outside $B$ carry zero mass, $Q(B^{c})=0$. Concretely the two outcomes in $B$, $HH$ and $TH$, each get reweighted to $P(\{o\})/P(B)=\tfrac{1/4}{1/2} =\tfrac12$, summing to $1$. Check: Q = conditional_measure(m,B); is_probability_measure(Q) $=$ True; prob(Q,B) $=1.0$; Q[('H','H')] $=0.5$; prob(Q,[o for o in m if o not in set(B)]) $=0.0$.

Toss three independent fair bits and let $A,B,C$ be "bit $2$ is $1$," "bit $3$ is $1$," "bit $1$ is $1$." Inside the world where $C$ holds, $A$ and $B$ are still independent fair bits, so $P(A\cap B\mid C)=P(A\mid C)P(B\mid C)$. Compute $P(A\mid C)=P(B\mid C)=\tfrac12$ and $P(A\cap B\mid C)=\tfrac14=\tfrac12\cdot \tfrac12$ — conditional independence holds. Check: build the space outs=[(c,a,b) for c in (0,1) for a in (0,1) for b in (0,1)], mm=normalize({o:1.0 for o in outs}), with C,A,B the bit-indicators; then conditionally_independent_check(mm,A,B,C) $=$ True, P(A|C)=0.5, P(A∩B|C)=0.25.

Let $T\sim\text{Exp}(\lambda)$ with $\lambda=0.7$, survival $P(T>t)=e^{-\lambda t}$. Show the conditional "extra wait" forgets the elapsed time: for $s=2,\,t=3$, $$P(T>s+t\mid T>s)=\frac{P(T>5)}{P(T>2)}=\frac{e^{-0.7\cdot5}}{e^{-0.7\cdot2}} =e^{-0.7\cdot3}=e^{-2.1}\approx0.122456=P(T>t),$$ independent of $s$. The exponential is the unique continuous law with this property (continuous cousin of the geometric, ~ST-08). Check: memoryless_conditional(2.0,3.0,0.7) $\approx 0.122456$ $=$ exponential_survival(3.0,0.7) $= e^{-2.1}$.

Let $X,Y$ be independent $\text{Uniform}(0,1)$, joint density $f(x,y)=1$ on the unit square. Show the product event factors, $$P(X<a,\,Y<b)=\int_0^a\!\!\int_0^b 1\,dy\,dx=ab=P(X<a)\,P(Y<b),$$ so for $a=0.6,\,b=0.4$ the joint is $0.24$, and the conditional equals the marginal, $P(X<a\mid Y<b)=ab/b=a=0.6$. This factorization is independence of random variables (~ST-13). Check: prob_rect(uniform_square_joint,0.0,0.6,0.0,0.4) $\approx0.24$ and conditional(0.24, 0.4) $=0.6$.

Three factories supply a parts bin: factory 1 makes $50\%$ of the parts, factory 2 makes $30\%$, factory 3 makes $20\%$, with defect rates $1\%$, $2\%$, $3\%$ respectively. (a) Find $P(\text{defective})$ by the law of total probability $P(B)=\sum_i P(B\mid A_i)P(A_i)$. (b) A part is found defective — which factory most likely made it? Use $P(A_i\mid B)=P(B\mid A_i)P(A_i)/P(B)$. Note the posterior is not ordered like the defect rates, because factory 1's large share offsets its low defect rate. Check: total_probability([0.5,0.3,0.2],[0.01,0.02,0.03]) $=0.017$; bayes_posterior([0.5,0.3,0.2],[0.01,0.02,0.03]) $=[0.2941,0.3529,0.3529]$ (factories 2 and 3 tie as most likely).

Urn 1 holds 2 red and 3 blue balls; urn 2 holds 4 red and 1 blue. You pick an urn at random (prior $\tfrac12,\tfrac12$) and draw a red ball. Show $P(\text{red}\mid U_1)=\tfrac25$, $P(\text{red}\mid U_2)=\tfrac45$, so by total probability $P(\text{red})=\tfrac12\cdot\tfrac25+\tfrac12\cdot\tfrac45=\tfrac35$, and by Bayes $P(U_1\mid\text{red})=\dfrac{\tfrac12\cdot\tfrac25}{\tfrac35}=\tfrac13$. The red draw shifts belief from $\tfrac12$ toward the red-richer urn 2. Check: total_probability([0.5,0.5],[0.4,0.8]) $=0.6$; bayes_posterior([0.5,0.5],[0.4,0.8]) $=[0.3333,0.6667]$.

A disease has prevalence $\pi=0.001$. A test has sensitivity $se=P(+\mid D)=0.99$ and specificity $sp=P(-\mid H)=0.99$. (a) Compute the positive predictive value $\mathrm{PPV}=P(D\mid +)=\dfrac{se\,\pi}{se\,\pi+(1-sp)(1-\pi)}$ and explain why a $99\%$-accurate test gives only $\approx 9\%$ — the base-rate effect: false positives from the huge healthy majority swamp the few true positives. (b) Show PPV $=\tfrac12$ exactly when $\pi=1-sp$ (here $\pi=0.01$). (c) Compute the negative predictive value and note it is essentially $1$. Check: ppv(0.001,0.99,0.99) $=0.0902$; ppv(0.01,0.99,0.99) $=0.5$; npv(0.001,0.99,0.99) $\approx 0.99999$.

Redo P3 in odds. The prior odds of disease are $\pi/(1-\pi)=0.001/0.999\approx 0.001$; the positive likelihood ratio is $LR^+=se/(1-sp)=0.99/0.01=99$. Multiply: posterior odds $=0.001\times 99\approx 0.099$, hence $P=\dfrac{0.099}{1.099}\approx 0.090$ — the same PPV, with the evidence $P(+)$ never computed because it cancels in the odds. Check: lr_positive(0.99,0.99) $=99$; odds_to_prob(prob_to_odds(0.001)*lr_positive(0.99,0.99)) $=0.0902$ $=$ ppv(0.001,0.99,0.99).

The patient retests. Treating the tests as conditionally independent given disease status, each positive multiplies the odds by $LR^+=99$. (a) After two positives the odds are $0.001\times 99^2\approx 9.81$, so $P\approx 0.908$; after three, $0.001\times 99^3\approx 971$, so $P\approx 0.999$ — accumulating evidence overturns the tiny prior. (b) Show a positive then a negative returns belief to the base rate, because $LR^+\!\cdot LR^-=\dfrac{se}{1-sp}\cdot\dfrac{1-se}{sp}=1$ here. Check: with L = np.array([[0.01,0.99],[0.99,0.01]]), sequential_update([0.001,0.999], L, [1,1,1])[:,0] $=[0.001,0.0902,0.9075,0.999]$; sequential_update([0.001,0.999], L, [1,0])[:,0] $=[0.001,0.0902,0.001]$ (a $+$ and a $-$ cancel).

A coin of unknown bias $\theta$ gets a uniform prior $\mathrm{Beta}(1,1)$. You flip $8$ heads in $10$. Using conjugacy, the posterior is $\mathrm{Beta}(1+8,\,1+2)=\mathrm{Beta}(9,3)$, with mean $\tfrac{9}{12}=0.75$ — between the prior mean $0.5$ and the MLE $\tfrac{8}{10}=0.8$ (Bayesian shrinkage). The evidence $P(8\mid 10)=\int_0^1 1\cdot\binom{10}{8}\theta^8(1-\theta)^2\,d\theta =\binom{10}{8}B(9,3)=\tfrac{1}{11}$, the uniform-prior result that all $11$ counts are equally likely a priori. Check: posterior_beta_binomial(1,1,8,10) $=(9,3)$; beta_mean(9,3) $=0.75$; evidence_beta_binomial(1,1,8,10) $=0.090909=\tfrac{1}{11}$ (and equals the midpoint integral of beta_pdf(t,1,1)*binomial_likelihood(t,8,10)).

For the rare disease of P3–P5, find the number $m$ of independent positive tests that first pushes $P(D)$ above $\tfrac12$, i.e. posterior odds $>1$. Solve $\pi/(1-\pi) \cdot (LR^+)^m>1$: $m>\dfrac{\ln\!\big((1-\pi)/\pi\big)}{\ln LR^+}=\dfrac{\ln 999} {\ln 99}\approx 1.50$, so $m=2$ positives are needed (one is not enough). This is the weight-of-evidence view: each positive adds $\log LR^+\approx 4.6$ nats of log-odds, and you need to overcome $\ln 999\approx 6.9$ of prior skepticism. Check: with L as in P5, sequential_update([0.001,0.999], L, [1,1])[:,0] $=[0.001,0.0902,0.9075]$ — below $0.5$ after one positive, above after two.

Let $X\sim\text{Uniform}\{1,\dots,6\}$, so $f(x)=\tfrac16$ on its support. Verify it is a valid pmf ($\sum f=1$); compute the cdf value $F(3)=P(X\le3)$; the mean $\mu=\tfrac16(1+\cdots+6)=\tfrac{21}{6}$; the second raw moment $E[X^2]=\tfrac16(1+4+\cdots+36)=\tfrac{91}{6}$; and hence $\sigma^2=E[X^2]-\mu^2=\tfrac{91}{6}-\tfrac{49}{4}=\tfrac{35}{12}$. Check: with x=np.arange(1,7); p=np.full(6,1/6), pmf_is_valid(x,p) $=$ True; cdf_from_pmf(x,p)(3.0) $=0.5$; expectation(x,p) $=3.5$; raw_moment(x,p,2) $=15.1\overline{6}$; variance(x,p) $=2.91\overline{6}$; std(x,p) $\approx1.7078$.

For the same die, compute $E[(X-\mu)^2]$ directly by LOTUS with $g(x)=(x-3.5)^2$ — summing $g(x)\,f(x)$ over $x=1,\dots,6$ — and confirm it equals the variance from P1. This is the definitional route $\sigma^2=E[(X-\mu)^2]$; P1 used the computational route $E[X^2]-\mu^2$. (They are forced equal by the algebra of §5.) Check: expectation_of(lambda t:(t-3.5)**2, x, p) $=2.91\overline{6}=$ variance(x,p); equivalently central_moment(x,p,2) $=2.91\overline{6}$.

Let $X$ be Bernoulli$(p)$: $f(1)=p$, $f(0)=1-p$. Show $\mu=E[X]=0\cdot(1-p)+1\cdot p =p$, and since $X^2=X$ here, $E[X^2]=p$ too, so $\sigma^2=E[X^2]-\mu^2=p-p^2=p(1-p)$ — maximal at $p=\tfrac12$. This single trial is the atom the binomial of ~ST-07 is built from. Check: with bx=np.array([0.,1]); bp=np.array([.7,.3]) ($p=0.3$), expectation(bx,bp) $=0.3$ and variance(bx,bp) $=0.21=p(1-p)$.

Let $X$ take values $\{-1,0,2\}$ with $f=\{0.2,0.5,0.3\}$. Compute $\mu=E[X]= -0.2+0+0.6=0.4$ and $E[X^2]=0.2(1)+0.5(0)+0.3(4)=1.4$, hence $\sigma^2=1.4-(0.4)^2=1.24$ and $\sigma=\sqrt{1.24}\approx1.1136$. Verify the computational formula $E[X^2]-\mu^2$ agrees with the central moment $E[(X-\mu)^2]$. Check: with cx=np.array([-1.,0,2]); cp=np.array([.2,.5,.3]), expectation(cx,cp) $=0.4$; raw_moment(cx,cp,2) $=1.4$; variance(cx,cp) $=1.24$ $=$ central_moment(cx,cp,2); std(cx,cp) $\approx1.1136$.

For the $X$ of P4, let $Y=3X-4$. Use $E[aX+b]=aE[X]+b$ to get $E[Y]=3(0.4)-4=-2.8$, and $\mathrm{Var}(aX+b)=a^2\mathrm{Var}(X)$ to get $\mathrm{Var}(Y)=9(1.24)=11.16$, so $\sigma_Y=3\sigma_X\approx3.341$. Note the shift $b=-4$ moves the mean but leaves the variance untouched. Check: yx,yp = linear_transform(3.,-4.,cx,cp); expectation(yx,yp) $=-2.8$; variance(yx,yp) $=11.16$; std(yx,yp) $\approx3.341$. Changing $b$ alone (e.g. linear_transform(3.,96.,cx,cp)) leaves variance at $11.16$.

Show the first central moment always vanishes, $\mu_1=E[X-\mu]=0$, and the zeroth raw moment is $\mu'_0=1$. Then compare skewness $\gamma_1=\mu_3/\sigma^3$: the fair die is symmetric about $3.5$, so $\gamma_1=0$; the RV of P4 is right-skewed (mass piled at $-1,0$ with a far atom at $2$), giving $\gamma_1>0$. Compute $\mu_3=\sum(x-0.4)^3 f=0.648$ and $\sigma^3=1.24^{3/2}\approx1.3808$, so $\gamma_1\approx0.469$. Check: central_moment(cx,cp,1) $=0$; raw_moment(cx,cp,0) $=1$; skewness(x,p) $=0.0$ (die); skewness(cx,cp) $\approx0.4693$.

Let $X$ take values $\{0,1,2,3\}$ with $f=\{0.1,0.2,0.3,0.4\}$. Compute the mean two ways: directly $\mu=\sum x f(x)=0.2+0.6+1.2=2.0$; and by the tail-sum identity for a non-negative integer RV, $E[X]=\sum_{k\ge0}\big(1-F(k)\big)$ with $1-F(0)=0.9$, $1-F(1)=0.7$, $1-F(2)=0.4$ (and $0$ thereafter), summing to $2.0$. Check: with x=np.array([0.,1,2,3]); p=np.array([.1,.2,.3,.4]), expectation(x,p) $=2.0$ and mean_via_survival(x,p) $=2.0$; cdf_from_pmf(x,p)(1.0) $=0.3$.

For Bernoulli$(p)$ the moment-generating function is $M(t)=E[e^{tX}]=(1-p)+p\,e^{t}$. Check the normalization $M(0)=1$; differentiate to get $M'(t)=p\,e^t$, so $M'(0)=p=\mu$; and $M''(0)=p=E[X^2]$, recovering $\sigma^2=M''(0)-M'(0)^2=p-p^2$. This is the §8 preview made explicit, and the method ~ST-06 generalizes. Check: with the $p=0.3$ Bernoulli, mgf(bx,bp,0.0) $=1.0$; the finite difference (mgf(bx,bp,1e-6)-mgf(bx,bp,-1e-6))/2e-6 $\approx0.3=$ expectation(bx,bp).

From $M(t)=\sum_x e^{tx}p(x)$ show $M(0)=\sum_x p(x)=1$ for any distribution. Then expand $e^{tX}=\sum_k(tX)^k/k!$ inside the expectation to get $M(t)=\sum_k E[X^k]\,t^k/k!$, and conclude $M'(0)=E[X]=\mu$ and, in general, $M^{(k)}(0)=E[X^k]$. Check: mgf(0.0, binomial_dist(10,0.3)) $=1.0$ and mean_from_mgf(binomial_dist(10,0.3)) $\approx 3.0$ ($=np$).

For $X\sim\text{Bin}(n,p)$ the mgf is $M(t)=\big((1-p)+pe^t\big)^n$ (a product of $n$ Bernoulli mgfs). Differentiate: $M'(t)=n((1-p)+pe^t)^{n-1}pe^t$ so $M'(0)=np$; compute $M''(0)=n(n-1)p^2+np$ and hence $$\sigma^2=M''(0)-[M'(0)]^2=np-np^2=np(1-p).$$ Check: with d = binomial_dist(10,0.3), mean_from_mgf(d) $\approx 3.0$ and var_from_mgf(d) $\approx 2.1=np(1-p)$ (both agree with mean(d), variance(d)).

For $X\sim\text{Poi}(\lambda)$, $M(t)=e^{\lambda(e^t-1)}$, so the cumulant generator is $K(t)=\ln M(t)=\lambda(e^t-1)$. Show $K^{(k)}(0)=\lambda$ for every $k\ge1$, so mean $=$ variance $=\lambda$ and indeed all cumulants equal $\lambda$. Using $m_k=\sum$(cumulant compositions), derive the third raw moment $E[X^3]=\lambda^3+3\lambda^2+\lambda$. Check: with dp = poisson_dist(2.0), [cumulant_from_mgf(k,dp) for k in (1,2,3)] $\approx[2,2,2]$ and moment_from_mgf(3,dp) $\approx 22.0=2^3+3\cdot2^2+2$.

For independent $X,Y$ prove $M_{X+Y}(t)=M_X(t)M_Y(t)$ by factoring the expectation of $e^{tX}e^{tY}$. Apply it to two independent $\text{Bernoulli}(0.4)$ variables: their sum has mgf $\big(0.6+0.4e^t\big)^2$, the $\text{Bin}(2,0.4)$ mgf. Confirm the distribution of the sum (the convolution $[0.36,0.48,0.16]$ on $\{0,1,2\}$) carries that same mgf. Check: at $t=0.5$, with b = bernoulli_dist(0.4), mgf_of_sum(0.5,[b,b]), binomial_mgf(0.5,2,0.4), and mgf(0.5, *convolve_dists(b,b)) are all $\approx 1.5863113$.

Let $X$ be the trial of the first success, $\text{Geom}(p)$ on $\{1,2,\dots\}$. Sum the geometric series to get $M(t)=\dfrac{pe^t}{1-(1-p)e^t}$ for $t<-\ln(1-p)$. Differentiate at $0$ to obtain $\mu=1/p$ and $\sigma^2=(1-p)/p^2$. Check: with g = geometric_dist(0.25), mean_from_mgf(g) $\approx 4.0=1/p$ and var_from_mgf(g) $\approx 12.0 =(1-0.25)/0.25^2$.

Using $M_N(t)=e^{\mu t+\frac12\sigma^2t^2}$, show that for independent $N(\mu_1,\sigma_1^2)$ and $N(\mu_2,\sigma_2^2)$ the product of mgfs is again a normal mgf with mean $\mu_1+\mu_2$ and variance $\sigma_1^2+\sigma_2^2$ — so by uniqueness the sum is normal. Separately verify the rule $M_{aX+b}(t)=e^{bt}M_X(at)$. Check: at $t=0.3$, normal_mgf(0.3,1,2)normal_mgf(0.3,-0.5,1.5) equals normal_mgf(0.3,0.5,math.sqrt(22+1.52)) ($\approx 1.5391803$); and with d = binomial_dist(6,0.5), v,p=d, mgf(0.3, 2v+3, p) equals math.exp(30.3)mgf(0.6,*d) ($\approx 19.41497$).

For $X\sim\text{Exp}(\lambda)$ with $f(x)=\lambda e^{-\lambda x}$ ($x\ge0$), evaluate $M(t)=\int_0^\infty e^{tx}\lambda e^{-\lambda x}dx=\dfrac{\lambda}{\lambda-t}$ for $t<\lambda$, and read off $\mu=M'(0)=1/\lambda$, $\sigma^2=1/\lambda^2$. Note $M(t)=\mathcal{L}\{f\}(-t)$: the mgf is the (two-sided) Laplace transform of the density at $s=-t$ (~MA-10). Check: mgf_continuous(0.5, lambda x:1.5math.exp(-1.5x), 0.0, 80.0) $\approx 1.5=$ exponential_mgf(0.5,1.5); differentiating the closed form gives mean $\approx 0.6667=1/1.5$.

Show that a sum of $\alpha$ independent $\text{Exp}(\lambda)$ variables has mgf $\big(\lambda/(\lambda-t)\big)^{\alpha}$, the $\text{Gamma}(\alpha,\lambda)$ mgf, so by uniqueness the sum is $\text{Gamma}(\alpha,\lambda)$ (the Erlang law). Hence its mean is $\alpha/\lambda$ and variance $\alpha/\lambda^2$ — $\alpha$ times the exponential's, as additivity demands. Check: gamma_mgf(0.7,4,1.5) equals exponential_mgf(0.7,1.5)**4 ($\approx 12.359619$).

For a single Bernoulli($p$) trial with indicator $X\in\{0,1\}$, derive $E[X]=p$ and $\mathrm{Var}(X)=p(1-p)$ directly from the ~ST-05 definitions, and show the variance is maximized at $p=\tfrac12$ (differentiate $p(1-p)$) and vanishes at $p=0,1$. Explain in one sentence why a fair coin is the "most uncertain" trial. Check: bernoulli_var(0.5) $=0.25$ (the maximum); bernoulli_var(0.3) $=0.21<0.25$; and bernoulli_mean(p) == binom_mean(1, p), bernoulli_var(p) == binom_var(1, p) for any $p$ (the $n=1$ binomial).

A biased coin lands heads with $p=0.3$; flip it $n=5$ times. Using $P(X=k)=\binom nk p^kq^{n-k}$, compute $P(X=2)$ by hand: there are $\binom52=10$ arrangements, each of probability $0.3^2\,0.7^3$. Then write $P(X\le2)$ as a sum of three terms. Confirm $\sum_{k=0}^5 P(X=k)=1$ is the binomial theorem $(0.3+0.7)^5$. Check: binom_pmf(2, 5, 0.3) $=0.3087$; binom_cdf(2, 5, 0.3) $=0.83692$; sum(binom_pmf(k, 5, 0.3) for k in range(6)) $=1.0$.

Write $X=\sum_{i=1}^{12}X_i$ as a sum of $n=12$ independent Bernoulli($p=0.25$) indicators. Use linearity of expectation to get $E[X]=np$ and independence (variances add) to get $\mathrm{Var}(X)=np(1-p)$, without summing $k\,P(X=k)$. Then verify the answers equal the definitional sums. Check: binom_mean(12, 0.25) $=3.0$, binom_var(12, 0.25) $=2.25$, binom_std(12, 0.25) $=1.5$; and the brute sums sum(kbinom_pmf(k,12,0.25) for k in range(13)) $=3.0$, sum((k-3)**2binom_pmf(k,12,0.25) for k in range(13)) $=2.25$.

For $X\sim\mathrm{Bin}(6,\tfrac12)$ the mgf is $M(t)=(\tfrac12+\tfrac12e^t)^6$. Show $M(0)=1$; differentiate to get $M'(0)=np=3$; and from $M''(0)=E[X^2]=np(1-p)+(np)^2$ recover $\mathrm{Var}(X)=M''(0)-M'(0)^2=np(1-p)=1.5$. Check: with M = lambda t: binom_mgf(t, 6, 0.5), _deriv1(M, 0.0) $=3.0$, _deriv2(M, 0.0) $=10.5=E[X^2]$, and _deriv2(M,0.0) - _deriv1(M,0.0)**2 $=1.5=$ binom_var(6, 0.5).

Let $X_1\sim\mathrm{Bin}(3,0.4)$ and $X_2\sim\mathrm{Bin}(4,0.4)$ be independent. Using the mgf product $(q+pe^t)^3(q+pe^t)^4=(q+pe^t)^7$, argue $X_1+X_2\sim \mathrm{Bin}(7,0.4)$, and confirm the pmf-level statement is Vandermonde's identity $\sum_j\binom3j\binom4{5-j}=\binom75$. Why does the argument fail if $p_1\ne p_2$? Check: sum_two_binomials_pmf(5, 3, 4, 0.4) $=0.0774144=$ binom_pmf(5, 7, 0.4); and the means/variances add: binom_mean(3,0.4)+binom_mean(4,0.4) == binom_mean(7,0.4).

Using the successive-ratio test $\dfrac{P(X=k)}{P(X=k-1)}=\dfrac{n-k+1}{k}\dfrac pq$, show the pmf rises while $k<(n+1)p$ and find the mode $\lfloor(n+1)p\rfloor$. (a) For $\mathrm{Bin}(10,0.35)$ compute the mode. (b) For $\mathrm{Bin}(9,\tfrac12)$ note $(n+1)p=5$ is an integer, so $k=4$ and $k=5$ are both modes — verify their pmfs are equal. Check: binom_mode(10, 0.35) $=3$; binom_mode(9, 0.5) $=5$ with binom_pmf(4, 9, 0.5) == binom_pmf(5, 9, 0.5) $=0.246094$ (the tie).

Take $n\to\infty$ with $p=\lambda/n$ fixed at mean $\lambda=2.5$. Using $\big(1-\tfrac\lambda n\big)^n\to e^{-\lambda}$ and $\tfrac{n!}{(n-k)!\,n^k}\to1$, show $\binom nk p^kq^{n-k}\to \dfrac{\lambda^k e^{-\lambda}}{k!}$. Evaluate the $k=3$ term for $n=2500$ and compare to $\mathrm{Poisson}(2.5)$. Check: binom_pmf(3, 2500, 0.001) $=0.213881$ vs _poisson_pmf(3, 2.5) $=0.213763$ (difference $\sim10^{-4}$, shrinking like $1/n$).

For $\mathrm{Bin}(100,\tfrac12)$ (so $np=50$, $npq=25$, $\sigma=5$), use the de Moivre–Laplace formula $P(X=k)\approx\frac{1}{\sqrt{2\pi npq}}e^{-(k-np)^2/2npq}$ to estimate $P(X=50)$ and compare to the exact binomial value. State the continuity correction you would use to approximate $P(48\le X\le 52)$. Check: binom_pmf(50, 100, 0.5) $=0.079589$ vs _normal_pdf(50, 50.0, 5.0) $=0.079788$ (agreement to $\sim10^{-3}$ at the peak).

For independent Bernoulli$(p)$ trials, argue $P(X=x)=(1-p)^{x-1}p$ ("$x-1$ failures, then a success") and verify $\sum_{x\ge1}f(x)=1$ using the geometric series $\sum_{k\ge0}q^k=1/(1-q)$. With $p=0.3$ evaluate $f(1)$ and $f(3)=q^2p$, and note that $f$ is strictly decreasing, so the mode is always $x=1$. Check: geometric_pmf(1, 0.3) $=0.3$; geometric_pmf(3, 0.3) $=0.147$ ($=0.7^2\cdot0.3$); sum(geometric_pmf(x,0.3) for x in range(1,2000)) $\approx1$.

Sum the pmf to get $P(X\le x)=1-q^{x}$, hence the survival $P(X>x)=q^{x}=(1-p)^{x}$ — interpret the latter directly as "the first $x$ trials all fail." For $p=0.2$ compute $P(X>5)$ and $P(X\le 5)$, and confirm $P(X\le 5)+P(X>5)=1$. Check: geometric_sf(5, 0.2) $=0.32768$ ($=0.8^5$); geometric_cdf(5, 0.2) $=0.67232$; their sum $=1.0$; geometric_cdf(10, 0.2) $=0.8926258176$.

Differentiate $\sum_{x\ge1}q^{x}=q/(1-q)$ to get $\sum x q^{x-1}=1/p^2$, hence $E[X]=1/p$; use the second factorial moment $\sum x(x-1)q^{x-2}=2/p^3$ to get $\operatorname{Var}(X)=q/p^2=(1-p)/p^2$. Evaluate both at $p=0.25$ and note the mean is the reciprocal of the success probability. Check: geometric_mean(0.25) $=4.0$; geometric_var(0.25) $=12.0$ ($=0.75/0.0625$).

Using the survival form of P2, prove $P(X>m+n\mid X>m)=q^{m+n}/q^{m}=q^{n}=P(X>n)$: a coin that has already shown $m$ tails is, going forward, as fresh as a new coin. Take $p=0.2$, $m=10$, $n=4$ and show both sides equal $0.8^4$. Check: memoryless_check(10, 4, 0.2) $=(0.4096,\,0.4096)$ (equal); both equal $0.8^4=0.4096$.

Sum the geometric series in $e^t$ to get $M(t)=pe^t/(1-qe^t)$ for $t<-\ln q$, and show $M(0)=1$ and $M'(0)=1/p=E[X]$. Evaluate $M(0.1)$ at $p=0.4$ and confirm the numerical derivative at $0$ reproduces the mean $2.5$. Check: geometric_mgf(0.0, 0.4) $=1.0$; geometric_mgf(0.1, 0.4) $\approx1.31217$; (geometric_mgf(1e-6,0.4)-geometric_mgf(-1e-6,0.4))/2e-6 $\approx2.5=$ geometric_mean(0.4).

For the trial of the $r$-th success, argue $f(x)=\binom{x-1}{r-1}(1-p)^{x-r}p^{r}$ (the $x$-th trial succeeds; $r-1$ of the first $x-1$ succeed). With $r=3$, $p=0.25$ evaluate $f(3)$ (the fastest possible: SSS) and $f(5)=\binom{4}{2}q^{2}p^{3}$. Check: negbinom_pmf(3, 3, 0.25) $=0.015625$ ($=0.25^3$); negbinom_pmf(5, 3, 0.25) $\approx0.052734$ ($=6\cdot0.75^2\cdot0.25^3$).

Write $X=G_1+\cdots+G_r$ with $G_i$ i.i.d. geometric$(p)$ (the process restarts after each success). Because means and variances of independent sums add, deduce $E[X]=r/p$ and $\operatorname{Var}(X)=r(1-p)/p^2$. Verify at $r=4$, $p=0.5$ that these equal $r$ times the geometric values. Check: negbinom_mean(4, 0.5) $=8.0=4\cdot$geometric_mean(0.5); negbinom_var(4, 0.5) $=8.0=4\cdot$geometric_var(0.5).

Show $\binom{x-1}{0}q^{x-1}p=q^{x-1}p$, so the negative binomial with $r=1$ is the geometric. Then, since the mgf of a sum is the product of mgfs, show $M_{\text{NB}}(t)=\big(pe^t/(1-qe^t)\big)^r$ equals the geometric mgf raised to $r$. Verify both with $p=0.3$. Check: negbinom_pmf(4, 1, 0.3) $=0.1029=$ geometric_pmf(4, 0.3); negbinom_mgf(0.1, 2, 0.3) $\approx2.144983=$ geometric_mgf(0.1, 0.3)**2.

A call center receives an average of $\lambda=3$ calls per minute, the calls independent. Write the pmf of the number $X$ of calls in a one-minute window and evaluate $P(X=0),P(X=1),P(X=2),P(X=3)$ by hand. Verify the recurrence $P(X{=}k{+}1)=\frac{\lambda}{k+1}P(X{=}k)$ links your four numbers, and note why $P(X{=}2)=P(X{=}3)$ here (the ratio $\lambda/(k+1)=3/3=1$ at $k=2$). Check: poisson_pmf(0,3.0)$=0.049787$, poisson_pmf(1,3.0)$=0.149361$, poisson_pmf(2,3.0)$=0.224042$, poisson_pmf(3,3.0)$=0.224042$ (equal, as argued).

Flaws occur in an optical fiber at rate $\lambda=1.5$ per km. Find the probability that a 1-km length has at least one flaw. Use the complement $P(X\ge1)=1-P(X=0)=1-e^{-\lambda}$ rather than summing the whole tail — the single most common Poisson manoeuvre. Check: poisson_pmf(0,1.5)$=0.223130$ so $P(X\ge1)=1-0.223130=0.776870$ $=$ 1 - poisson_pmf(0,1.5).

Particles strike a detector at $\lambda=4$ per second. Find $P(X\le2)$ and $P(X\ge3)$ for a one-second count. Sum $P(0)+P(1)+P(2)$ for the first; use $P(X\ge3)=1-P(X\le2)$ for the second. Check: poisson_cdf(2,4.0)$=0.238103$ and 1 - poisson_cdf(2,4.0)$=0.761897$.

Starting from $M(t)=\exp\!\big(\lambda(e^{t}-1)\big)$, show $M'(t)=\lambda e^{t}M(t)$, hence $E[X]=M'(0)=\lambda$, and $M''(0)=\lambda+\lambda^{2}$, hence $\mathrm{Var}(X)=M''(0)-M'(0)^{2}=\lambda$. Confirm the universal fact $M(0)=1$. Take $\lambda=5$. Check: poisson_mgf(0.0,5.0)$=1.0$; the central-difference derivatives of poisson_mgf give $M'(0)=5.0=$poisson_mean(5.0) and $M''(0)-M'(0)^2=5.0=$poisson_var(5.0) (see test_mgf_value_and_moment_derivatives).

A page has $n=100$ characters, each mistyped independently with probability $p=0.02$. The number of typos is $\mathrm{Bin}(100,0.02)$; approximate it by $\mathrm{Poisson}(\lambda)$ with $\lambda=np=2$ and compare $P(X=3)$ exactly and approximately. Explain why the approximation is good (large $n$, small $p$, moderate $\lambda$). Check: binomial_pmf(3,100,0.02)$=0.182276$ vs poisson_pmf(3,2.0)$=0.180447$; the worst-case gap poisson_limit_of_binomial(100,0.02)$=2.743\times10^{-3}$, and it shrinks $\approx10\times$ each time $n$ grows $10\times$ at fixed $\lambda$.

Red cars pass at $\lambda_1=3$ per minute and blue cars at $\lambda_2=4$ per minute, independently. Show via the mgf that the total $X_1+X_2$ is $\mathrm{Poisson}(\lambda_1+\lambda_2)=\mathrm{Poisson}(7)$, and find the probability of exactly $4$ cars in a minute. Confirm the pmf-level convolution $\sum_{j=0}^{4}P(X_1{=}j)P(X_2{=}4{-}j)$ gives the same value (binomial theorem). Check: sum_of_poissons([3.0,4.0])$=7.0$; poisson_convolution_pmf(4,3.0,4.0)$=0.091226=$poisson_pmf(4,7.0).

For $\lambda=9$ give $\sigma$, the skewness $1/\sqrt{\lambda}$ and the excess kurtosis $1/\lambda$. Comment on why both shape measures $\to0$ as $\lambda$ grows and what that says about a normal approximation (~ST-18). Check: poisson_var(9)$=9.0$, poisson_std(9)$=3.0$, poisson_skewness(9)$=0.3333$, poisson_excess_kurtosis(9)$=0.1111$; both shape numbers vanish like $\lambda^{-1/2}$ and $\lambda^{-1}$.

For a unit-rate process, "fewer than $k$ events by time $\lambda$" equals "the $k$-th arrival is after $\lambda$." Use this to argue the identity $$\sum_{j=0}^{k-1}\frac{e^{-\lambda}\lambda^{j}}{j!} =\int_{\lambda}^{\infty}\frac{t^{\,k-1}e^{-t}}{(k-1)!}\,\mathrm{d}t,$$ the survival function of a $\mathrm{Gamma}(k,1)$ arrival time. Verify it numerically for $\lambda=3,k=4$ (left side $=P(X\le3)$). Check: poisson_cdf(3,3.0)$=0.647232$ equals the midpoint integral $\int_{3}^{\infty}t^{3}e^{-t}/6\,\mathrm{d}t$ to $\sim10^{-4}$ (the continuous _integrate test test_poisson_process_gamma_integral).

Let $f(x)=c\,x^2$ on $[0,3]$ (and $0$ elsewhere). Use $\int_0^3 f=1$ to find $c$, then compute $E[X]$ and $\operatorname{Var}[X]$. Solution: $\int_0^3 c x^2\,dx=c\cdot\tfrac{27}{3}=9c=1\Rightarrow c=\tfrac19$. Then $E[X]=\tfrac19\int_0^3 x^3\,dx=\tfrac19\cdot\tfrac{81}{4}=\tfrac94=2.25$, and $E[X^2]=\tfrac19\int_0^3 x^4\,dx=\tfrac19\cdot\tfrac{243}{5}=\tfrac{27}{5}=5.4$, so $\operatorname{Var}[X]=5.4-2.25^2=0.3375$. Check: with f = lambda x: x*x/9, pdf_is_normalized(f, 0, 3) $=$ True; expectation_continuous(f, 0, 3) $=2.25$; variance_continuous(f, 0, 3) $=0.3375$.

For the uniform $U(2,8)$ ($f=\tfrac16$ on $[2,8]$), find $F(5)$ and $P(3<X<5)$ two ways: as $\int f$ and as a difference of the cdf. Solution: $F(5)=\int_2^5\tfrac16\,dx=\tfrac{5-2}{6}=\tfrac12$; and $P(3<X<5)=F(5)-F(3)=\tfrac12-\tfrac16=\tfrac13$. Check: with U = lambda x: uniform_pdf(x,2,8), cdf_from_pdf(U, 5, 2) $=0.5$ and prob_between(U, 3, 5) $=0.3333\ldots=$ uniform_cdf(5,2,8)-uniform_cdf(3,2,8).

For the triangular density $f(x)=2x$ on $[0,1]$ (cdf $F(x)=x^2$), show $P(X=\tfrac12)=0$ and that the endpoints are irrelevant, then evaluate $P(0.3<X<0.7)$. Solution: $P(X=\tfrac12)=F(\tfrac12)-F(\tfrac12)=0$, so $P(0.3<X<0.7)=P(0.3\le X\le 0.7)=F(0.7)-F(0.3)=0.49-0.09=0.40$. Check: with g = lambda x: 2*x, prob_between(g, 0.5, 0.5) $=0.0$ and prob_between(g, 0.3, 0.7) $=0.40$.

For $f(x)=2x$ on $[0,1]$ compute $E[X]$ and $\operatorname{Var}[X]$. Solution: $E[X]=\int_0^1 x\cdot2x\,dx=\tfrac23$; $E[X^2]=\int_0^1 x^2\cdot2x\,dx=\tfrac12$; hence $\operatorname{Var}[X]=\tfrac12-\big(\tfrac23\big)^2=\tfrac12-\tfrac49=\tfrac1{18}\approx0.05556$. Check: expectation_continuous(g, 0, 1) $=0.6667$ and variance_continuous(g, 0, 1) $=0.05556\ (=1/18)$.

For $F(x)=x^2$ on $[0,1]$ find the median and the 25th percentile. Solution: solve $F(\pi_p)=p$, i.e. $\pi_p^2=p$, so $\pi_p=\sqrt p$. The median is $\pi_{0.5}=\sqrt{0.5}\approx0.7071$ and the lower quartile is $\pi_{0.25}=\sqrt{0.25}=0.5$. Check: with G = lambda x: x*x, quantile(G, 0.5, 0, 1) $=0.7071$ and quantile(G, 0.25, 0, 1) $=0.5$.

For $U(2,8)$ verify mean $\tfrac{a+b}2$, variance $\tfrac{(b-a)^2}{12}$, and find the 90th percentile. Solution: mean $=\tfrac{2+8}{2}=5$; variance $=\tfrac{(8-2)^2}{12}=\tfrac{36}{12}=3$; percentile $\pi_p=a+p(b-a)$ gives $\pi_{0.9}=2+0.9\cdot6=7.4$. Check: uniform_mean(2,8) $=5.0$, uniform_var(2,8) $=3.0$, uniform_quantile(0.9, 2, 8) $=7.4$, and bisection agrees: quantile(lambda x: uniform_cdf(x,2,8), 0.9, 2, 8) $=7.4$.

Smear a unit point mass at $c=4$ over $[c-\varepsilon,c+\varepsilon]$ with the box density $f_\varepsilon=\tfrac1{2\varepsilon}$. Show $E[X]=c$ for all $\varepsilon$ and $\operatorname{Var}[X]=\tfrac{\varepsilon^2}{3}\to0$ as $\varepsilon\to0$ — the density collapses to $\delta(x-4)$ (~MA-15), recovering a discrete point mass (~ST-05). Solution: by symmetry $E[X]=c=4$; and $\operatorname{Var}[X]=\tfrac1{2\varepsilon}\int_{-\varepsilon}^{\varepsilon}u^2\,du =\tfrac1{2\varepsilon}\cdot\tfrac{2\varepsilon^3}{3}=\tfrac{\varepsilon^2}{3}$. At $\varepsilon=0.1$ this is $0.00333\overline{3}$. Check: with pm = lambda x: point_mass_pdf(x, 4, 0.1), expectation_continuous(pm, 3.9, 4.1) $=4.0$ and variance_continuous(pm, 3.9, 4.1) $=0.003333\ (=0.1^2/3)$.

From $\Gamma(\alpha)=\int_0^\infty t^{\alpha-1}e^{-t}dt$ derive the recursion $\Gamma(\alpha+1)=\alpha\Gamma(\alpha)$ by parts, and use $\Gamma(1)=1$ to get $\Gamma(n)=(n-1)!$. Then, starting from $\Gamma(\tfrac12)=\sqrt\pi$ (the Gaussian integral), evaluate $\Gamma(\tfrac72)=\tfrac52\cdot\tfrac32\cdot\tfrac12\,\Gamma(\tfrac12) =\tfrac{15}{8}\sqrt\pi$. Check: gamma_function(5) $=24.0=4!$; gamma_function(0.5) $=1.772454=\sqrt\pi$; gamma_function(3.5) $=3.323351 =\tfrac{15}{8}\sqrt\pi$; and gamma_lanczos(5.7) matches gamma_function(5.7) to $\sim10^{-13}$.

A device's lifetime is exponential with mean $5$ years, so $\lambda=1/5=0.2$. Find $P(X>10)$ and $P(X\le 10)$. The survival function is $S(x)=e^{-\lambda x}$, so $P(X>10)=e^{-2}=0.13534$ and $P(X\le 10)=1-e^{-2}=0.86466$. Confirm the mean and variance are $1/\lambda=5$ and $1/\lambda^2=25$. Check: exponential_survival(10, 0.2) $=0.135335$; exponential_cdf(10, 0.2) $=0.864665$; exponential_mean(0.2) $=5.0$, exponential_var(0.2) $=25.0$.

The same device has already run $3$ years. Show the chance it lasts $2$ more is the same as for a brand-new one: $P(X>5\mid X>3)=P(X>5)/P(X>3)=e^{-\lambda\cdot 5}/ e^{-\lambda\cdot 3}=e^{-2\lambda}=P(X>2)$. With $\lambda=0.4$ both equal $e^{-0.8}=0.449329$. This is the continuous mirror of the geometric law's memorylessness (~ST-08). Check: exp_memoryless_check(3.0, 2.0, 0.4) $=(0.449329,\,0.449329)$ — the two probabilities are equal.

Customers arrive as a Poisson process at rate $\lambda=1/3$ per minute, so each inter-arrival is exponential with scale $\theta=3$. The wait for the 2nd arrival is gamma with $\alpha=2,\theta=3$. Write its density $f(x)=x\,e^{-x/3}/(\Gamma(2)\,3^2)=\tfrac{x}{9}e^{-x/3}$, and from $M(t)=(1-3t)^{-2}$ read off mean $\alpha\theta=6$ and variance $\alpha\theta^2=18$. Evaluate $f(3)=\tfrac{3}{9}e^{-1}=0.12263$ and the probability the 2nd customer has arrived by $6$ min, $F(6)$. Check: gamma_pdf(3.0, 2.0, 3.0) $=0.122626$; gamma_mean(2.0, 3.0) $=6.0$, gamma_var(2.0, 3.0) $=18.0$; gamma_cdf(6.0, 2.0, 3.0) $=0.593994$.

Let $X_1,\dots,X_\alpha$ be iid exponential with scale $\theta$. Using that mgfs of independent sums multiply, show $M_{\sum X_i}(t)=\big[(1-\theta t)^{-1}\big]^{\alpha}=(1-\theta t)^{-\alpha}$, the gamma$(\alpha,\theta)$ mgf — so by uniqueness (~ST-06) the sum is gamma. For $\alpha=2$ verify by convolution: $f_{X_1+X_2}(x)=\int_0^x\lambda e^{-\lambda y} \lambda e^{-\lambda(x-y)}dy=\lambda^2 x e^{-\lambda x}$. Check: with $\lambda=0.8$ ($\theta=1.25$), convolve_two_exponentials_pdf(2.5, 0.8) $=0.216536=$ gamma_pdf(2.5, 2.0, 1.25); and the mgf identity exponential_mgf(0.2, 1/1.5)**4 $=4.164931=$ gamma_mgf(0.2, 4.0, 1.5).

The chi-square with $r$ degrees of freedom is gamma$(\alpha=r/2,\theta=2)$, so its mean is $\alpha\theta=r$ and its variance $\alpha\theta^2=2r$. Show that for $r=2$ (shape $1$, scale $2$) it collapses to the exponential of mean $2$, i.e. $\lambda=\tfrac12$: $f(x)=\tfrac12 e^{-x/2}$. Check: chi2_mean(10) $=10.0$, chi2_var(10) $=20.0$; chi2_pdf(3.0, 2) $=0.111565=$ exponential_pdf(3.0, 0.5); and chi2_cdf(3.0, 2) $=1-e^{-1.5}=0.776870$.

In a gas at temperature $T$ a molecule's velocity has three independent normal components, so its kinetic energy $E=\tfrac12 m v^2$ has $2E/kT=\sum_{i=1}^3 (v_i/\sigma)^2\sim\chi^2_3$ — a gamma of shape $\tfrac32$ (~SM-06). Hence $E[2E/kT]=3$ gives the equipartition mean energy $\tfrac32 kT$, with variance $\mathrm{Var}(2E/kT)=2\cdot 3=6$. (The same "sum of squared standard normals" construction, $\chi^2_r=\sum_{i=1}^r Z_i^2$, is the launch point of sampling theory ~ST-17; note $P(\chi^2_1\le 1)=0.6827$, the one-$\sigma$ probability.) Check: chi2_mean(3) $=3.0$, chi2_var(3) $=6.0$, chi2_cdf(3.0, 3) $=0.608375$; and chi2_cdf(1.0, 1) $=0.682689$.

Square $I=\int_{-\infty}^\infty e^{-x^2/2}dx$ and pass to polar coordinates to show $I^2=2\pi$, hence $I=\sqrt{2\pi}$. Conclude that $f(x)=\frac{1}{\sigma\sqrt{2\pi}} e^{-(x-\mu)^2/2\sigma^2}$ integrates to one, and that the peak of the standard normal is $f(0)=1/\sqrt{2\pi}$. Check: gaussian_integral_check() $=2.5066283$ $=$ SQRT_2PI; float(standard_normal_pdf(0.0)) $=0.39894228=1/\sqrt{2\pi}$.

Heights are $X\sim N(\mu=3,\sigma=2)$ (in some units). Find the $Z$-score of $x=5$ and verify the density factorizes as $f_X(5)=\varphi(z)/\sigma$. Why does this reduce every normal probability to a question about $N(0,1)$? Check: standardize(5, 3, 2) $=1.0$; float(normal_pdf(5,3,2)) equals float(standard_normal_pdf(1.0))/2.

Show $E[X]=\mu$ by writing $x=(x-\mu)+\mu$ (the first piece integrates to zero by oddness), and $\operatorname{Var}[X]=\sigma^2$ by the substitution $u=(x-\mu)/\sigma$ and $\int u^2e^{-u^2/2}du=\sqrt{2\pi}$. So $N(\mu,\sigma^2)$ is named by its own mean and variance. Check: for $N(3,2^2)$, normal_mean_by_integration(3,2) $=3.000000$ and normal_variance_by_integration(3,2) $=4.000000$.

Complete the square to derive $M(t)=E[e^{tX}]=\exp(\mu t+\tfrac12\sigma^2 t^2)$. Differentiate at $t=0$ to get $M'(0)=\mu$ and $M''(0)=\mu^2+\sigma^2$, and confirm $\operatorname{Var}=M''(0)-M'(0)^2=\sigma^2$. Then show $aX+b\sim N(a\mu+b,a^2\sigma^2)$ from $e^{bt}M(at)$. Check: for $N(3,2^2)$, normal_mgf(0,3,2) $=1.0$, mgf_moment(1,3,2) $\approx3.0$, and mgf_moment(2,3,2) $\approx13.0=\mu^2+\sigma^2$.

Using $P(|Z|\le k)=2\Phi(k)-1$, fill in the probabilities of $\mu\pm1\sigma$, $\mu\pm2\sigma$, $\mu\pm3\sigma$. For IQ scores $X\sim N(100,15^2)$, what fraction of people score between $85$ and $115$ (i.e. within one $\sigma$)? Check: empirical_rule(1), empirical_rule(2), empirical_rule(3) $=(0.6827,0.9545,0.9973)$; float(normal_interval_prob(85,115,100,15)) $=0.6826895$.

SAT-section scores are modeled $X\sim N(500,100^2)$. Standardize $x=650$ to $z=1.5$ and find $P(X>650)=1-\Phi(1.5)$. Separately, what is $P(X>130)$ for the IQ model $N(100,15^2)$ (a two-$\sigma$ event)? Check: 1-float(normal_cdf(650,500,100)) $=0.0668072$ (so $z=1.5\Rightarrow\Phi=0.9331928$); 1-float(normal_cdf(130,100,15)) $=0.0227501$ (the upper-$2\sigma$ tail, half of $1-0.9545$).

The probit values $\Phi^{-1}(0.90)=1.2816$, $\Phi^{-1}(0.95)=1.6449$, $\Phi^{-1}(0.975)=1.9600$ are the multipliers behind one-sided $90\%/95\%$ limits and two-sided $95\%$ intervals. Find the $90$th percentile of the SAT model $N(500,100^2)$ from $x_p=\mu+\sigma\Phi^{-1}(p)$, and check $\Phi^{-1}(\tfrac12)=0$. Check: probit(0.5) $=0$; probit(0.975) $=1.9599640$; normal_quantile(0.90,500,100) $=628.1552=500+100(1.2816)$.

Prove $\Phi(-z)=1-\Phi(z)$ from the evenness of $\varphi$, and use it to show the central $95\%$ of $N(0,1)$ lies in $(-1.96,1.96)$ — i.e. $\Phi(1.96)\approx0.975$, not the round $\Phi(2)=0.9772$ of the empirical rule. Check: float(standard_normal_cdf(1.96)) $=0.9750021$ while float(standard_normal_cdf(2.0)) $=0.9772499$; and float(standard_normal_cdf(-1.96)) $=1-0.9750021=0.0249979$.

Verify $f(x,y)=\tfrac{x+y}{32}$ is a valid joint pmf: $f\ge0$ and the $2\times4$ cell values $x+y$ sum to $32$. Then compute the marginals by summing out the other variable, $f_X(x)=\sum_y f(x,y)$ and $f_Y(y)=\sum_x f(x,y)=\tfrac{3+2y}{32}$, and confirm each is itself a one-variable pmf (sums to $1$). Check: joint_pmf_is_valid(P) $=$ True; marginal_x(P) $=[0.4375,\,0.5625]$ $=[\tfrac{14}{32},\tfrac{18}{32}]$; marginal_y(P) $=[0.15625,0.21875,0.28125,0.34375]$ and marginal_y(P).sum() $=1.0$.

Test $f(x,y)=f_X(x)\,f_Y(y)$. Show it fails at $(1,1)$: $f(1,1)=\tfrac{2}{32}=0.0625$ but $f_X(1)\,f_Y(1)=\tfrac{14}{32}\cdot\tfrac{5}{32} =0.06836$. Conclude $X,Y$ are dependent (a sum $x+y$ never factors as $g(x)h(y)$). Contrast with the outer-product table $P_i=\big[\begin{smallmatrix}.2&.2\\.3&.3\end{smallmatrix}\big] =\mathrm{outer}([.4,.6],[.5,.5])$, which does factor. Check: independent_rv_check(P) $=$ False; independent_rv_check(np.outer([0.4,0.6],[0.5,0.5])) $=$ True.

Using $E[g(X,Y)]=\sum_{x,y}g(x,y)f(x,y)$, compute $E[X]=\tfrac{50}{32}$, $E[Y]=\tfrac{90}{32}$, $E[XY]=\tfrac{140}{32}$. Verify linearity $E[X+Y]=E[X]+E[Y]$ (it needs no independence), and that the gap $E[XY]-E[X]E[Y]=-\tfrac{5}{256}$ is nonzero — the covariance of ~ST-14. Check: expectation_joint(lambda x,y:x, xv,yv,P) $=1.5625$; ...:y... $=2.8125$; ...:x*y... $=4.375$; expectation_joint(lambda x,y:x+y,...) $=4.375\stackrel{?}{=}1.5625+2.8125=4.375$; $E[XY]-E[X]E[Y]=-0.01953$.

Compute $F(x,y)=\sum_{x_i\le x,\,y_j\le y}f(x_i,y_j)$. Show $F(1,2)=f(1,1)+f(1,2)=\tfrac{5}{32}$, the corner $F(2,4)=1$, and that pushing $y\to\max$ recovers the marginal cdf $F_X(1)=F(1,4)=f_X(1)=\tfrac{14}{32}$. Argue $F$ is nondecreasing in each argument. Check: joint_cdf(P,xv,yv,1,2) $=0.15625$; joint_cdf(P,xv,yv,2,4) $=1.0$; joint_cdf(P,xv,yv,1,4) $=$ marginal_cdf_x(P,xv,1) $=0.4375$.

A density has the shape $f\propto xy$ on the unit square. Find the constant $c$ with $\iint c\,xy\,dx\,dy=1$: since $\iint_0^1 xy\,dx\,dy=\tfrac14$, $c=4$, so $f=4xy$. Confirm it is a valid pdf ($f\ge0$, integral $1$). Check: normalize_joint_pdf(lambda x,y:xy, 0,1,0,1) $\approx4.0$; joint_pdf_is_valid(lambda x,y:4x*y, 0,1,0,1) $=$ True.

For $f=4xy$ integrate out $y$: $f_X(x)=\int_0^1 4xy\,dy=2x$, similarly $f_Y(y)=2y$, and note $4xy=(2x)(2y)$ — independent. For $f=x+y$: $f_X(x)=\int_0^1(x+y)\,dy=x+\tfrac12$, which gives $f_X(x)f_Y(y)=(x+\tfrac12)(y+\tfrac12) \ne x+y$ — dependent (they coincide only at the centre, so the check must sample several points). Check: marginal_pdf_x(lambda x,y:4xy, 0.5, 0,1) $=1.0$ $(=2\cdot0.5)$; marginal_pdf_x(lambda x,y:x+y, 0.3, 0,1) $=0.8$ $(=0.3+\tfrac12)$; independent_pdf_check(...4xy...) $=$ True, independent_pdf_check(...x+y...) $=$ False.

For $f=x+y$ compute $E[X]=\iint x(x+y)\,dx\,dy=\tfrac{7}{12}$ and $E[XY]=\iint xy(x+y)=\tfrac13$. For the independent $f=4xy$ show the product rule $E[XY]=E[X]E[Y]=\tfrac23\cdot\tfrac23=\tfrac49$. Finally get the joint cdf of $f=4xy$, $F(x,y)=\iint_0^{x,y}4uv=x^2y^2$, so $F(\tfrac12,\tfrac12)=\tfrac1{16}$ and $F(1,1)=1$. Check: expectation_joint_continuous(lambda x,y:x, lambda x,y:x+y, 0,1,0,1) $\approx0.5833\,(=\tfrac{7}{12})$; ...xy... $\approx0.3333\,(=\tfrac13)$; for $4xy$ expectation_joint_continuous(lambda x,y:xy,...) $\approx0.4444=\tfrac49$; joint_cdf_continuous(lambda x,y:4xy, 0,0, 0.5,0.5) $\approx0.0625$, ...,1,1) $\approx1.0$.

Let $f(x,y)=2$ on the triangle $0<x<y<1$ (area $\tfrac12$, so it integrates to $1$). Because the support is not a rectangle, $X,Y$ cannot be independent. Verify by marginalizing: $f_X(x)=\int_x^1 2\,dy=2(1-x)$ and $f_Y(y)=\int_0^y 2\,dx=2y$, whose product $4y(1-x)\ne2$. (The diagonal makes the 2-D midpoint integral of the indicator slightly inexact, but the marginals — 1-D integrals at fixed $x$ or $y$ — are clean.) Check: with f = lambda x,y: 2.0*(np.asarray(x) < np.asarray(y)), marginal_pdf_x(f, 0.25, 0,1) $=1.5=2(1-0.25)$; marginal_pdf_y(f, 0.6, 0,1) $=1.2=2\cdot0.6$; independent_pdf_check(f, 0,1,0,1) $=$ False.

Using the shortcut $\operatorname{Cov}(X,Y)=E[XY]-E[X]E[Y]$, show for the demo joint that $E[X]=\tfrac{25}{16}=1.5625$, $E[Y]=\tfrac{45}{16}=2.8125$, and $E[XY]=\tfrac{35}{8}=4.375$, hence $$\operatorname{Cov}(X,Y)=\frac{35}{8}-\frac{25}{16}\cdot\frac{45}{16} =\frac{1120-1125}{256}=-\frac{5}{256}.$$ The negative sign says larger $X$ goes (weakly) with smaller $Y$. Check: covariance(demo_joint()) $=-0.01953125=-5/256$; mean_x, mean_y, mean_xy give $1.5625,\,2.8125,\,4.375$.

Rescale the covariance to the dimensionless $\rho=\operatorname{Cov}/(\sigma_X\sigma_Y)$. With $\operatorname{Var}X=\tfrac{63}{256}$ and $\operatorname{Var}Y=\tfrac{295}{256}$, $$\rho=\frac{-5/256}{\sqrt{(63/256)(295/256)}}=\frac{-5}{\sqrt{18585}}\approx-0.03668,$$ comfortably inside $[-1,1]$ — a weak negative linear association. Explain why Cauchy–Schwarz ($\det\Sigma=\sigma_X^2\sigma_Y^2(1-\rho^2)\ge0$) guarantees the bound. Check: correlation(demo_joint()) $\approx-0.036677$ and abs(correlation(...)) <= 1; covariance_matrix(...) has determinant $\sigma_X^2\sigma_Y^2-\operatorname{Cov}^2\approx0.2832>0$.

Prove that $|\rho|=1$ iff $Y=aX+b$ almost surely, with $\rho=\operatorname{sgn}(a)$ (the equality case of $0\le E[(U-tV)^2]$, $U=X-\mu_X$, $V=Y-\mu_Y$). Confirm on two deterministic lines $Y=2X+1$ and $Y=-3X+5$. Check: correlation(linear_dependence_example(2.0, 1.0)) $=1.0$ and correlation(linear_dependence_example(-3.0, 5.0)) $=-1.0$.

First show that independence makes $E[XY]=E[X]E[Y]$, so $\rho=0$. Then take $X\sim U\{-1,0,1\}$ and $Y=X^2$: compute $E[X]=0$, $E[XY]=E[X^3]=0$, so $\operatorname{Cov}=0$ and $\rho=0$ — yet $Y$ is a function of $X$ (maximally dependent). Explain why $\rho$ misses this: the dependence is purely quadratic, invisible to a linear measure. Check: with u = uncorrelated_dependent_example(), covariance(u) $=0.0$, correlation(u) $=0.0$, but is_independent(u) $=$ False; and the conditional means differ, conditional_expectation(u, -1) $=1$, conditional_expectation(u, 0) $=0$, conditional_expectation(u, 1) $=1$.

For the demo joint, condition on $X=1$: $f(y\mid 1)=f(1,y)/f_X(1)=\big(\tfrac{1+y}{32}\big) /\big(\tfrac{14}{32}\big)=\tfrac{1+y}{14}$, a valid pmf ($\sum_{y=1}^4\tfrac{1+y}{14}=1$). Its mean is the regression function at $x=1$: $$E[Y\mid X=1]=\sum_{y=1}^{4}y\,\frac{1+y}{14}=\frac{40}{14}=\frac{20}{7}\approx2.857, \qquad E[Y\mid X=2]=\frac{50}{18}=\frac{25}{9}\approx2.778 .$$ Check: conditional_pmf(demo_joint(), 1) returns probabilities $\big(\tfrac{2}{14},\tfrac{3}{14},\tfrac{4}{14},\tfrac{5}{14}\big)$ summing to $1$; conditional_expectation(demo_joint(), 1) $\approx2.857143=20/7$ and ...(demo_joint(), 2) $\approx2.777778=25/9$.

Verify the tower property $E[E[Y\mid X]]=E[Y]$ for the demo joint by averaging the conditional means over the marginal of $X$: $$\frac{20}{7}\cdot\frac{14}{32}+\frac{25}{9}\cdot\frac{18}{32} =\frac{20\cdot14}{7\cdot32}+\frac{25\cdot18}{9\cdot32}=1.25+1.5625=2.8125=E[Y].$$ Do the same for the $Y=X^2$ counterexample: $1\cdot\tfrac13+0\cdot\tfrac13+1\cdot\tfrac13 =\tfrac23=E[Y]$. Check: law_of_total_expectation(demo_joint()) $=2.8125=$ mean_y(demo_joint()); law_of_total_expectation(uncorrelated_dependent_example()) $\approx0.6667=2/3$.

The best linear predictor has slope $b=\operatorname{Cov}/\operatorname{Var}X= \rho\,\sigma_Y/\sigma_X$ and intercept $a=\mu_Y-b\mu_X$. For the demo joint $b=\tfrac{-5/256}{63/256}=-\tfrac{5}{63}\approx-0.07937$ and $a\approx2.93651$, and the line passes through $(\mu_X,\mu_Y)=(1.5625,2.8125)$. Now take the symmetric law on $\{-1,1\}^2$ with $\rho=0.6$ (regression_to_mean_example()): since $\mu_X=\mu_Y=0$ and $\sigma_X=\sigma_Y=1$, an extreme $X=+1$ predicts $\hat Y=\rho\cdot1=0.6$ — closer to the mean than $X$, the regression-to-the-mean effect forced by $|\rho|\le1$. Check: regression_line(demo_joint()) $\approx(2.93651,\,-0.07937)$; best_linear_predictor(demo_joint(), 1.5625) $\approx2.8125$; for r = regression_to_mean_example(), correlation(r) $=0.6$ and best_linear_predictor(r, 1.0) $=0.6<1$.

For $f(x,y)=x+y$ on the unit square, $f_X(x)=x+\tfrac12$, so $$E[Y\mid X=x]=\frac{\int_0^1 y(x+y)\,dy}{x+\tfrac12}=\frac{\tfrac{x}{2}+\tfrac13}{x+\tfrac12},$$ a nonlinear regression function ($\tfrac23$ at $x=0$, $\tfrac59$ at $x=1$). Show also $E[X]=\tfrac{7}{12}$, $\operatorname{Cov}=-\tfrac{1}{144}$, and $\rho=-\tfrac{1}{11}$, and confirm the tower $\int_0^1 E[Y\mid X{=}x]f_X(x)\,dx= E[Y]=\tfrac{7}{12}$. Check: with f = lambda x,y: x+y, bd = (0,1,0,1): cont_covariance(f, bd) $\approx-0.006944=-1/144$, cont_correlation(f, bd) $\approx-0.0909=-1/11$, cont_conditional_expectation(f, 0.0, bd) $\approx0.6667=2/3$, cont_conditional_expectation(f, 1.0, bd) $\approx0.5556=5/9$, and cont_total_expectation(f, bd) $\approx0.5833=7/12$.

Write the bivariate normal density at $(x,y)=(1.5,3.0)$ two ways: (i) the explicit five-parameter form $f=\frac{1}{2\pi\sigma_X\sigma_Y\sqrt{1-\rho^2}} e^{-z/2(1-\rho^2)}$ with $z=u^2-2\rho uv+v^2$, and (ii) the matrix form $\frac{1}{2\pi\sqrt{\det\Sigma}}e^{-\frac12 w^\top\Sigma^{-1}w}$. Using $u=v=0.5$, get $z=0.25-2(0.6)(0.25)+0.25=0.2$, exponent $-z/2(1-\rho^2)= -0.2/1.28=-0.15625$, prefactor $1/(2\pi\cdot1\cdot2\cdot0.8)=0.099472$, so $f=0.099472\,e^{-0.15625}$. Show the matrix form gives the same number because $w^\top\Sigma^{-1}w=z/(1-\rho^2)=0.3125$ and $\det\Sigma=2.56$. Check: bivariate_normal_pdf(1.5,3.0,1,2,1,2,0.6) $=0.08508277$ equals bivariate_normal_pdf_cov(1.5,3.0,(1,2),covariance_matrix(1,2,0.6)); and mahalanobis_sq(1.5,3.0,1,2,1,2,0.6) $=0.3125$.

Show $\det\Sigma=\sigma_X^2\sigma_Y^2(1-\rho^2)$ and that, by Sylvester's criterion, $\Sigma$ is positive-definite iff its leading minors $\sigma_X^2$ and $\det\Sigma$ are both positive — i.e. iff $|\rho|<1$. Conclude the density exists exactly on $-1<\rho<1$, and that at $\rho=\pm1$ the matrix is singular and the mass collapses onto a line (degenerate case). Verify you can read $\rho$ back off $\Sigma$ as $\Sigma_{xy}/(\sigma_X\sigma_Y)$. Check: with S=covariance_matrix(1,2,0.6), numpy.linalg.det(S) $=2.56$, correlation(S) $=0.6$, is_positive_definite(S) $=$ True; while is_positive_definite(covariance_matrix(1,2,1.0)) $=$ False.

By completing the square $z=(v-\rho u)^2+(1-\rho^2)u^2$ and integrating out $y$, show $f_X(x)=\frac{1}{\sigma_X\sqrt{2\pi}}e^{-(x-\mu_X)^2/2\sigma_X^2}$, i.e. $X\sim N(1,1)$. Evaluate at the mean $x=\mu_X=1$: $f_X(1)=1/\sqrt{2\pi}=0.398942$. (The inner $y$-integral is a normal density and equals $1$.) Check: marginal_pdf_x(1,1,1) $=0.398942$, and the numerical integral _integrate(lambda y: bivariate_normal_pdf(1,y,1,2,1,2,0.6), -14,18) $\approx0.398942$ (this is exactly the assertion in test_marginal_x_is_normal).

Divide the joint by $f_X$ to show $Y\mid X=x\sim N\!\big(\mu_Y+\rho\frac{\sigma_Y} {\sigma_X}(x-\mu_X),\ \sigma_Y^2(1-\rho^2)\big)$. For our numbers the conditional mean is $2+1.2\,(x-1)$ (slope $\rho\sigma_Y/\sigma_X=1.2$) and the conditional variance is $4(1-0.36)=2.56$, independent of $x$. Evaluate at $x=1.5$ and at $x=3.0$. Check: conditional_params(1.5,1,2,1,2,0.6) $=(2.6,\,2.56)$ and conditional_params(3.0,1,2,1,2,0.6) $=(4.4,\,2.56)$ — same variance, mean on the line.

Standardize: if $X$ sits $k$ standard deviations above its mean, show the predicted $Y$ sits only $\rho k$ standard deviations above its mean, since $\frac{E[Y\mid X=x]-\mu_Y}{\sigma_Y}=\rho\,\frac{x-\mu_X}{\sigma_X}$. Take $x=\mu_X+2\sigma_X=3$ (a $+2\sigma$ outlier in $X$): the predicted $Y$ is only $\rho\cdot2=1.2$ standard deviations above $\mu_Y$ — regression to the mean. Also note the conditional variance is the marginal variance shrunk by $1-\rho^2$, so $X$ explains a fraction $\rho^2=0.36$ of $Y$'s variance. Check: conditional_params(3.0,1,2,1,2,0.6)[0] $=4.4$, and $(4.4-\mu_Y)/\sigma_Y= (4.4-2)/2=1.2=\rho\cdot2$.

For a general joint distribution $\mathrm{Cov}=0$ does not imply independence. Show that for the bivariate normal it does: setting $\rho=0$ kills the cross term and $\sqrt{1-\rho^2}=1$, so $f(x,y)=f_X(x)f_Y(y)$. Verify at $(x,y)=(0.5,3.0)$ that the $\rho=0$ joint equals the product of the marginals $N(1,1)$ and $N(2,4)$. Check: bivariate_normal_pdf(0.5,3.0,1,2,1,2,0.0) $=0.06197500$ equals marginal_pdf_x(0.5,1,1)*marginal_pdf_y(3.0,2,2) $=0.06197500$; with $\rho=0.6$ the joint no longer equals that product.

The joint mgf is $M(t_1,t_2)=\exp\!\big(\mu_X t_1+\mu_Y t_2+\tfrac12(\sigma_X^2 t_1^2 +2\rho\sigma_X\sigma_Y t_1t_2+\sigma_Y^2 t_2^2)\big)$. Confirm $M(0,0)=1$, then show $\partial_{t_1}M|_0=\mu_X$, $\partial^2_{t_1}M|_0-\mu_X^2=\sigma_X^2$, and the mixed partial $\partial_{t_1}\partial_{t_2}M|_0-\mu_X\mu_Y=\mathrm{Cov}(X,Y)=\rho\sigma_X \sigma_Y=1.2$. As a numeric anchor evaluate the exponent at $t_1=t_2=0.1$: $0.1+0.2+\tfrac12(0.01+0.024+0.04)=0.337$, so $M=e^{0.337}=1.40074$. Check: mgf(0,0,1,2,1,2,0.6) $=1.0$ and mgf(0.1,0.1,1,2,1,2,0.6) $=1.400739$; the finite-difference moment/covariance identities are asserted in test_mgf_derivatives_give_moments and test_mgf_mixed_partial_is_covariance.

The constant-density curves solve $w^\top\Sigma^{-1}w=c^2$. Diagonalize $\Sigma=V\Lambda V^\top$; the ellipse has semi-axes $c\sqrt{\lambda_1}, c\sqrt{\lambda_2}$ along the eigenvectors, so its area is $\pi c^2\sqrt{\lambda_1 \lambda_2}=\pi c^2\sqrt{\det\Sigma}$. For our $\Sigma$ the eigenvalues are $\lambda\approx\{0.579,4.421\}$; the $c=1$ ellipse has semi-axes $\approx\{0.761,2.103\}$ tilted $\approx70.7^\circ$ from the $x$-axis, with $a\,b=\sqrt{\det\Sigma}=\sqrt{2.56}=1.6$. Confirm every point on it has Mahalanobis distance squared equal to $1$. Check: ellipse_axes(covariance_matrix(1,2,0.6),1.0) $\approx(2.1026,\,0.7610,\,1.2334\,\text{rad})$ with product $2.1026\times0.7610=1.6$; and for points from ellipse_points(1,2,1,2,0.6,c=1), mahalanobis_sq(...) $=1.0$ at each (test_contour_ellipse_level_set).

Let $X\sim U(0,1)$ and $Y=X^2$. Using the distribution-function method, find the inverse image $\{x:x^2\le y\}$ for $0\le y\le1$, integrate $f_X$ over it to get $F_Y(y)=\sqrt y$, then differentiate to $f_Y(y)=\tfrac{1}{2\sqrt y}$. Note the density diverges at $y\to0^+$ yet $\int_0^1\tfrac{1}{2\sqrt y}\,dy=1$. Check: cdf_of_Y(uniform_pdf, lambda x:xx, 0.64, 0.0, 1.0) $=0.8=\sqrt{0.64}$ and pdf_of_Y_cdf_method(uniform_pdf, lambda x:xx, 0.25, 0.0, 1.0) $=1.0=\tfrac{1}{2\sqrt{0.25}}$.

Let $X\sim U(0,1)$ and $Y=-\tfrac1\lambda\ln X$ with $\lambda=2$. Show $g$ is decreasing with $g^{-1}(y)=e^{-\lambda y}$ and $dx/dy=-\lambda e^{-\lambda y}$, so the change-of-variables formula $f_Y(y)=f_X(g^{-1}(y))\,|dx/dy|=\lambda e^{-\lambda y}$ gives $Y\sim\mathrm{Exp}(\lambda)$. Confirm this equals $\tfrac{d}{dy}F_Y$ with $F_Y(y)=1-e^{-\lambda y}$. Check: change_of_variables_1d(uniform_pdf, lambda y: math.exp(-2y), lambda y: -2math.exp(-2*y), 0.5) $=0.73576=$ exp_pdf(0.5, 2.0).

Let $Z\sim N(0,1)$ and $Y=Z^2$. Because $g(z)=z^2$ is two-to-one, the event $\{Z^2\le y\}=\{-\sqrt y\le Z\le\sqrt y\}$; show $F_Y(y)=2\Phi(\sqrt y)-1$ and $f_Y(y)=\dfrac{1}{\sqrt{2\pi y}}\,e^{-y/2}$, i.e. $Y\sim\chi^2_1=\mathrm{Gamma}(\tfrac12,\tfrac12)$ — the first member of the sampling-distribution family of ~ST-17. Check: pdf_of_Y_cdf_method(normal_pdf, lambda x:x*x, 1.0, -12.0, 12.0) $\approx0.2410$, matching gamma_pdf(1.0, 0.5, 0.5) $=0.24197=e^{-1/2}/\sqrt{2\pi}$.

Let $X,Y\sim U(0,1)$ independent and $Z=X+Y$. From the convolution $f_Z(z)=\int f_X(t)f_Y(z-t)\,dt$, show the integrand is $1$ only on the overlap of $[0,1]$ and $[z-1,z]$, whose length gives the triangular density: $f_Z(z)=z$ on $[0,1]$, $2-z$ on $[1,2]$, peaking at $z=1$. Check: convolution(uniform_pdf, uniform_pdf, 0.7, 0.0, 1.0) $=0.7$ and convolution(uniform_pdf, uniform_pdf, 1.5, 0.0, 1.0) $=0.5$, equal to triangular_pdf(0.7) and triangular_pdf(1.5).

Let $X_1,\dots,X_n\sim\mathrm{Exp}(\lambda)$ independent. Convolve two: $f_{X_1+X_2}(z)=\int_0^z\lambda e^{-\lambda t}\lambda e^{-\lambda(z-t)}\,dt =\lambda^2 z\,e^{-\lambda z}=\mathrm{Gamma}(2,\lambda)$. By induction the sum of $n$ is $\mathrm{Gamma}(n,\lambda)=\tfrac{\lambda^n}{(n-1)!}z^{n-1}e^{-\lambda z}$ — the Erlang waiting time for the $n$-th Poisson event. Check: with $\lambda=1$, convolution(exp_pdf, exp_pdf, 2.0, 0.0, 2.0) $=0.27067=$ gamma_pdf(2.0, 2, 1); convolving that Gamma(2) once more, convolution(lambda z: gamma_pdf(z,2,1), exp_pdf, 3.0, 0.0, 3.0) $=0.22404=$ gamma_pdf(3.0, 3, 1).

Let $X,Y\sim\mathrm{Exp}(1)$ independent, $U=X+Y$, $V=\dfrac{X}{X+Y}$. Invert to $x=uv,\ y=u(1-v)$, compute the Jacobian $J=\det\!\begin{pmatrix}v&u\\1-v&-u\end{pmatrix}=-u$, and show $f_{UV}(u,v)=e^{-uv}e^{-u(1-v)}\,|{-u}|=u\,e^{-u}$, which factorizes into $\mathrm{Gamma}(2,1)$ for $U$ and $U(0,1)$ for $V$ — so the total and the ratio are independent. Check: with fXY = lambda x,y: float(exp_pdf(x,1))float(exp_pdf(y,1)) and inv = lambda u,v: (uv, u(1-v)), jacobian_det_2d(inv, 2.0, 0.7) $=-2.0$ and jacobian_transform_2d(fXY, inv, 2.0, 0.7) $=0.27067=2e^{-2}=$ gamma_pdf(2.0,2,1)uniform_pdf(0.7).

(a) Show that for any continuous $X$ with cdf $F_X$, the variable $U=F_X(X)$ has density $f_U(u)=f_X(F^{-1}(u))\,|dF^{-1}/du|=1$ on $(0,1)$, i.e. $U\sim U(0,1)$. (b) Run it backwards: if $U\sim U(0,1)$ then $X=F_X^{-1}(U)$ has cdf $F_X$, so drawing $X=-\tfrac1\lambda\ln(1-U)$ samples $\mathrm{Exp}(\lambda)$. Check: (a) probability_integral_transform_pdf(lambda x: exp_pdf(x,2.0), lambda u: exp_quantile(u,2.0), 0.6) $=1.0$ and the same for the linear density, probability_integral_transform_pdf(linear_pdf, linear_quantile, 0.6) $=1.0$. (b) exp_cdf(exp_quantile(0.8, 2.0), 2.0) $=0.8$ and the median exp_quantile(0.5, 2.0) $=0.34657=\ln 2/2$.

Let $X_1\sim N(1,4)$, $X_2\sim N(-0.5,2.25)$, $X_3\sim N(2,1)$ be independent. Use $M_{\sum X_i}(t)=\prod_iM_{X_i}(t)$ to show $S=X_1+X_2+X_3$ is normal, and find its mean and variance. Why does the product of three exponential-quadratics stay an exponential-quadratic? Check: sum_of_normals_via_mgf([1,-0.5,2],[4,2.25,1]) $=(2.5,7.25)$; product_of_normal_mgfs(0.3,[1,-0.5,2],[4,2.25,1]) $=2.9336577$ equals normal_mgf(0.3, 2.5, math.sqrt(7.25)).

Heights of men and women are modeled $X_1\sim N(5,4)$, $X_2\sim N(3,9)$ (independent, some units). Find the distribution of the difference $D=X_1-X_2$. Students often expect $\operatorname{Var}(D)=4-9<0$ — explain why the correct answer adds the variances. Check: linear_combination_of_normals([1,-1],[5,3],[4,9]) $=(2.0, 13.0)$.

Independent $V_1\sim\chi^2_3$, $V_2\sim\chi^2_5$, $V_3\sim\chi^2_2$. Show $V_1+V_2+V_3\sim\chi^2_{10}$ by multiplying mgfs, and explain how this is the same fact as "a sum of $n$ independent squared standard normals is $\chi^2_n$." Check: sum_of_chisquares([3,5,2]) $=10.0$; product_of_chi2_mgfs(0.2,[3,5,2]) equals chi2_mgf(0.2,10).

For $Z\sim N(0,1)$, derive the density of $W=Z^2$ by the change of variables $w=z^2$ (two branches), and identify it as $\chi^2_1$. Confirm the special case $\chi^2_2$ is the exponential of mean $2$. Check: square_of_standard_normal_pdf(2.0) $=0.103777=$ chi2_pdf(2.0,1); and chi2_pdf(2.0,2) $=0.18394=\tfrac12e^{-1}$.

A normal population has $\mu=100$, $\sigma^2=225$ (so $\sigma=15$). A random sample of $n=25$ is drawn. Find the distribution of $\bar X$ and its standard error, and show the standardized mean $Z=(\bar X-\mu)/(\sigma/\sqrt n)$ is exactly $N(0,1)$. Check: sampling_dist_of_mean(100,225,25) $=(100.0, 9.0)$ (so SE $=\sqrt9=3$); standardized_mean_mgf(0.7,100,225,25) $=1.2776213=e^{0.7^2/2}$.

Define $T=Z/\sqrt{V/r}$ with $Z\sim N(0,1)$, $V\sim\chi^2_r$ independent. State the density, give $E[T]$ and $\operatorname{Var}[T]$, and explain why estimating $\sigma$ by $S$ turns the standardized mean into $t_{n-1}$. Why is the $t_1$ (Cauchy) so pathological? Check: students_t_pdf(0.0,3) $=0.367553$; students_t_var(5) $=1.6667=5/3$; students_t_var(2) $=\infty$; the $t_1$ peak students_t_pdf(0.0,1) $=0.31831=1/\pi$.

Show $f_T(t)\to\varphi(t)$ as $r\to\infty$ by taking the limits of its two pieces, and confirm numerically that $t_{2000}$ is essentially standard normal. Why does $\operatorname{Var}(t_r)=r/(r-2)\to1$? Check: students_t_pdf(1.0,2000) $=0.241910$ vs standard_normal_pdf(1.0) $=0.241971$; students_t_var(1000) $<$ students_t_var(10).

Define $F=(U/r_1)/(V/r_2)$ with $U\sim\chi^2_{r_1}$, $V\sim\chi^2_{r_2}$ independent. State its density and mean, then prove that the square of a $t_r$ variable is $F_{1,r}$. Check: f_mean(10) $=1.25=10/8$; f_pdf(1.0,5,10) $=0.49548$; and the squared-$t$ density t_squared_pdf(4.0,6) $=0.032018$ equals f_pdf(4.0,1,6).

A fair die $X$ is uniform on $\{1,\dots,6\}$ — flat, nothing like a bell. Let $S_n=X_1+\dots+X_n$. State the CLT for the standardized sum and explain why $\max_k|F_{S_n}(k)-\Phi(\tfrac{k+0.5-\mu_n}{\sigma_n})|$ should fall to zero with $n$, where $\mu_n=3.5n$, $\sigma_n=\sqrt{(35/12)\,n}$. Check: clt_cdf_max_error(die_pmf(),1) $=0.05424$ (one die: visibly non-normal), while clt_cdf_max_error(die_pmf(),16) $=0.00163$ and clt_cdf_max_error(die_pmf(),64) $=0.00041$ (the sum is the bell to four places).

For i.i.d. $X_i$ with mean $\mu$ and variance $\sigma^2$, show $E[S_n]=n\mu$ and $\operatorname{Var}(S_n)=n\sigma^2$ (so $\operatorname{SD}=\sigma\sqrt n$), and that $\bar X=S_n/n$ has mean $\mu$ and variance $\sigma^2/n$ — the standard error. Why does the pmf of $S_n$ equal the $n$-fold convolution of the parent? Check: for five dice, nfold_pmf(die_pmf(),5) sums to $1.0$, with pmf_mean $=17.5=5\cdot3.5$ and pmf_var $=14.5833=5\cdot\tfrac{35}{12}$; and convolve_pmf(die_pmf(),die_pmf()) equals nfold_pmf(die_pmf(),2).

Let $Y_i=(X_i-\mu)/\sigma$, so $E[Y]=0$, $E[Y^2]=1$, and $Z_n=\frac1{\sqrt n}\sum_iY_i$. Using $M_{Z_n}(t)=[M_Y(t/\sqrt n)]^n$ and the expansion $M_Y(s)=1+\tfrac12 s^2+o(s^2)$, show $M_{Z_n}(t)\to e^{t^2/2}$, the mgf of $N(0,1)$, and conclude $Z_n\to N(0,1)$ by uniqueness. Which features of the parent survive the limit? Check: the limit target is the standard normal: float(standard_normal_cdf(0.0)) $=0.5$, float(standard_normal_cdf(1.96)) $=0.975002$, and (continuous check) the integral of standard_normal_pdf up to $1$ equals float(standard_normal_cdf(1.0)) $=0.841345$.

State the Berry–Esseen bound $\sup_x|F_{S_n}(x)-\Phi(\tfrac{x-n\mu}{\sigma\sqrt n})| \le C\rho/(\sigma^3\sqrt n)$ with $\rho=E|X-\mu|^3$, and read off the rate. For the die, compute $\rho$ and verify the bound at $n=1$. Check: pmf_third_abs_moment(die_pmf()) $=6.375$; berry_esseen_bound(die_pmf(),1) $=0.9797$ and kolmogorov_cdf_error(die_pmf(),1) $=0.1434\le0.9797$; the bound at $n=4$ is $0.4899=0.9797/\sqrt4$.

For $X\sim\mathrm{Bin}(20,\tfrac12)$ ($\mu=10$, $\sigma=\sqrt5$), approximate $P(X\le13)$ with and without the continuity correction $P(X\le k)\approx\Phi(\tfrac{k+0.5-\mu}{\sigma})$, and compare to the exact value. Why is the correction needed? Check: binom_cdf(13,20,0.5) $=0.942341$ (exact), normal_approx_binomial(13,20,0.5,False) $=0.910144$ (no correction), normal_approx_binomial(13,20,0.5) $=0.941238$ (corrected).

Two refinements: (a) the local de Moivre–Laplace theorem $P(X=k)\approx \tfrac1\sigma\varphi(\tfrac{k-\mu}\sigma)$ approximates a single bar; (b) an interval uses $P(a\le X\le b)\approx\Phi(\tfrac{b+0.5-\mu}\sigma)-\Phi(\tfrac{a-0.5-\mu}\sigma)$. For $\mathrm{Bin}(20,\tfrac12)$ estimate $P(X=10)$ and $P(8\le X\le12)$. Check: de_moivre_laplace_pmf(10,20,0.5) $=0.178412$ vs exact binom_pmf(10,20,0.5) $=0.176197$; normal_approx_binomial_interval(8,12,20,0.5) $=0.736448$ vs exact $0.736824$ (binom_cdf(12,20,0.5)-binom_cdf(7,20,0.5)).

Since $\mathrm{Poisson}(\lambda)$ is a sum of $m$ independent $\mathrm{Poisson}(\lambda/m)$ pieces, the CLT gives $\mathrm{Poisson}(\lambda)\approx N(\lambda,\lambda)$ as $\lambda\to\infty$. Argue the max cdf error decays like $1/\sqrt\lambda$, and state the binomial rule of thumb $np\ge5,\ n(1-p)\ge5$. Check: poisson_approx_max_error(4.0) $=0.03218$, poisson_approx_max_error(16.0) $=0.01648$, poisson_approx_max_error(64.0) $=0.00829$ (each $\times4$ in $\lambda$ halves the error); normal_approx_applicable(20,0.5) $=$ True but normal_approx_applicable(8,0.5) $=$ False ($np=4<5$).