The secretary problem: also called the matching or marriage problem.
You will see \(T\) items, one at a time, from a distribution of types \(F(x)\).
You want the maximum. After you pass, you cannot go back.
\(t\) is the slot you are in: \(t=1,\ldots,T\).
Strategy: skip a prefix, then take the next candidate who is better than everyone you have seen.
For large \(T\) the prefix is a fraction \(1/e\) of the list: \[
\frac{1}{e} \approx 0.368.
\]
A Simple Example
Say I have a budget to consider \(T=10\) candidates, and they have the following scores (which I do not know)
1 7 4 5 3 8 10 3 0 11
I decide to screen 3 candidates before making a decision.
My best candidate from the screening pool is 7
Next is 5 –> dump
Next is 3 –> dump
Next is 8 –> accept
Probability and Decision - Strategy
What’s your best strategy?
Turns out its insensitive to the choice of distribution.
Although there is the random sample i.i.d. assumption lurking.
You’ll not doubt get married between \(18\) and \(60\). Waiting\(\frac{1}{e}\) along this sequence gets you to age \(33\). Then, pick the next best person!
If you take the best so far at \(t\)
\(T\) is the length of the list. Ten candidates, \(T=10\).
\(t\) is the slot you are looking at. In the example you accepted at \(t=6\).
Best so far at \(t\) means better than the first \(t-1\). Candidate 8, after a screen whose best was 7.
The overall best is equally likely to sit in any of the \(T\) positions.
Overall best implies best of the first \(t\), so \[
P(\text{overall best}\mid \text{best so far at }t)
= \frac{P(\text{overall best})}{P(\text{best so far at }t)}
= \frac{1/T}{1/t}
= \frac{t}{T}.
\]
Take that candidate and you win with probability \(t/T\). In the example, \(6/10\).
If \(t\) is the overall best, it is automatically the best of the first \(t\). The earlier ranks do not change the ratio.
Being the best of the first \(t\) has probability \(1/t\).
Success under the \(r\)-rule
Reject the first \(r\). Then take the first later candidate who is the best so far.
You win only if you take the overall best. Put that person at slot \(a\). Probability \(1/T\).
You never take them if \(a\le r\): still screening.
If \(a>r\) you still miss them whenever someone in \(r+1,\ldots,a-1\) was the best so far. That happens when the best of the first \(a-1\) sits after \(r\). You would have taken that person and gone home.
Success under the \(r\)-rule
So you need the best of the first \(a-1\) to sit in the skipped prefix. Those \(a-1\) ranks are in random order, so their max is equally likely in any of those \(a-1\) slots: \[
\frac{r}{a-1}.
\]
In the example, \(r=3\). If the overall best is at \(a=6\), the first five are \(1,7,4,5,3\). Their max is \(7\), in the prefix. Chance \(3/5\).
Denominator \(a-1\), not \(a\): slot \(a\) is already the overall best. The comparison is among the others. \[
P(\text{success})
= \frac{1}{T}\sum_{a=r+1}^{T}\frac{r}{a-1}
= \frac{r}{T}\sum_{b=r}^{T-1}\frac{1}{b}.
\]
Need \(r\ge 1\). If \(r=0\) you take the first candidate and win with probability \(1/T\).
If nobody after \(r\) is the best so far, you take the last one. That path never wins.
Then \(r=T/e\)
For large \(T\) the harmonic sum is a log: \[
\frac{r}{T}\sum_{b=r}^{T-1}\frac{1}{b}
\approx \frac{r}{T}\log\frac{T}{r}.
\]
Differentiate in \(r\) and set the derivative to zero: \[
\frac{1}{T}\log\frac{T}{r}-\frac{1}{T}=0
\qquad\Rightarrow\qquad
r=\frac{T}{e}.
\]
Plug back in: \[
P(\text{success})\approx\frac{1}{e}\log e=\frac{1}{e}.
\]
About \(37\) percent for large \(T\), approached from above.
This \(r\) is best among skip-then-take-the-next-best-so-far rules. The Bellman step says the optimum is of that form.
The distribution of the scores never entered. Only the random order of the ranks.
Bellman check
\(p^*(x_{t-1}=0)\): chance of eventual success if you have not hired by \(t-1\).
At \(t\) the candidate is the best so far with probability \(1/t\). Take them if \(t/T\) beats the value of continuing.
Continuing after you pass on a best-so-far has the same value as continuing after a weaker candidate: the remaining permutation is still uniform. $$ p^(x_{t-1}=0) = ,p^(x_t=0)
(t/T,,p^*(x_t=0)). $$
Once \(t/T\) wins the max it keeps winning: the hire payoff rises, continuation falls.
The recursion telescopes from the end, \(p^*(x_T=0)=0\), to \[
p^*(x_{t-1}=0)=\frac{t-1}{T}\sum_{\tau=t-1}^{T-1}\frac{1}{\tau},
\]
That is the success probability of the rule that starts accepting at \(t-1\). Same \(r\).
Bellman: after the first decision the rest of the policy is still optimal.
\(1/e\) in a sample of \(1000\)
The peak is near \(T/e\) and height \(1/e\) (Monte Carlo, so the top wiggles).
Average quality, not shown, peaks much earlier.
Maximizing the chance of the absolute best is not the same as maximizing expected rank.
Decision Trees
Decision Tree
Catastrophe Modeling
You live in a house that is somewhat prone to mud slides.
Each rainy season there is a 1% chance of a mud slide occurring.
You estimate that a mud slide would do $1 million in damage.
You have the option of building a retaining wall that would help reduce the chance of a devastating mud slide. The wall costs $40,000 to build, and if the slide occurs, the wall will hold with a 95% probability.
You also have the option of a Geologist’s opinion.
Should you build the wall?Should you use the Geologist’s Test and Bayes Rule?
Decision Tree Setup
Let’s formally solve this as follows:
Build a decision tree.
The tree will list the probabilities at each node. It will also list any costs there are you going down a particular branch.
Finally, it will list the expected cost of going down each branch, so we can see which one has the better risk/reward characteristics.
There’s also the possibility of a further test to see if the wall will hold.
Decision Tree
Mudslide Decision Tree
Testing Option
Let’s include the testing option
You also have the option of having a test done to determine whether or not a slide will occur in your location.
The test costs $3000 and has the following accuracies. \[P(T \mid \text{Slide}) = 0.90 \quad \text{and} \quad P(\bar{T} \mid \text{No Slide}) = 0.85\]
If you choose the test, then should you build the wall?
Bayes Rule for Testing
Bayes Rule is as follows:
The initial prior probabilities are \[P(\text{Slide}) = 0.01 \quad \text{and} \quad P(\text{No Slide}) = 0.99\]
We’ll use this to find our optimal course of action.
Bayes Probabilities - Positive Test
The Bayes probability given a positive test is \[\begin{align*}
P(\text{Slide} \mid T) &= \frac{P(T \mid \text{Slide}) P(\text{Slide})}{P(T)} \\
&= \frac{0.90 \times 0.01}{0.1575} = 0.0571
\end{align*}\]
Bayes Probabilities - Negative Test
The Bayes probability given a negative test is \[\begin{align*}
P(\text{Slide} \mid \bar{T}) &= \frac{P(\bar{T} \mid \text{Slide}) P(\text{Slide})}{P(\bar{T})} \\
&= \frac{0.1 \times 0.01}{0.8425} \\
&=0.001187
\end{align*}\]
Compare this to the initial base rate of a \(1\)% chance of having a mud slide.
Probability Lose Everything - Without Testing
You build the wall without testing, what’s the probability that you lose everything?
With the given situation, there is one path (or sequence of events and decisions) that leads to losing everything:
You choose the test, what’s the probability that you’ll lose everything?
There are two paths that lead to losing everything:
First Path: There are three things that have to happen to lose everything Test +ve (\(P=0.1575\)), Build, Slide (\(P= 0.0571\)), Doesn’t Hold (\(P=0.05\))
Second Path: Now you lose everything if Test -ve (\(P=0.8425\)), Don’t Build, Slide given negative (\(P=0.001187\))
Conditional Probabilities
For the first path \[P(\text{first path}) = 0.1575 \times 0.0571 \times 0.05 = 0.00045\]
For the second path \[P(\text{second path}) = 0.8425 \times 0.001187 = 0.00101\]
Hence putting it all together \[P(\text{losing everything} \mid \text{testing}) = 0.00045 + 0.00101 = 0.00146\]
Risk and Reward
Risk-Return Trade-off
Choice
Expected Cost
Risk
P
Don’t Build
$10,000
0.01
1 in 100
Build w/o testing
$40,500
0.0005
1 in 2000
Test
$10,760
0.00146
1 in 700
Expected Cost: Fee + Build + Loss
Expected cost: \(3 + 40\times 0.1575+ 1000\times 0.00146\) or $10,760
What do you choose?
Estimation
Three situations
Map for the hour
No data
No sample. You post a price.
Coin, a bet, coherence.
Secretary. Mudslide prior.
Good data
n/N you can trust.
More of the same sample, tighter around the right number.
Alpha particles: the worked example.
Bad data
Looks like 2. Most dangerous.
Tiny interval, wrong sample.
Literary Digest, WHI, Spirit.
The trap is mixing 2 and 3. The machine sounds like 2.
1. No data
No sample from the process you will decide for.
Probability is a language for uncertainty. Usually lack of information, not a property of the coin.
Operational answer: a price you would bet at.
The number has to cohere. Otherwise someone books you.
How many coins in my pocket?
Probability is a property of what you know.
Probability as a bet
Ramsey and de Finetti: your probability is the odds at which you will bet.
de Finetti: “Probability does not exist.” Uncertainty is personal. It is not mass.
Different news, a trade, a price.
2016: bookmakers posting prices on the same event. Those prices are probabilities.
This course uses both languages. Frequency is how we talk about data. The organizing view is belief, updated by evidence.
Presidential odds, 2016
The coin is risk
Most people agree on \(P(\text{Heads})=1/2\). We elicit that number as a willingness to play.
Will you pay $1 to receive $2 if the coin lands tails, and $0 if it lands heads? If you are indifferent, your \(P(\text{Tails})\) is \(1/2\).
Probability and payoff arrive together. You cannot separate them.
Litowitz asks the same question, then raises the payoff. Pay $1 for a fair coin that pays $2 on heads. Most people pause. Make it $3 on heads. Almost everyone says yes.
We know the possibilities and we know the odds, fifty-fifty. That is risk.
The mudslide tree was the same object: \(0.01\) from historical seasons, \(0.95\) treated as known.
The odds are not posted
Often you will decide on something that does not look like a coin. A career, which model to ship, whether the merger closes.
You know the possibilities and not the probabilities. That is ambiguity, not risk: Ellsberg’s second urn.
Knight called the same contrast uncertainty. That is the word Litowitz used.
Sometimes you do not even know the possibilities. Then you are past Ellsberg.
The tools for the coin do not transfer.
Three bad moves: freeze, pretend the odds are knowable, or hand the decision to a machine that sounds sure.
The machine is superb at the coin. The rest stay yours.
The move that still works: a strong provisional view, name what would change your mind, run the cheapest experiment, update.
Stay the author of the decision. That is Bayes as a habit, not a formula.
If you set \(P(A)=0.7\) and \(P(\bar A)=0.2\), someone can lock in a loss against you.
Pay $0.70 on \(A\) and $0.20 on \(\bar A\). You pay $0.90 and receive $1 either way.
Coherence means no Dutch book exists. In finance that is no-arbitrage.
The leak is \(P(A)+P(\bar A)=0.9\neq 1\).
Kolmogorov’s axioms
Coherence forces the same rules Kolmogorov wrote down.
Non-negativity: \(P(A)\ge 0\)
Normalization: \(P(\Omega)=1\)
Additivity: if \(A\) and \(B\) cannot both happen, \(P(A\text{ or }B)=P(A)+P(B)\)
The complement rule \(P(A)+P(\bar A)=1\) is 2 and 3 together. Violate any of them and someone can lock in a loss.
How do we estimate a probability?
Assigning the number is the hard part.
No data: a price you would bet at.
Good data: a long-run frequency from the process you will decide for.
Bad data: same plot, wrong process. Looks trustworthy. That is the dangerous one.
2. Good data
A sample from the process you will decide for.
Then you estimate: n/N, a mean, a Poisson \(\lambda\).
Trust in the usual sense: more good data, smaller variance around the right number.
Frequency, then a sampling distribution. Alpha particles after the warning.
Frequency
Event \(A\) is heads.
Each grey path is one sequence of tosses. The blue path is one of them. The red line is \(1/2\).
As \(N\) grows, \(n/N\) settles near \(1/2\).
That is the law of large numbers.
von Mises wrote the limit as the definition: \[
P(\text{Heads}) = \lim_{N \to \infty} \frac{n}{N}.
\]
Fitting a probability law
Fit a Poisson or a Normal to observed data.
Most laws depend on a few parameters: \(\lambda\), or \(\mu\) and \(\sigma\).
Estimate those from the sample. Then the distribution is specified.
Then ask whether the fit is any good. Later.
3. Bad data
Same shape as 2: rows, a histogram, a tight number.
Blue is the sample. Orange is the law you fit. They already miss.
The red target is the population you will decide for.
Most dangerous because it looks like good data. Tiny interval, wrong sample.
Bias, not variance. We name this after sampling distributions: Literary Digest, WHI.
Wrong law is the same trap: a tight Normal on the sample, a crash the model called impossible.
The 1987 crash
October 1987: S&P 500 dropped \(-21.76\%\).
Prior monthly returns: mean \(1.2\%\), volatility \(4.3\%\).
The model fitted to earlier months: \(X \sim N(0.012, 0.043^2)\).
Standardize: \[
Z = \frac{-0.2176 - 0.012}{0.043} = -5.34.
\]
A 5-sigma month. Under that Normal, \(P(Z < -5.34) \approx 4.6 \times 10^{-8}\).
The 1987 crash: the model
We assumed returns are Normal. The middle of the histogram looks like 2.
The tails are not. Extreme months happen. The crash is one of them.
Same picture as the orange curve: a law fitted to the sample, a target the law cannot see.
Week 1 computed the z-score. The lesson here is the assumption.
S&P 500 returns vs Normal
The limit is not always a method
Same plot as 2: \(n/N\) settling toward a number. That is what makes 3 dangerous.
The von Mises limit needs an infinite sequence of repeatable trials. You never observe it.
Ville pointed this out. Bernoulli’s workable version is the finite ratio \(n/N\) in the experiment you actually ran.
It requires the same conditions. Change the coin, and the number changes.
Blue is fair. Orange is \(P(\text{Heads})=0.7\). Same \(N\). Two different samples. Two different limits.
Data-driven AI
Modern AI is data-driven. The rest of the course is this approach: estimate from samples, then predict.
A model trained on a limited sample is a frequency from that sample. The interval can be tiny around the wrong number.
Very common in modern AI: confident when it should not be. The machine sounds sure. The sample was not the population.
Data is a currency
Self-driving is sold as a data race: the fleet that logs the most miles trains the better driver. Tesla’s pitch. Waymo’s counter: unsupervised miles are not the same as a human still in the seat.
Data sells. Spirit Airlines, in bankruptcy in 2026, auctioned an operational archive. Google’s bid was $10 million: emails, chats, operations, for AI.
The named passenger list was carved out of that Google deal. The estate kept the right to sell it to travel buyers.
Whenever we have a sample we should estimate from it.
First a pattern with no named law: correlation. Two series, one number \(r\).
Then a law with a named parameter. The alpha counts: we estimate \(\lambda\) in the open.
Good data if the sample is from the process you will decide for. Bad data is the standing risk: same formulas, wrong population.
The mudslide \(0.01\) was a frequency from historical data, then treated as known.
Working mechanism after the pattern: frequency. \(n/N\), likelihood, sampling distributions.
Unsupervised patterns
Why Pattern Matching?
Traditional approach: estimating parameters when functional form is known
Linear relationships (Earth-shape example)
Multiplicative relationships (Newton’s \(F=ma\))
Modern challenge: relationship itself is unknown
Human behavior, natural language
Complex, non-mathematical relationships
Pattern matching: identifying recurring sequences, relationships, or structures
Like finding puzzle pieces that complete a picture
Early Example: Mammography
Albert Solomon (1913)
German surgeon, pioneer of X-ray imaging
“Contributions to the Pathology and Clinic of Breast Cancers”
First to recognize calcium deposits-cancer association
Challenges
Subjective visual assessment
Prone to errors in dense breasts
Subtle abnormalities easily missed
Pattern recognition
Dense tissue areas (masses)
Small white calcium deposits (microcalcifications)
Impact
Foundation for mammography screening
Saved countless lives
Molecular mechanisms still not fully understood
Feynman on Pattern Matching
Chess Analogy
“Let’s say a chess game. And you don’t know the rules… you try to figure out what the rules are… You might discover that when there’s only one bishop around on the board, that the bishop maintains its color. Later on you might discover the law for the bishop is that it moves on a diagonal…”
Experts develop pattern recognition abilities
Chess players don’t calculate every move
They recognize advantageous patterns
Skill transfers to other domains
Machine vs Human Intelligence
Feynman on AI
“When we make a bird, the airplanes don’t fly like a bird… It’s different, right? So, there’s no question that the later machines are not going to think like people think… they’re not going to do arithmetic the same way as we do arithmetic, but they’ll do it better.”
Evolution in Chess:
1996: Computers stronger than grandmasters
AlphaZero: Self-learning through competition
Neural networks learn features vs human heuristics
Preference for initiative, speed, momentum over material
Correlations: Simplest Pattern Matching
Correlation coefficient: measures strength of linear relationship
Plot the distribution of the simulated electoral votes for Obama.
Code
library(lattice)# Lattice Graphdensityplot(sim.EV, plot.points ="rug", xlab ="Electoral Votes for Obama", panel =function(x, ...) {panel.densityplot(x, ...)panel.abline(v =270)panel.text(x =285, y =0.01, "270 EV to Win")panel.abline(v =332)panel.text(x =347, y =0.01, "Actual Obama")}, main ="Electoral College Results Probability")
Beta distribution
A law on a single probability, \(\theta \in (0,1)\). Coin, poll share, conversion rate. \[
f(\theta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\,\theta^{\alpha-1}(1-\theta)^{\beta-1}.
\]
\(\alpha\) and \(\beta\) are positive. Mean \(\alpha/(\alpha+\beta)\).
\(\mathrm{Beta}(1,1)\) is uniform. Larger \(\alpha+\beta\): tighter around the mean.
Two shares that sum to 1: Obama vs not Obama is a Beta. Three shares need Dirichlet.
Dirichlet Distribution
Same idea for \(k\) shares that sum to 1. Poll: Republican, Democrat, other.
\(k=2\) is Beta: \(\mathrm{Dirichlet}(\alpha,\beta)=\mathrm{Beta}(\alpha,\beta)\).
Parameters \(\alpha=(\alpha_1,\ldots,\alpha_k)\), each \(\alpha_i>0\).
Support: \(x_i\ge 0\) and \(\sum_i x_i=1\). \[
f(x\mid\alpha)=\frac{\Gamma(\sum_i\alpha_i)}{\prod_i\Gamma(\alpha_i)}\prod_{i=1}^k x_i^{\alpha_i-1}.
\]
Dirichlet Distribution
Application in Election Forecasting:
Models uncertainty in polling percentages
Three components: Republican, Democrat, Undecided
Parameters reflect poll sample sizes and percentages
Generates probability distributions for each candidate’s support
Why Use Dirichlet for Polling?
Naturally handles the constraint that percentages sum to 100%
Incorporates polling uncertainty
Allows for Monte Carlo simulation of election outcomes
Conjugate prior for multinomial likelihood (Bayesian convenience)
2008 Election Example
Results of recent state polls in the 2008 United States Presidential Election between Barack Obama and John McCain.
Code
# Source: LearnBayes library#| fig-height: 6election.2008=read.csv("data/election2008.csv")attach(election.2008)## Dirichlet simulationprob.Obama =function(j) { p=rdirichlet(5000,500*c(M.pct[j],O.pct[j],100-M.pct[j]-O.pct[j])/100+1)mean(p[,2]>p[,1]) }## sapply function to compute Obama win prob for all statesObama.win.probs=sapply(1:51,prob.Obama)## sim.EV functionsim.election =function() { winner =rbinom(51,1,Obama.win.probs)sum(EV*winner) }sim.EV =replicate(1000,sim.election())## histogram of simulated electionhist(sim.EV,min(sim.EV):max(sim.EV),col="blue",prob=T)abline(v=365,lwd=3) # Obama received 365 votestext(375,30,"Actual \n Obama \n total")
Estimation
Case Study: Emissions of Alpha Particles
Berkson (1966) conducted a careful analysis of data obtained from the National Bureau of Standards.
The source of the alpha particles was americium 241.
n is the number of counts that were observed in 1207 intervals, each of length 10 sec.
In 18 of the 1207 intervals, there were 0, 1, or 2 counts; in 28 of the intervals there were 3 counts, etc.
lmb =sum(d$Observed*d$n)/sum(d$Observed)/10# mean emission ratelmb
[1] 0.8354598
Case Study: Emissions of Alpha Particles
To fit a Poisson distribution, we need to estimate its single parameter, λ.
\[ P(X=k) = \frac{\lambda^k e^{-\lambda}}{k!} \]
A natural choice for the estimate, \(\hat{\lambda}\), is the sample mean of the observed counts.
Estimate: A value for an unknown parameter based on data.
Sampling Distribution: The probability distribution of an estimate if the experiment were repeated. The spread of this distribution indicates the estimate’s variability.
We will explore the rationale for this choice and how to assess its reliability.
Fitting a Poisson Distribution
In order to fit the Poisson distribution, we must estimate a value for λ from the observed data.
Since the average count in a 10-second interval was 8.352, we take this as an estimate of λ (recall that the E(X) = λ) and denote it by \(\hat{\lambda}\).
Let’s plot the observed and expected frequencies.
Code
#| fig-width: 8#| fig-height: 5#| code-fold: trueexpected_counts =dpois(0:17, lmb*10)*1207# expected countexpected_counts =c(sum(expected_counts[1:3]), expected_counts[4:18])# Create a data frame for plottingresults <-data.frame(Particles = d$n,Observed = d$Observed,Expected = expected_counts)# Plottinglibrary(ggplot2)ggplot() +geom_bar(data = results, aes(x = Particles, y = Observed), stat ="identity", fill ="skyblue", alpha =0.7, width=0.8) +geom_line(data = results, aes(x = Particles, y = Expected), color ="red", linewidth=1) +geom_point(data = results, aes(x = Particles, y = Expected), color ="red", size=3) +scale_x_continuous(breaks =0:17) +labs(title ="Alpha Particle Emissions: Observed vs. Expected Frequencies",subtitle =paste("Poisson Fit with lambda =", round(lmb, 3)),x ="Number of Alpha Particles per Interval",y ="Frequency") +theme_minimal()
Sampling Distribution of Estimates
If the experiment were to be repeated, the counts would be different and the estimate of λ would be different
It is thus appropriate to regard the estimate of λ as a random variable which has a probability distribution referred to as its sampling distribution
The situation is entirely analogous to tossing a coin 10 times and regarding the number of heads as a binomially distributed random variable
Understanding Variability of Estimates
Observing 6 heads in 10 coin tosses generates one realization of a random variable
In the same sense, 8.392 is a realization of a random variable
Key Question: What is the sampling distribution?
This is of practical interest, since the spread of the sampling distribution reflects the variability of the estimate
We could ask: to what decimal place is the estimate 8.392 accurate?
Later we will discuss the rationale for choosing to estimate λ as we have done
Normal Distribution
Code
# Data from https://www.openintro.org/data/index.php?data=satgpa&utm_source=chatgpt.comsat =read.csv("data/satgpa.csv")hist(sat$sat_sum, breaks =40, freq=F)
Normal Distribution
The density of normal is \[
f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}
\]
The mean and variance of the normal distribution are: \[
E(X) = \mu
\]
Now Plot the histogram and the fitted normal curve.
Code
# Draw bell curve on top of histogramx =seq(min(sat$sat_sum), max(sat$sat_sum), length.out =100)y =dnorm(x, mean = mu_hat, sd = sigma_hat)hist(sat$sat_sum, breaks =40, freq=F)lines(x, y, col ="red", lwd=3)
Parameter Estimation: The Core Idea
We treat observed data (\(x_1, x_2, \dots, x_n\)) as realizations of a random variable (\(X_1, X_2, \dots, X_n\)).
The joint distribution of these random variables depends on an unknown parameter, \(\theta\).
Often, we assume the \(X_i\) are independent and identically distributed (i.i.d.) from a distribution \(f(x|\theta)\).
An estimate of \(\theta\) is a function of the data, and thus is itself a random variable with a sampling distribution.
We will look at two general methods for creating estimates: 1. The Method of Moments 2. The Method of Maximum Likelihood
The Method of Moments
This method equates sample moments with the theoretical moments of the distribution, which are functions of the parameters.
Procedure:
Calculate the first few theoretical moments of the distribution, \(E(X^k)\), in terms of the parameters.
Express the parameters as functions of these moments.
Replace the theoretical moments with the sample moments, \(\hat{\mu}_k = \frac{1}{n}\sum_{i=1}^n X_i^k\), to get the estimates.
Method of Moments: Poisson Example
For a Poisson distribution, the first theoretical moment is the mean: \[ E(X) = \lambda \]
The first sample moment is the sample mean: \[ \hat{\mu}_1 = \bar{X} \]
Equating them gives the method of moments estimate for λ: \[ \hat{\lambda} = \bar{X} \]
This is precisely what was done in the alpha particle example.
Method of Moments: Normal Distribution
For a Normal distribution, the first two theoretical moments are: \[E(X) = \mu\]\[E(X^2) = \mu^2 + \sigma^2\] Solving for the parameters gives: \[\mu = E(X)\]\[\sigma^2 = E(X^2) - [E(X)]^2\] The method of moments estimates are found by replacing theoretical moments with sample moments: \[\hat{\mu} = \bar{X}\]\[\hat{\sigma}^2 = \frac{1}{n}\sum_{i=1}^n X_i^2 - \bar{X}^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2\]
Method of Moments: Body Temperature Example
Let’s fit a normal distribution to human body temperature data. We’ll use the method of moments to estimate μ and σ and visualize the fit.
Code
# Load the body temperature databodytemp <-read.csv("data/bodytemp.txt")temperatures <- bodytemp$temperature# Method of Moments estimatesmu_hat <-mean(temperatures)# Note: R's var() function uses the n-1 denominator for an unbiased estimate.# The strict method of moments uses n. The difference is negligible for large n.sigma2_hat <-var(temperatures) * (length(temperatures)-1) /length(temperatures)sigma_hat <-sqrt(sigma2_hat)print(paste("Estimated mu:", round(mu_hat, 2)))
# Plotting the histogram and the fitted normal curvelibrary(ggplot2)ggplot(data.frame(x = temperatures), aes(x = x)) +geom_histogram(aes(y =after_stat(density)), binwidth =0.2, fill ="skyblue", color ="white", alpha =0.7) +stat_function(fun = dnorm, args =list(mean = mu_hat, sd = sigma_hat), color ="red", size =1) +labs(title ="Fit of Normal Distribution to Body Temperatures",subtitle =paste("mu_hat =", round(mu_hat, 2), ", sigma_hat =", round(sigma_hat, 2)),x ="Body Temperature (F)",y ="Density") +theme_minimal()
Error = bias + variance
An estimate is wrong for two reasons.
Bias: even with infinite data you are pointed at the wrong place. \(E(\hat\theta)-\theta\) does not go to zero.
Variance: a new sample would move the answer. That is the sampling distribution we just discussed.
The \(p\)-th percentile is \(x_p\) such that \(F(x_p) = p/100\)
Median: \(F^{-1}(0.5)\)
3. Random Number Generation:
Inverse transform method: If \(U \sim \text{Uniform}(0,1)\), then \(X = F^{-1}(U)\) has CDF \(F\)
Inverse CDF Method
Algorithm for generating random samples:
Step 1: Draw \(U \sim \text{Uniform}[0,1]\)
Step 2: Set \(X = F^{-1}(U)\)
Where \(F^{-1}(u) = \inf\{x : F(x) = u\}\) is the inverse CDF (quantile function).
Why it works: If \(U\) is uniform on \([0,1]\), then \(F^{-1}(U)\) has the desired distribution \(F\).
This method is particularly useful when:
The inverse CDF can be computed analytically
We need exact samples from a known distribution
Normal Distribution CDF
For a normal random variable \(X \sim N(\mu, \sigma^2)\), the CDF is: \[F(x) = P(X \leq x) = \int_{-\infty}^x \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{(t-\mu)^2}{2\sigma^2}} dt\]
No closed-form solution exists - must use numerical methods or tables.
Standard Normal CDF: When \(\mu = 0, \sigma = 1\): \[\Phi(z) = \int_{-\infty}^z \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} dt\]
Hypothesis testing (p-values): Calculate probabilities of observing data as extreme or more extreme than what was observed, assuming a null hypothesis is true
Confidence intervals: Determine ranges of plausible values for population parameters based on sample data and desired confidence levels
Process control: Monitor manufacturing or business processes by identifying when measurements fall outside expected ranges
Risk analysis: Assess probabilities of extreme events in finance, insurance, and engineering applications
Quality control: Evaluate whether products meet specifications by calculating defect rates and acceptance probabilities
Medical diagnostics: Determine reference ranges for biomarkers and assess probability of disease given test results
A/B testing: Compare conversion rates between different website versions or marketing strategies
Portfolio optimization: Calculate value-at-risk and expected returns in financial modeling
Quantile Function: Definition and Intuition
The quantile function (also called the inverse CDF) gives the value \(x\) such that the probability of the random variable being less than or equal to \(x\) is a specified probability \(p\).
Definition: For a random variable \(X\) with CDF \(F(x)\), the quantile function \(Q(p)\) is: \[
Q(p) = \inf \{ x : F(x) \geq p \}
\]
Interpretation:\(Q(p)\) is the value below which a proportion \(p\) of the data falls.
Special cases:
\(Q(0.5)\) is the median
\(Q(0.25)\) and \(Q(0.75)\) are the first and third quartiles
Why is it useful? - Summarizes distributions (e.g., percentiles, medians) - Used in constructing confidence intervals and tolerance intervals - Important in risk analysis (e.g., Value-at-Risk in finance)
Quantile Function in R: Examples
R function:qnorm() for the normal distribution (and q* for other distributions)
Code
qnorm(0.975) # 97.5th percentile of standard normal
[1] 1.959964
Code
qnorm(0.95, mean =1000, sd =200) # 95th percentile of N(1000, 200²)
[1] 1328.971
Code
curve(qnorm(x), from =0, to =1, n =200,xlab ="Probability (p)", ylab ="Quantile Q(p)",main ="Quantile Function of Standard Normal (qnorm)")abline(h =0, v =0.5, col ="red", lty =2)
Summary: PDF, CDF, Quantile Function
PDF
Definition: Describes the relative likelihood for a random variable to take on a given value.
Notation:\(f(x)\)
Properties:
\(f(x) \geq 0\)
\(\int_{-\infty}^{\infty} f(x) dx = 1\)
Interpretation: Height of the curve at \(x\) (not a probability itself for continuous RVs).
CDF
Definition: Probability that the random variable is less than or equal to \(x\).
Notation:\(F(x) = P(X \leq x)\)
Properties:
\(F(x)\) is non-decreasing
\(F(-\infty) = 0\), \(F(\infty) = 1\)
Relation to PDF:\(F(x) = \int_{-\infty}^x f(t) dt\)
Quantile
Definition: Value \(x\) such that \(P(X \leq x) = p\)
Assume we observe \(y_1, \ldots, y_n\) which are independent draws from a random variable \(Y\) with probability density function \(p(y | \theta)\) (normal, binomial, etc.).
Likelihood is a function of \(\theta\) that measures how likely the observed data is given the parameter \(\theta\).
Maximum Likelihood Estimate
The maximum likelihood estimate of \(\theta\) is the value of \(\theta^*\) that maximizes the likelihood function.
\[\theta^* = \arg\max_{\theta} L(\theta)\] Sometimes we use notations \(\theta_{\text{MLE}}\) or simply \(\hat \theta\) to denote the maximum likelihood estimate of \(\theta\).
Or the log-likelihood function. \[\theta^* = \arg\max_{\theta} l(\theta)\]
Logarithm is a monotonically increasing function, so the maximum likelihood estimate of \(\theta\) is the same as the maximum log-likelihood estimate of \(\theta\).
Take the derivative of \(l(\mu)\) with respect to \(\mu\) and set it to 0: \[\frac{\partial l(\mu)}{\partial \mu} = 2\sum_{i=1}^n (y_i - \mu) = 0\]\[\sum_{i=1}^n (y_i - \mu) = 0\]\[\sum_{i=1}^n y_i = n\mu\]\[\hat \mu = \frac{1}{n}\sum_{i=1}^n y_i\]
Example: MLE for Normal Mean with Real Data
Dataset: SAT score for 1000 students
Code
# Data from https://www.openintro.org/data/index.php?data=satgpa&utm_source=chatgpt.comsat =read.csv("data/satgpa.csv")hist(sat$sat_sum, breaks =40, freq=F)
Maximum Likelihood Estimates for Normal Distribution
Parameter
Value
μ (mean)
103.33
σ² (variance)
203.92
σ (std dev)
14.28
n (sample size)
1000.00
Draw bell curve on top of histogram
Code
# Draw bell curve on top of histogramx =seq(min(sat$sat_sum), max(sat$sat_sum), length.out =100)y =dnorm(x, mean = mu_mle, sd = sigma_mle)hist(sat$sat_sum, breaks =40, freq=F)lines(x, y, col ="red", lwd=3)
Use Optimization to find MLE
Code
# Define negative log-likelihood functionneg_log_likelihood <-function(params, data) { mu <- params[1] sigma <- params[2]# Ensure sigma is positiveif (sigma <=0) return(Inf)# Calculate negative log-likelihood n <-length(data)-sum(dnorm(data, mean = mu, sd = sigma, log =TRUE))}# Optimize using built-in optimizerresult <-optim(c(50, 10), neg_log_likelihood, data = sat$sat_sum)mu_optim <- result$par[1]sigma_optim <- result$par[2]# Create results table for optimizationoptim_table <-data.frame(Parameter =c("μ", "σ"),Analytical =c(round(mu_mle, 2), round(sigma_mle, 2)),Optimized =c(round(mu_optim, 2), round(sigma_optim, 2)),Difference =c(round(abs(mu_mle - mu_optim), 6), round(abs(sigma_mle - sigma_optim), 6)))kable(optim_table, caption ="Comparison of Analytical vs Numerical Optimization Results")
Comparison of Analytical vs Numerical Optimization Results
Parameter
Analytical
Optimized
Difference
μ
103.33
103.33
0.000542
σ
14.28
14.28
0.002433
Likelihood Function Visualization
Create log-likelihood surface
# Create grid of mu values around the MLElibrary(ggplot2)mu_range <-seq(mu_mle -10, mu_mle +10, length.out =100)log_likelihood_values <-sapply(mu_range, function(mu) {-neg_log_likelihood(c(mu, sigma_mle), sat$sat_sum)})# Plot log-likelihoodlikelihood_df <-data.frame(mu = mu_range, log_likelihood = log_likelihood_values)ggplot(likelihood_df, aes(x = mu, y = log_likelihood)) +geom_line(color ="darkblue", linewidth =1) +geom_vline(xintercept = mu_mle, color ="red", linetype ="dashed", linewidth =1) +annotate("point", x = mu_mle, y =max(log_likelihood_values), color ="red", size =3) +labs(title ="Log-Likelihood Function for μ",subtitle =paste("Maximum at μ =", round(mu_mle, 2)),x ="Mean Score (μ)", y ="Log-Likelihood") +theme_minimal()
# Show how CI width changes with sample sizesample_sizes =c(100, 400, 1600, 6400)ci_widths =2* z_score * true_sd /sqrt(sample_sizes)cat("\nCI widths for different sample sizes:\n")
n = 100 → CI width = 0.0392 grams
n = 400 → CI width = 0.0196 grams
n = 1600 → CI width = 0.0098 grams
n = 6400 → CI width = 0.0049 grams
The Cost of Precision: - To get 4x more precise (0.08 → 0.02 grams), need 16× more data! - This is why large-scale experiments are expensive but necessary
Central Limit Theorem
The Most Powerful Theorem in Statistics
CLT states: Given sufficiently large sample size, sample means are approximately normally distributed, regardless of population distribution shape.
De Moivre-Laplace Theorem (1738): For independent Rademacher variables \(x_i\) with \(P(x_i = \pm 1) = 1/2\):
cat("Margin of error:", round(z_critical * se_mu, 2))
Margin of error: 0.89
Bernoulli Likelihood
For an i.i.d. sample \(X_1, \dots, X_n\) from a Bernoulli(p) distribution, the likelihood is: \[ l(p) = \prod_{i=1}^n p^{x_i} (1-p)^{1-x_i} \]
The log-likelihood is: \[ l(p) = \sum_{i=1}^n \ln\left(p^{x_i} (1-p)^{1-x_i}\right) \]\[ = \sum_{i=1}^n \left( x_i \ln(p) + (1-x_i) \ln(1-p) \right) \]
Bernoulli MLE
To find the maximum, we take the derivative with respect to p and set it to zero: \[ \frac{\partial l}{\partial p} = \frac{1}{p}\sum_{i=1}^n x_i - \frac{1}{1-p}\sum_{i=1}^n (1-x_i) = 0 \]\[ \implies \hat{p}_{MLE} = \frac{1}{n}\sum_{i=1}^n x_i = \bar{X} \]
MLE: Poisson Example
For an i.i.d. sample \(X_1, \dots, X_n\) from a Poisson(λ) distribution, the log-likelihood is: \[ l(\lambda) = \sum_{i=1}^n \ln\left(\frac{\lambda^{X_i} e^{-\lambda}}{X_i!}\right) = (\ln \lambda) \sum X_i - n\lambda - \sum \ln(X_i!) \]
To find the maximum, we take the derivative with respect to λ and set it to zero: \[ \frac{dl}{d\lambda} = \frac{1}{\lambda}\sum X_i - n = 0 \]\[ \implies \hat{\lambda}_{MLE} = \frac{\sum X_i}{n} = \bar{X} \]
For the Poisson distribution, the Method of Moments and the Maximum Likelihood estimates are the same.
Gamma Distribution
The density of the Gamma distribution is: \[
f(x) = \frac{1}{\Gamma(\alpha)\lambda^\alpha} x^{\alpha-1} e^{-\frac{x}{\lambda}}
\]
The mean and variance of the Gamma distribution are: \[
E(X) = \alpha\lambda
\]\[
Var(X) = \alpha\lambda^2
\]
Examples of Gamma Distribution
The Gamma distribution is extremely versatile and appears in many real-world applications:
Examples:
Waiting times: Time until the k-th event in a Poisson process (e.g., time until the 5th customer arrives)
Reliability engineering: Time to failure of components or systems
Insurance: Claim amounts and time between claims
Weather: Rainfall amounts, wind speeds
Finance: Stock price movements, portfolio losses
Medical: Time to recovery, survival times
Queue theory: Service times in complex systems
Let’s examine a real dataset: inter-arrival times of cosmic ray detections.
MLE: Gamma Distribution Example
For the Gamma distribution, the log-likelihood is: \[l(\alpha, \lambda) = n\alpha\ln(\lambda) - n\ln(\Gamma(\alpha)) + (\alpha-1)\sum\ln(X_i) - \lambda\sum X_i\] The MLEs for α and λ don’t have a simple closed-form solution and must be found numerically by maximizing this function.
Let’s fit a Gamma distribution to the inter-arrival times of gamma-rays.
MLE: Gamma Distribution Example
Code
#| fig-width: 8#| fig-height: 5# Load the datagamma_arrivals <-read.csv("data/gamma-arrivals.txt", header =FALSE)$V1# Log-likelihood function for the Gamma distributionlog_likelihood_gamma <-function(params, data) { alpha <- params[1] lambda <- params[2]if (alpha <=0|| lambda <=0) return(NA) # Parameters must be positive n <-length(data) loglik <- n * alpha *log(lambda) - n *lgamma(alpha) + (alpha -1) *sum(log(data)) - lambda *sum(data)return(-loglik) # optim minimizes by default, so we return negative log-likelihood}# Use Method of Moments for starting valuesmean_g <-mean(gamma_arrivals)var_g <-var(gamma_arrivals)alpha_start <- mean_g^2/ var_glambda_start <- mean_g / var_g# Find MLEs using numerical optimizationmle_fit <-optim(par =c(alpha_start, lambda_start), fn = log_likelihood_gamma, data = gamma_arrivals)alpha_mle <- mle_fit$par[1]lambda_mle <- mle_fit$par[2]print(paste("MLE for alpha:", round(alpha_mle, 4)))
[1] "MLE for alpha: 1.0263"
Code
print(paste("MLE for lambda:", round(lambda_mle, 4)))
[1] "MLE for lambda: 0.0128"
Code
# Plottingggplot(data.frame(x = gamma_arrivals), aes(x = x)) +geom_histogram(aes(y =after_stat(density)), bins =50, fill ="skyblue", color ="white", alpha =0.7) +stat_function(fun = dgamma, args =list(shape = alpha_mle, rate = lambda_mle), color ="red", linewidth =1) +labs(title ="Fit of Gamma Distribution to Gamma-Ray Inter-Arrival Times",subtitle =paste("alpha_mle =", round(alpha_mle, 3), ", lambda_mle =", round(lambda_mle, 3)),x ="Inter-Arrival Time (seconds)",y ="Density") +xlim(0, 1000) +theme_minimal()
Properties of Maximum Likelihood Estimators
Under general “smoothness” conditions, MLEs have desirable large-sample properties:
Consistency: As the sample size \(n \to \infty\), the MLE \(\hat{\theta}_n\) converges in probability to the true parameter value \(\theta_0\).
Asymptotic Normality: The distribution of \(\sqrt{n}(\hat{\theta}_n - \theta_0)\) approaches a Normal distribution with mean 0.
Asymptotic Efficiency: For large samples, the MLE achieves the smallest possible variance among a large class of unbiased estimators (the Cramér-Rao lower bound).
This provides a strong theoretical justification for using MLEs.
Confidence Intervals from MLEs
The asymptotic normality of MLEs allows us to construct approximate confidence intervals. The variance of the limiting normal distribution is given by \(1/I(\theta_0)\), where \(I(\theta)\) is the Fisher Information.
An approximate \(100(1-\alpha)\%\) confidence interval for \(\theta\) is: \[ \hat{\theta} \pm z_{\alpha/2} \frac{1}{\sqrt{nI(\hat{\theta})}} \] The term \(\frac{1}{\sqrt{nI(\hat{\theta})}}\) is the estimated standard error of the MLE.
MLE: Negative Binomial Example
The Negative Binomial distribution is often used for count data that is “overdispersed” (variance is greater than the mean), making it a good alternative to the Poisson. Let’s fit it to data on potato beetle counts.
Code
#| message: false#| warning: false#| fig-width: 8#| fig-height: 5library(MASS) # For fitdistr# Load data and expand from frequency tableevans <-read.csv("data/evans1953.txt")beetle_data <-rep(evans$Count, evans$PotatoBeetles)# Use fitdistr to find MLEs for the negative binomial distribution# 'size' is r (number of successes), 'mu' is the meanneg_bin_fit <-fitdistr(beetle_data, "negative binomial")r_hat <- neg_bin_fit$estimate["size"]mu_hat <- neg_bin_fit$estimate["mu"]print(paste("MLE for r (size):", round(r_hat, 2)))
[1] "MLE for r (size): 2.24"
Code
print(paste("MLE for mean (mu):", round(mu_hat, 2)))
[1] "MLE for mean (mu): 4.7"
Code
# Calculate expected frequenciesmax_count <-max(evans$Count)counts_range <-0:max_countexpected_freq <-length(beetle_data) *dnbinom(counts_range, size = r_hat, mu = mu_hat)# Create a data frame for plottingplot_data <-data.frame(Count = counts_range,Observed =sapply(counts_range, function(c) sum(beetle_data == c)),Expected = expected_freq)# Plotting observed vs. expectedggplot(plot_data, aes(x = Count)) +geom_bar(aes(y = Observed), stat ="identity", fill ="skyblue", alpha =0.7, width =0.8) +geom_line(aes(y = Expected), color ="red", size =1) +geom_point(aes(y = Expected), color ="red", size =3) +scale_x_continuous(breaks =seq(0, max_count, by =2)) +labs(title ="Potato Beetle Counts: Observed vs. Expected Frequencies",subtitle =paste("Negative Binomial Fit with r_hat =", round(r_hat, 2), ", mu_hat =", round(mu_hat, 2)),x ="Number of Beetles per Plant",y ="Frequency") +theme_minimal()
MLE For LLM
Large Language Models (LLMs) like GPT, BERT, and others are trained using Maximum Likelihood Estimation as a core principle.
The Setup: - Given a sequence of tokens \(w_1, w_2, \ldots, w_n\) (words, subwords, or characters) - Model parameterized by \(\theta\) (neural network weights) - Goal: Find \(\theta\) that maximizes the likelihood of observing the training text
The Likelihood: For autoregressive models (like GPT), we factorize the joint probability: \[P(w_1, w_2, \ldots, w_n | \theta) = \prod_{i=1}^n P(w_i | w_1, \ldots, w_{i-1}, \theta)\]
The MLE objective is to maximize this likelihood over the training corpus.
