Statistics for Applications - MIT 18.650

MIT 18.650 Statistics for Applications, Fall 2016 - YouTube

MIT 18.650, taught by Professor Philippe Rigollet, provides a rigorous, mathematical introduction to the fundamental principles of statistics. The course is designed not as a guide to applied data analysis, but as a theoretical foundation that explains why statistical methods work, what their guarantees are, and under what assumptions they are valid. It rigorously distinguishes statistics as the "inverse problem"—inferring the properties of an unknown data-generating process from observed data—from probability, the "forward problem" of predicting data from a known process. The central theme is the development of a principled framework for estimation, decision-making, and modeling in the face of uncertainty.
Part I: The Foundations of Parametric Inference
The course begins by establishing the core concepts of statistical modeling, assuming data is drawn IID (Independent and Identically Distributed) from a parametric family of distributions. The initial focus is on estimation and confidence intervals. Key foundational tools are introduced:

• The Central Limit Theorem (CLT): The workhorse of asymptotic statistics, used to construct confidence intervals for means and proportions by providing a "pivotal quantity" whose distribution is known and independent of unknown parameters.
• Modes of Convergence: A formal review of convergence in probability, almost sure convergence, and convergence in distribution provides the theoretical language to justify statistical procedures.
• Slutsky's Theorem and the Delta Method: These crucial theorems allow for the practical application of asymptotic results. Slutsky's Theorem justifies the "plug-in" principle (replacing unknown parameters in a variance term with their estimates), while the Delta Method provides a CLT for functions of estimators, dramatically expanding the scope of inference.

Part II: Principles of Estimation
The course then delves into the theory of finding "good" estimators. Moving beyond simple sample averages, it introduces two major principled frameworks for deriving estimators:

• Maximum Likelihood Estimation (MLE): Presented as the premier method of estimation, the MLE is rigorously derived from the fundamental goal of finding a probability distribution that is "close" to the true data-generating distribution. Using Kullback-Leibler (KL) divergence as a surrogate for statistical distance, the principle of maximizing the likelihood function is shown to be equivalent to minimizing the estimated KL divergence.
• The Method of Moments (MoM): An alternative and often simpler method based on the principle of equating population moments to their empirical (sample) counterparts and solving for the model parameters.
• Evaluating Estimators: The performance of estimators is formally evaluated using the concepts of bias, variance, and quadratic risk (Mean Squared Error - MSE). The fundamental bias-variance decomposition (MSE = Bias² + Variance) is introduced as the central trade-off in estimator design. The Fisher Information is defined as a key measure of the "informativeness" of a statistical model, directly linking the curvature of the log-likelihood function to the asymptotic variance and precision of the MLE.

Part III: Hypothesis Testing and Model Validation
The focus shifts from estimation to formal decision-making through hypothesis testing.

• The Neyman-Pearson Framework: The course establishes the classic framework of null (H₀) and alternative (H₁) hypotheses, emphasizing their asymmetric roles ("innocent until proven guilty"). The concepts of Type I and Type II errors, the level (α), and the power of a test are defined. The paradigm of fixing the level and maximizing the power is established as the standard for designing optimal tests. The p-value is introduced as the primary tool for reporting the strength of evidence against the null hypothesis.
• A Catalog of Tests: Several key parametric tests are developed:
  • Student's t-test: For inference on means when the variance is unknown, relying on the t-distribution which accounts for the extra uncertainty in estimating the variance.
  • Wald Test and Likelihood Ratio Test (LRT): Two powerful and general tests based on the MLE for testing hypotheses in multi-parameter settings. Their test statistics are shown to have asymptotic Chi-Squared (χ²) distributions.
• Goodness-of-Fit: To validate the underlying distributional assumptions of a model, goodness-of-fit tests are explored. The Chi-Squared test is used for discrete data, while the Kolmogorov-Smirnov (KS) test is developed for continuous data, based on measuring the maximum distance between the Empirical CDF and the hypothesized CDF. The Lilliefors test is introduced as a necessary modification for testing normality when parameters are estimated. Quantile-Quantile (QQ) plots are presented as an essential visual diagnostic tool.

Part IV: Modeling Relationships and Discovering Structure
The course applies these foundational principles to two major areas of statistics.

• Linear Regression: This section covers the fundamental model for relating a response variable y to one or more predictor variables x. The derivation proceeds from the simple univariate case to the full multivariate model using matrix algebra.
  • Estimation: The method of least squares is derived and shown to be equivalent to the MLE under Gaussian noise. The geometric interpretation of least squares as an orthogonal projection is heavily emphasized.
  • Inference: Statistical inference for regression coefficients is developed, including t-tests for the significance of individual predictors and the F-test for joint hypotheses about multiple coefficients. The problem of multiple testing and the Bonferroni correction are addressed.
• Principal Component Analysis (PCA): As the primary example of unsupervised learning, PCA is introduced as a technique for dimension reduction and data visualization. The goal is not to predict a response, but to find the inherent structure in a high-dimensional dataset. PCA is derived as a procedure to find the directions of maximal variance in the data through the eigen-decomposition of the sample covariance matrix. The lecture also touches upon the challenges of interpretation and the modern extension to Sparse PCA.

Part V: Unifying Framework and Alternative Paradigms
The final part of the course introduces an alternative statistical philosophy and a powerful framework that unifies many of the preceding topics.

• Bayesian Statistics: This section presents the Bayesian paradigm as a fundamentally different approach to inference. The core idea is to treat the unknown parameter as a random variable and update a prior belief into a posterior distribution using Bayes' Rule. The concepts of conjugate priors (e.g., Beta-Bernoulli), uninformative priors (including Jeffreys prior), and Bayesian confidence regions are explored, highlighting the different interpretation compared to their frequentist counterparts. The Bernstein-von Mises theorem is mentioned to show the asymptotic convergence of Bayesian and frequentist results.
• Generalized Linear Models (GLMs): This capstone topic provides a powerful unification of regression modeling. GLMs extend linear regression by:
  1. Allowing the response variable's distribution to be any member of the exponential family (e.g., Bernoulli, Poisson), connecting directly back to the theory of MLE and Fisher Information.
  2. Introducing a link function that relates the mean of the response to the linear predictor, enabling the modeling of non-linear relationships.
     The course concludes by detailing the estimation procedure for GLMs, showing how the Newton-Raphson and Fisher Scoring algorithms lead to the elegant and standard Iteratively Reweighted Least Squares (IRLS) method, which solves a complex non-linear problem through a sequence of familiar weighted least squares regressions.

L1: Introduction

The lecture begins with an introduction to the course, 18.650 Fundamentals of Statistics, by the instructor, Professor Philippe Rigollet. He clarifies that the course is a mathematical introduction to the field, focusing on providing theoretical guarantees for statistical methods rather than on applied data analysis or "statistical thinking."
The primary goals of the course are for students to:

1. Formulate real-world questions in precise mathematical terms. For example, translating "is a drug effective?" into a question about a specific parameter of a statistical model.
2. Learn a catalog of fundamental statistical methods to apply to these formulated problems.
3. Understand the assumptions and limitations of each method to make sound and critical conclusions from data.
4. Gain the foundational knowledge necessary for more advanced, graduate-level courses in statistical machine learning.

Course Logistics and Requirements
The instructor outlines the practical details and expectations for the course.

• Lectures and Recitations: Lectures are Tuesdays and Thursdays. Recitations, held on Wednesdays, are mandatory and will focus on problem-solving techniques.
• Homework: There will be 11 weekly problem sets, with the 10 best counting for 30% of the final grade. Submissions are due on Mondays via Stellar as readable PDF files. Students are allowed two 24-hour late submissions without penalty.
• Exams:
  • Midterms: There will be two in-class midterms. To reduce pressure, only the best of the two will be kept, counting for 30% of the grade. They are closed-book, but a two-sided, letter-sized cheat sheet is permitted.
  • Final Exam: A three-hour exam counting for 40% of the grade. It will be open-notes but not open-book.
• Prerequisites: A strong background in probability (equivalent to 18.600 or 6.041B) is a firm prerequisite. Calculus and basic linear algebra are also required.
• Textbook and Materials: There is no required textbook. The recommended text is "All of Statistics" by Wasserman. All necessary material will be provided through lecture slides, annotations, and video recordings available on Stellar.

Motivation: Why Study Statistics?
Professor Rigollet motivates the study of statistics by highlighting its ubiquity in science, media, and industry.

• Interpreting Studies: News headlines often report findings from scientific studies (e.g., "study finds traffic is bad for your health"). Statistics provides the tools to critically evaluate these claims, understand what they truly mean (e.g., they are statements about populations, not individuals), and recognize their underlying assumptions.
• Avoiding Statistical Flukes: The famous "dead salmon" fMRI experiment is used as a cautionary tale. Without proper statistical methods to correct for searching through massive amounts of data (multiple comparisons), it is easy to find meaningless patterns in random noise. The phenomenon of "p-hacking" is also mentioned, where researchers might misuse data to find publishable (but false) results.
• Real-World Applications: Statistics is foundational to many fields:
  • Hydrology and Engineering: The historical example of Dutch dikes shows how statistics is used to make decisions under uncertainty based on past data.
  • Insurance: Actuarial science uses statistics to calculate premiums based on the expected risk associated with events like car accidents or broken phones.
  • Clinical Trials: The FDA relies on rigorous statistical protocols to determine the efficacy and safety of new drugs.
  • Genetics and Medicine: Statistics is essential for identifying associations between genes and diseases like Alzheimer's from large datasets.

The Core Distinction: Probability vs. Statistics
The central theme of the lecture is the fundamental difference between probability and statistics.

• Probability (The "Forward Problem"): In probability, you are given the truth—a complete description of a random process and its parameters (e.g., "a fair six-sided die is rolled"). Your task is to predict the properties of the data that will be generated (e.g., "the probability of rolling a 7 with two dice is 1/6").
• Statistics (The "Inverse Problem"): In statistics, you work backward. You are given the data—a set of observations (e.g., "a die was rolled 100 times, resulting in these outcomes"). Your task is to infer the truth—the underlying process or its unknown parameters (e.g., "is the die fair? What is the probability of rolling a 6?").
• The Role of Modeling: Since the "truth" is unknown, statistics relies on making simplifying modeling assumptions. A complex process is often modeled as a simple process plus random noise. The validity of any statistical conclusion depends heavily on the quality of these assumptions. The instructor emphasizes the "golden rule": Statistics is about replacing expectations with averages.

A Case Study: The Kissing Couples Experiment
A study from the journal Nature on whether kissing couples tend to turn their heads to the right is used as a running example.

• The Question: Is there a natural preference for turning right? This is translated into a mathematical question about an unknown parameter p, the true proportion of all couples that turn right.
• The Data: The study observed n=124 couples, and 80 of them turned right.
• The Estimate: The statistical estimate for p is p̂ = 80/124 ≈ 64.5%. p̂ (p-hat) is the estimate, while p is the true, unknown parameter.
• The Statistical Challenge: Is this 64.5% result "significant," or could it have arisen by chance if the true proportion p was actually 50%? The sample size (n=3 vs. n=124) is crucial in determining how much confidence we have in the estimate.

The Statistical Modeling Process
The kissing problem is used to demonstrate the formal steps of statistical modeling.

1. Define Random Variables: The outcome for each couple is modeled as a random variable. R_i is an indicator variable: 1 if the i-th couple turns right, 0 otherwise.
2. Define the Estimator: The estimator for p is the sample average, p̂ = (1/n) * Σ R_i. An estimator is a function of the random data and is therefore itself a random variable. The specific number it produces from a given dataset is the estimate.
3. State the Modeling Assumptions (The IID Model):
   • Each R_i is a Bernoulli(p) random variable.
   • The parameter p is the same for all i (the data is Identically Distributed). This is a crucial assumption that allows us to pool information across observations.
   • The R_i variables are mutually Independent (one couple's choice doesn't affect another's).
   • Together, these form the IID (Independent and Identically Distributed) model, a cornerstone of classical statistics.

Practical Exercise: Estimating Mean and Variance
The lecture concludes with a quick, practical exercise to solidify the concept of estimation.

• Problem: Given 15 student grades assumed to be drawn from a Gaussian (Normal) distribution with unknown mean μ and variance σ², find estimates for μ and σ.
• Estimator for the Mean The sample average of the grades.
• Estimator for the Variance The sample variance, which is the average of the squared differences between each grade and the sample mean. This is another application of replacing the true expectation (μ) with its average (the sample mean μ̂).

L2: Introduction cont.

The lecture begins with a review of how confidence intervals were constructed in the previous session using the Central Limit Theorem (CLT), taking the Bernoulli "kissing couples" example as a reference.

• The Setup: For IID Bernoulli(p) observations, the sample mean X̄n is used to estimate the unknown probability p.
• The Pivotal Quantity: The CLT states that the standardized sample mean converges in distribution to a standard normal (Gaussian) distribution:
  √n * (X̄n - p) / √(p(1-p)) → N(0, 1)
• The Key Idea: Pivotal Distributions: The instructor emphasizes that the power of the CLT lies in producing a pivotal quantity. This is a function of the data and the unknown parameter whose limiting distribution (N(0,1)) is known and does not depend on any unknown parameters. This allows us to use standard tables (like the quantiles q_(α/2) of a Gaussian) to make probabilistic statements.
• Deriving the Interval: By setting the probability of this pivotal quantity being within [-q_(α/2), q_(α/2)] to 1-α, we can algebraically solve for p to obtain a confidence interval. This process, however, can lead to complex equations (e.g., a quadratic equation in p).

2. Alternatives to the CLT: Exact Calculations and Hoeffding's Inequality
The lecture explores alternatives to the asymptotic CLT approach.

• Exact Calculation (for small For a very small sample size (e.g., n=3), the distribution of the sample mean X̄n is simply a scaled Binomial distribution. One could, in principle, compute the exact probabilities for each possible outcome. However, these probabilities would still depend on the unknown parameter p, making it difficult to invert and form a confidence interval without further techniques.
• Hoeffding's Inequality (A Non-Asymptotic Bound): For bounded random variables (like Bernoulli, which are bounded in [0, 1]), Hoeffding's inequality provides a powerful, non-asymptotic bound on the deviation of the sample mean from its expectation:
  P(|X̄n - p| > ε) ≤ 2 * exp(-2nε²)
  • Key Advantage: This inequality holds for any sample size n, not just large ones.
  • Deriving a Confidence Interval: We can set the right-hand side equal to α and solve for ε. This gives an explicit width for a 1-α confidence interval: ε = √(log(2/α) / 2n). The interval is then p̂ ± ε.
  • Trade-off: The confidence interval derived from Hoeffding's inequality is generally wider (more conservative) than the one from the CLT. The CLT leverages more specific knowledge about the limiting distribution's shape (Gaussian), while Hoeffding's provides a more universal but looser bound.

3. The Theoretical Foundation: Modes of Convergence
To properly justify statistical procedures like the "plug-in" principle (replacing unknown parameters with their estimates), the lecture pivots to a more formal discussion on the different ways a sequence of random variables can converge.

• Four Modes of Convergence:
  1. Almost Sure (a.s.) Convergence: The strongest form. Xn → X means that for almost every possible underlying outcome, the sequence of real numbers Xn(ω) converges to X(ω). This is the basis for the Strong Law of Large Numbers.
  2. Convergence in Probability: A weaker form. Xn → X in probability means that the probability of Xn and X being more than a small distance ε apart goes to zero as n → ∞. This is the basis for the Weak Law of Large Numbers.
  3. Convergence in L_p: Xn → X in L_p means that the expected p-th power of the distance between Xn and X goes to zero. Convergence in L_2 is particularly important as it relates to variance.
  4. Convergence in Distribution: The weakest form, but the most important for statistics. Xn → X in distribution means that the cumulative distribution functions (CDFs) of Xn converge to the CDF of X. The Central Limit Theorem is a statement about convergence in distribution.
• The Hierarchy: These modes are related in a strict hierarchy: Almost Sure convergence and L_p convergence both imply convergence In Probability, which in turn implies convergence In Distribution.

4. Slutsky's Theorem: The Statistician's Best Friend
While convergence in distribution is very useful, it is too weak to allow for simple algebraic manipulation (e.g., the limit of a sum is not necessarily the sum of the limits). Slutsky's Theorem provides a critical exception that makes many statistical procedures valid.

• Slutsky's Theorem: If Xn converges in distribution to X, and Yn converges in probability to a constant , then:
  • Xn + Yn converges in distribution to X + c.
  • Xn * Yn converges in distribution to X * c.
• This theorem is essential because it allows us to combine a result from the CLT (which gives convergence in distribution) with an estimate that is consistent (converges in probability to a constant), such as a sample mean converging to the true mean.

5. New Example: Inter-arrival Times of the T
A new problem is introduced to apply these theoretical tools: modeling the inter-arrival times of the T at Kendall Square.

• The Model: The inter-arrival times T₁, T₂, ..., Tn are assumed to be IID Exponential(λ) random variables.
• The Goal: Estimate the rate parameter λ and construct a confidence interval for it.
• The Estimator: The expectation of an Exponential(λ) variable is 1/λ. Following the principle of "replacing expectations with averages," a natural estimator for 1/λ is the sample mean T̄n. Therefore, the estimator for λ is λ̂ = 1 / T̄n.

6. The Delta Method: CLT for Functions of Averages
The CLT applies directly to the sample mean T̄n, but we need a result for our estimator λ̂ = 1 / T̄n, which is a function of the sample mean. The Delta Method provides the tool to do this.

• The Idea: Using a first-order Taylor expansion, we can approximate a function g(T̄n) around the true mean θ = 1/λ.
  g(T̄n) ≈ g(θ) + (T̄n - θ) * g'(θ)
• The Result: The Delta Method states that if the CLT holds for √n(T̄n - θ), then a CLT also holds for √n(g(T̄n) - g(θ)). The limiting distribution is still Gaussian with mean 0, but the variance is scaled by the square of the derivative:
  √n * (g(T̄n) - g(θ)) → N(0, σ² * [g'(θ)]²)
• Application: For λ̂ = 1 / T̄n, the function is g(x) = 1/x, so g'(x) = -1/x². Applying the Delta Method yields a CLT for λ̂.

7. Justifying the "Plug-in" Principle with Slutsky's Theorem
The confidence interval derived from the Delta Method will still contain the unknown parameter λ in its variance term. The final step is to replace this unknown λ with its estimate λ̂.

• The "Plug-in" Step: The standardized quantity √n(λ̂ - λ) / λ converges in distribution to N(0, 1). We want to replace the λ in the denominator with λ̂.
• The Justification: This step is rigorously justified by Slutsky's Theorem. We can rewrite the new quantity as:
  [√n(λ̂ - λ) / λ] * [λ / λ̂]
  • The first term converges in distribution to a standard normal N(0, 1) (by the Delta Method).
  • The second term, λ / λ̂, converges in probability to the constant 1 (by the Law of Large Numbers and the Continuous Mapping Theorem).
• According to Slutsky's Theorem, the product of these two terms converges in distribution to N(0, 1) * 1 = N(0, 1). This proves that using the "plug-in" estimate in the denominator is asymptotically valid, allowing us to construct a computable confidence interval for λ.

L3: Parametric Inference

The lecture begins with a recap of the previous session's main topic: constructing confidence intervals using the Central Limit Theorem (CLT).

• The Setup: Given a set of IID (Independent and Identically Distributed) Bernoulli(p) random variables X₁, ..., Xn, the goal is to estimate the unknown parameter p.
• The Estimator: The natural estimator for p is the sample mean, p̂ = X̄n.
• The Pivotal Quantity: The power of the CLT is that it provides a "pivotal" quantity—a function of the data and the parameter whose distribution converges to something known and independent of any unknown parameters. For Bernoulli variables, this quantity is:
  √n * (X̄n - p) / √(p(1-p)) → N(0, 1)
• Constructing the Interval: By using the known quantiles of the standard normal distribution (N(0,1)), we can form an inequality and solve for the unknown p to get a 1-α confidence interval. This process relies on the asymptotic nature of the CLT.

2. The Concept of a "Pivotal Distribution"
The instructor emphasizes the importance of finding a pivotal quantity. The key is that the limiting distribution (N(0,1) in the case of the CLT) does not depend on any unknown parameters. This is what allows us to use standard tables or software to find critical values (like quantiles q_(α/2)) and construct intervals. If the limiting distribution depended on p, the process would be circular.
3. Exact vs. Asymptotic Approaches
The lecture contrasts the asymptotic CLT approach with what would be required for an exact, finite-sample analysis.

• Exact Calculation ( For a small sample like n=3, the distribution of n * X̄n is exactly a Binomial(3, p). While we can write down the exact probabilities for each outcome, these probabilities depend on the unknown . Creating a confidence interval from this requires more complex inversion techniques and does not yield a simple, universal formula.
• Hoeffding's Inequality: As an alternative that is non-asymptotic (valid for any n), Hoeffding's inequality is introduced. For bounded random variables, it provides a universal upper bound on the tail probabilities of the sample mean:
  P(|X̄n - E[X]| > ε) ≤ 2 * exp(-2nε² / (b-a)²)
  • This can be used to construct a valid confidence interval by setting the right-hand side to α and solving for ε.
  • The Trade-off: The resulting interval is generally wider (more conservative) than the one from the CLT because Hoeffding's is a more general bound and doesn't use the specific information that the limit is Gaussian. Statistical assumptions (like the CLT) allow for sharper, more precise statements.

4. Review of Convergence Modes for Random Variables
To justify more advanced statistical techniques, the lecture formally reviews different modes of convergence for sequences of random variables. The key takeaway is the hierarchy of these modes.

• Strongest: Almost Sure Convergence (Xn → X a.s.).
• Strong: Convergence in L_p.
• Weaker: Convergence in Probability.
• Weakest: Convergence in Distribution.
• The Hierarchy: Almost Sure and L_p convergence both imply convergence in Probability, which in turn implies convergence in Distribution.
• Importance for Statistics: Convergence in distribution is the weakest but most frequently used mode in statistics, as it underpins the CLT. It states that the probabilities associated with Xn become close to those of X, which is sufficient for constructing confidence intervals.

5. Slutsky's Theorem and the Delta Method
Two powerful theorems are introduced that allow for the manipulation of these convergence results.

• Slutsky's Theorem: This is a crucial tool that allows for combining different modes of convergence. If Xn → X in distribution and Yn → c in probability (where c is a constant), then Xn + Yn → X + c and Xn * Yn → X * c in distribution. This theorem will be essential for justifying the "plug-in" principle.
• The Delta Method: This theorem provides a CLT for functions of sample means. If we have a CLT for an estimator Zn (e.g., a sample mean), the Delta Method gives us a CLT for g(Zn) where g is a differentiable function. The variance of the new limiting Gaussian distribution is scaled by the square of the derivative, [g'(θ)]².

6. Application: Estimating the Rate of the T
A new example is introduced: estimating the arrival rate λ of the T, assuming the inter-arrival times T₁, ..., Tn are IID Exponential(λ).

• The Estimator: The mean of an Exponential(λ) is 1/λ. Following the principle of replacing expectations with averages, the estimator for 1/λ is the sample mean T̄n. Therefore, the estimator for λ is λ̂ = 1 / T̄n.
• Deriving a Confidence Interval for
  1. The CLT applies directly to T̄n, giving a limiting normal distribution for √n(T̄n - 1/λ).
  2. We need a limit theorem for λ̂ = 1 / T̄n. The Delta Method is used with the function g(x) = 1/x. This yields a CLT for √n(λ̂ - λ).
  3. The variance of the resulting limiting distribution depends on the unknown λ. To create a computable confidence interval, we must replace the λ in the standard deviation term with its estimate, λ̂.
  4. This "plug-in" step is formally justified by Slutsky's Theorem. The standardized quantity can be written as the product of a term that converges in distribution to N(0,1) and a term (λ/λ̂) that converges in probability to the constant 1. By Slutsky's theorem, the product also converges in distribution to N(0,1), validating the procedure.

7. Introduction to Statistical Modeling
The second half of the lecture formalizes the concept of a statistical model.

• Purpose of Modeling: The goal is to simplify reality by assuming that the observed data comes from a restricted family of probability distributions. This allows us to pool information and estimate a small number of parameters instead of an entire, complex distribution. For example, instead of estimating the probability of every possible number of siblings, we can assume the data follows a Poisson distribution and only estimate its single rate parameter, λ.
• Formal Definition: A statistical model is a pair (E, {P_θ}), where E is the sample space (the set of possible outcomes) and {P_θ} is a family of probability distributions indexed by a parameter θ from a parameter space Θ.
• Key Assumptions:
  • Well-specified model: We assume the true data-generating distribution is one of the P_θ in our family.
  • Identifiability: We assume the mapping from θ to P_θ is injective. This means that if θ₁ ≠ θ₂, then P_{θ₁} ≠ P_{θ₂}. This ensures that if we can learn the distribution P, we can uniquely identify the parameter θ.
• Examples of Models: The lecture lists several common parametric models: Bernoulli, Exponential, Poisson, and Normal (Gaussian).

8. Estimators and Their Properties
The lecture concludes by defining estimators and the properties used to evaluate their quality.

• Estimator: An estimator θ̂ is a function of the data X₁, ..., Xn. It is a random variable. A specific value computed from a dataset is an estimate.
• Consistency: An estimator is consistent if it converges to the true parameter θ as the sample size n → ∞.
• Bias: The bias of an estimator is E[θ̂] - θ. An estimator is unbiased if its bias is zero. The sample mean X̄n is an unbiased estimator of the population mean.
• Quadratic Risk (Mean Squared Error - MSE): A key measure of an estimator's quality is its quadratic risk, or MSE, defined as E[(θ̂ - θ)²].
• Bias-Variance Decomposition: A fundamental result is that the MSE can be decomposed into two components:
  MSE(θ̂) = (Bias(θ̂))² + Var(θ̂)
  This shows that a good estimator must have both low bias (it's centered on the right value on average) and low variance (it doesn't fluctuate wildly from sample to sample). An estimator with low MSE will be consistent.

L4: Parametric Inference cont. and MLE

The lecture begins by reviewing the key properties used to evaluate the performance of a statistical estimator θ̂.

• Bias: Bias(θ̂) = E[θ̂] - θ. An estimator is unbiased if its bias is zero. This means that, on average, over repeated experiments, the estimator would give the true parameter value.
• Variance: Var(θ̂) = E[(θ̂ - E[θ̂])²]. This measures the estimator's fluctuation or volatility from one sample to another.
• Quadratic Risk (Mean Squared Error - MSE): The primary measure of an estimator's overall quality is its MSE, defined as E[(θ̂ - θ)²].
• Bias-Variance Decomposition: A fundamental result shows that the MSE can be broken down into two components:
  MSE(θ̂) = (Bias(θ̂))² + Var(θ̂)
  This decomposition highlights the inherent trade-off in estimator design: a good estimator must balance having low bias (being accurate on average) and low variance (being stable and not fluctuating wildly).
• Example (Bernoulli Estimators):
  • The sample mean X̄n is unbiased, and its variance (and thus its risk) decreases like 1/n.
  • The constant estimator 0.5 has zero variance but has a bias of 0.5 - θ. Its risk is (0.5 - θ)².
  • The estimator X₁ (using only the first observation) is unbiased but has a high, constant variance that does not improve with more data.
  • A plot of the risks shows that for the true parameter θ near 0.5, the constant estimator can be better than X₁, but as n grows, the sample mean X̄n becomes superior across almost the entire parameter space.

2. The Goal of Estimation: Finding a "Close" Distribution
The lecture transitions from evaluating estimators of a single parameter θ to the more fundamental goal of finding an entire probability distribution P_θ̂ that is "close" to the true, unknown data-generating distribution P_{θ*}.

• Measuring Distance Between Distributions: Total Variation (TV) Distance
  • The Total Variation (TV) distance is introduced as a natural way to measure the closeness of two probability distributions, P and Q.
  • Definition 1 (Worst-case probability difference): TV(P, Q) = max_A |P(A) - Q(A)| over all possible events A. This has a clear statistical interpretation: if TV(P, Q) is small, then any probability computed under P will be very close to the same probability computed under Q.
  • Definition 2 (L₁ distance between densities/PMFs): An equivalent, more computable formula is TV(P, Q) = (1/2) * ∫ |p(x) - q(x)| dx (or a sum for discrete distributions), where p and q are the densities or probability mass functions.
  • It is proven that these two definitions are equivalent by showing that the maximum difference in probabilities occurs for the event A* = {x | p(x) > q(x)}.

3. A Principled Strategy for Estimation
A general strategy for finding a good estimator is proposed:

1. Define a notion of distance (like TV distance) between the true distribution P_{θ*} and a candidate distribution P_θ.
2. Find an estimator for this distance.
3. Choose the θ̂ that minimizes this estimated distance.

• The Problem with TV Distance: This strategy fails for TV distance because we cannot construct a good estimator for it. The TV distance TV(P_{θ*}, P_θ) does not easily decompose into an expectation with respect to P_{θ*} that can be replaced by a sample average.

4. A Surrogate Distance: Kullback-Leibler (KL) Divergence
To overcome this obstacle, a different measure of dissimilarity between distributions is introduced: the Kullback-Leibler (KL) divergence.

• Definition: For two distributions with densities p and q, the KL divergence is defined as:
  KL(p || q) = ∫ p(x) * log(p(x) / q(x)) dx
• Properties:
  • Not a true distance: KL divergence is not symmetric (KL(p || q) ≠ KL(q || p)) and does not satisfy the triangle inequality.
  • Key Property (Gibbs' Inequality): KL(p || q) ≥ 0, with equality if and only if p = q. This means that while not a distance, it still serves as a valid measure of dissimilarity, minimized when the distributions are identical.

5. From KL Divergence to Maximum Likelihood
The crucial advantage of KL divergence is that it can be estimated from data.

• KL as an Expectation: The KL divergence between the true distribution (P_{θ*}) and a candidate (P_θ) can be written as an expectation:
  KL(P_{θ*} || P_θ) = E_{θ*}[log(p_{θ*}(X) / p_θ(X))]
• Decomposition: Using the property log(a/b) = log(a) - log(b), this decomposes into two terms:
  KL(P_{θ*} || P_θ) = E_{θ*}[log(p_{θ*}(X))] - E_{θ*}[log(p_θ(X))]
• The Key Insight: We want to choose a θ that minimizes this KL divergence.
  • The first term, E_{θ*}[log(p_{θ*}(X))], is the negative entropy of the true distribution. It is an unknown constant with respect to our choice of θ.
  • Therefore, minimizing the KL divergence with respect to θ is equivalent to maximizing the second term: E_{θ*}[log(p_θ(X))].
• The Final Step (Replacing Expectation with Average): We don't know the true expectation E_{θ*}, but we have IID samples X₁, ..., Xn drawn from P_{θ*}. We can estimate the expectation by a sample average. Maximizing the estimated expectation (1/n) * Σ log(p_θ(X_i)) is equivalent to maximizing the sum without the 1/n factor.
• The Likelihood Function: The quantity Σ log(p_θ(X_i)) is known as the log-likelihood function. The principle of finding the θ that maximizes this function is the celebrated Method of Maximum Likelihood Estimation (MLE).

The lecture concludes by showing that this very common and powerful estimation method (MLE) can be rigorously derived from the fundamental goal of finding a distribution that is "close" to the true one, using KL divergence as a mathematically convenient surrogate for the more intuitive TV distance.
6. Formalizing Statistical Models
The lecture also provides a formal definition of a statistical model.

• Definition: A model is a pair (E, {P_θ}), where E is the sample space and {P_θ} is a family of distributions indexed by a parameter θ from a parameter space Θ.
• Parametric vs. Non-parametric: The course focuses on parametric models, where Θ is a finite-dimensional space (e.g., R or R²). Non-parametric models, like estimating an arbitrary density function, have an infinite-dimensional parameter space.
• Examples: The lecture provides formal definitions for several common models: Bernoulli, Exponential, Poisson, and Normal (Gaussian), specifying their sample spaces and parameter spaces.
• Identifiability: A model is identifiable if different parameter values θ always lead to different distributions P_θ. This is a crucial property, as it ensures that if we can learn the distribution, we can uniquely determine the parameter. An example is given where observing only whether a Gaussian variable is positive or negative makes the mean μ and standard deviation σ individually unidentifiable; only their ratio μ/σ can be identified.

L5: MLE cont.

The lecture begins with a recap of the key metrics used to evaluate the performance of a statistical estimator θ̂.

• Bias: Bias(θ̂) = E[θ̂] - θ. An estimator is unbiased if its bias is zero. This means it is correct on average.
• Variance: Var(θ̂). This measures the estimator's fluctuation or stability across different samples.
• Quadratic Risk (Mean Squared Error - MSE): MSE(θ̂) = E[(θ̂ - θ)²]. The primary measure of an estimator's overall quality.
• Bias-Variance Decomposition: The MSE can be decomposed into MSE(θ̂) = (Bias(θ̂))² + Var(θ̂). This highlights the fundamental trade-off: an ideal estimator needs to be both accurate on average (low bias) and stable (low variance).
• Example (Bernoulli The instructor compares three estimators for the proportion p:
  1. The sample mean X̄n: Unbiased, with variance p(1-p)/n that decreases with n.
  2. The constant 0.5: Has zero variance but is biased unless the true p is exactly 0.5.
  3. The single observation X₁: Unbiased, but has a high, constant variance p(1-p).
     A plot of their risks (MSEs) as a function of the true p shows that while the constant estimator might be better if p is very close to 0.5, the sample mean X̄n is superior for larger sample sizes n as its risk shrinks across the entire parameter space.

2. A Principled Approach to Estimation: Minimizing Distance Between Distributions
The lecture introduces a more fundamental, principled approach to finding good estimators. Instead of just estimating a parameter, the goal is to find a probability distribution P_θ̂ that is "close" to the true, unknown data-generating distribution P_{θ*}.

• Total Variation (TV) Distance: The TV distance is presented as the most natural measure of closeness, as it quantifies the maximum possible difference in probability for any event.
• The Challenge: It is difficult to construct an estimator for the TV distance itself, as it doesn't easily decompose into an expectation that can be replaced by a sample average.

3. Kullback-Leibler (KL) Divergence: A More Convenient Surrogate
To overcome this, a different measure of dissimilarity, the Kullback-Leibler (KL) divergence, is introduced as a mathematically convenient surrogate for TV distance.

• Definition: KL(p || q) = ∫ p(x) * log(p(x) / q(x)) dx.
• Key Properties:
  • It is not a true distance (it's not symmetric).
  • Crucially, KL(p || q) ≥ 0, with equality if and only if p = q (Gibbs' Inequality). This means it is a valid measure of dissimilarity, minimized when the distributions are identical.
• The Advantage: The KL divergence can be written as an expectation. KL(P_{θ*} || P_θ) = E_{θ*}[log(p_{θ*}(X) / p_θ(X))].

4. Deriving the Maximum Likelihood Principle
The lecture then shows how the strategy of minimizing the estimated KL divergence naturally leads to the principle of maximum likelihood.

1. Decompose KL: The KL divergence is decomposed: KL(P_{θ*} || P_θ) = E_{θ*}[log(p_{θ*}(X))] - E_{θ*}[log(p_θ(X))].
2. Identify the Objective: To find the best θ, we want to minimize this expression. The first term, E_{θ*}[log(p_{θ*}(X))], is the negative entropy of the true distribution. It is an unknown constant with respect to our choice of θ. Therefore, minimizing the KL divergence is equivalent to maximizing the second term: E_{θ*}[log(p_θ(X))].
3. Replace Expectation with Average: We don't know the true expectation E_{θ*}. The core statistical principle is to estimate it with a sample average over our data X₁, ..., Xn. So, we aim to maximize: (1/n) * Σ log(p_θ(X_i)).
4. Define the Likelihood: Maximizing the average of the logs is equivalent to maximizing the sum of the logs, which is equivalent to maximizing the product of the probabilities (or densities) themselves. This product, L(θ) = Π p_θ(X_i), is defined as the likelihood function.

• The Maximum Likelihood Principle: The strategy of choosing the parameter θ that maximizes the likelihood function is called Maximum Likelihood Estimation (MLE). This derivation shows that MLE is a principled method rooted in the idea of finding a distribution that is "close" to the true one.

5. Maximizing the Likelihood Function (Optimization Interlude)
The MLE principle turns an estimation problem into an optimization problem: maximizing the likelihood function.

• Concavity: A key property of many standard statistical models is that their log-likelihood functions are concave. A concave function has at most one local maximum, which is also the global maximum.
• Finding the Maximum: For a differentiable concave function, the maximum can be found by taking the derivative (or gradient in higher dimensions), setting it to zero, and solving for the parameter θ.
• Examples of Concave Functions: The lecture reviews basic calculus to identify concave functions (e.g., -x², √x, log(x)) by checking if their second derivative is non-positive.
• Boundary Solutions: It is noted that if the maximizer is not in the interior of the parameter space, the derivative-based method may not work, and the maximum may occur at the boundary.

6. Calculating MLEs for Standard Models
The lecture demonstrates how to apply the MLE principle to find estimators for common statistical models.

• Bernoulli Model:
  • The log-likelihood is derived: log L(p) = (Σ X_i) * log(p) + (n - Σ X_i) * log(1-p).
  • Taking the derivative with respect to p, setting it to zero, and solving yields the MLE: p̂_MLE = (Σ X_i) / n = X̄n. This confirms that the intuitive sample proportion is the MLE.
• Poisson Model:
  • The log-likelihood is derived: log L(λ) = (Σ X_i) * log(λ) - nλ - Σ log(X_i!).
  • Taking the derivative with respect to λ and setting to zero yields the MLE: λ̂_MLE = (Σ X_i) / n = X̄n. Again, the sample mean is the MLE for the Poisson rate parameter.
• Gaussian Model: Computing the MLE for the Gaussian parameters (μ, σ²) is left as an exercise.

7. Performance of the MLE: The Fisher Information
The lecture concludes with a preview of how the performance of the MLE is analyzed.

• The Question: How much does the MLE θ̂ fluctuate around the true θ*?
• Fisher Information: A key quantity called the Fisher Information, I(θ), is introduced. It is defined in terms of the derivatives of the log-likelihood function and measures the curvature of the log-likelihood function at its maximum.
• Intuition: A highly curved ("peaked") log-likelihood function means that the data provides a lot of information about the parameter; small changes in θ lead to large changes in likelihood. In this case, the MLE will have low variance. A flat log-likelihood implies low information and a high-variance MLE.
• Asymptotic Variance of MLE: A major result in statistics (to be discussed later) is that the asymptotic variance of the MLE is the inverse of the Fisher Information: Var(θ̂) ≈ 1 / (n * I(θ)). This provides a direct link between the shape of the likelihood function and the precision of the resulting estimator.

L6: MLE cont. and MOM

The lecture begins with a recap of how to evaluate statistical estimators using the concepts of bias, variance, and quadratic risk (Mean Squared Error - MSE).

• Bias: Bias(θ̂) = E[θ̂] - θ. An unbiased estimator is correct on average (Bias = 0).
• Variance: Var(θ̂). Measures the estimator's fluctuation.
• Quadratic Risk (MSE): MSE(θ̂) = E[(θ̂ - θ)²]. The primary measure of an estimator's overall quality.
• Bias-Variance Decomposition: A fundamental result shows that MSE(θ̂) = (Bias(θ̂))² + Var(θ̂). This highlights the trade-off in estimator design. An estimator with low MSE will be consistent (converges to the true parameter as sample size increases).
• Example (Bernoulli
  • The sample mean X̄n is shown to be unbiased with a variance of θ(1-θ)/n. Its risk decreases as n increases.
  • The constant estimator 0.5 has zero variance but is biased. Its risk (0.5 - θ)² is small only when the true θ is near 0.5.
  • The estimator X₁ (using only one observation) is unbiased but has a high, constant variance θ(1-θ).
  • A plot of the risks reveals that while a simple guess might be better for a small sample size if the true parameter is known to be in a certain range, the sample mean X̄n is a superior estimator for larger n because its risk shrinks to zero.

2. Maximum Likelihood Estimation (MLE) Review
The lecture revisits the principle of Maximum Likelihood Estimation (MLE) derived in the previous session.

• The Principle: MLE is a principled method for finding estimators. It is derived from the goal of minimizing the Kullback-Leibler (KL) divergence between the true data-generating distribution P_{θ*} and a candidate distribution P_θ from a statistical model.
• Likelihood Function: Minimizing the estimated KL divergence is equivalent to maximizing the likelihood function, L(θ) = Π p_θ(X_i), or more conveniently, its logarithm, the log-likelihood.
• The MLE: The Maximum Likelihood Estimator θ̂_MLE is the value of θ that maximizes the likelihood (or log-likelihood) function for the observed data.
• Finding the MLE: For many standard models, the log-likelihood function is concave. This is a critical property, as it guarantees that any local maximum is also the unique global maximum. Therefore, the MLE can often be found by taking the derivative (or gradient), setting it to zero, and solving for θ.

3. Fisher Information: Quantifying Curvature and Precision
A new and central concept is introduced: Fisher Information.

• Motivation: The precision of the MLE depends on the shape of the log-likelihood function around its maximum. A sharply peaked ("curvy") function implies high precision (low variance), while a flat function implies low precision (high variance).
• Definition: The Fisher Information I(θ) is a measure of this curvature. For a one-dimensional parameter, it has two equivalent definitions:
  1. I(θ) = Var[l'(θ)] (the variance of the first derivative of the log-likelihood for a single observation).
  2. I(θ) = -E[l''(θ)] (the negative expected value of the second derivative).
• Proof of Equivalence: A proof is provided showing that these two definitions are equal, relying on the fact that the integral of a probability density function is always 1, and thus its derivatives with respect to the parameter are zero.
• Interpretation: The Fisher Information quantifies the amount of "information" that the data provides about the unknown parameter θ. A high Fisher Information corresponds to a highly curved log-likelihood and a more precise estimator.

4. Asymptotic Properties of the MLE
A major theorem is presented that describes the asymptotic behavior of the Maximum Likelihood Estimator under certain regularity conditions.

• The Theorem: Under suitable conditions (identifiability, smoothness, etc.), the MLE θ̂ has two key properties as the sample size n → ∞:
  • Consistency: θ̂ converges in probability to the true parameter θ*.
  • Asymptotic Normality: The distribution of the standardized estimator converges to a standard normal distribution:
    √n * (θ̂ - θ*) → N(0, I(θ*)⁻¹)
• Significance: This is a powerful, general result.
  • It guarantees that the MLE is a good estimator for large samples (it's consistent and unbiased in the limit).
  • It provides a CLT-like result even when the MLE is not a simple average.
  • It reveals that the asymptotic variance of the MLE is 1 / (n * I(θ*)). This directly links the estimator's precision (variance) to the Fisher Information (curvature of the likelihood), confirming the intuition developed earlier.

5. The Method of Moments: An Alternative Estimation Principle
A second major estimation technique is introduced: the Method of Moments (MoM).

• Motivation: The method is motivated by the Weierstrass Approximation Theorem, which states that any continuous function can be approximated by polynomials. This suggests that the moments of a distribution (which are expectations of powers of the random variable, E[X^k]) contain a great deal of information about the distribution.
• The Principle: The Method of Moments is based on the familiar statistical trick of "replacing expectations with averages."
  1. Population Moments: For a given model, express the first d population moments m_k = E[X^k] as functions of the d-dimensional parameter θ. This gives d equations: m_k = ψ_k(θ).
  2. Empirical Moments: Calculate the first d empirical (sample) moments from the data: m̂_k = (1/n) * Σ X_i^k.
  3. Equate and Solve: Set the population moments equal to the empirical moments (m_k = m̂_k) and solve the resulting system of d equations for the d unknown components of θ.
• The Estimator: The solution to this system is the Method of Moments Estimator, θ̂_MoM.

6. Properties of the Method of Moments Estimator

• Consistency: The MoM estimator is consistent. This is a direct consequence of the Law of Large Numbers (which ensures that empirical moments converge to population moments) and the Continuous Mapping Theorem (which ensures that applying the inverse function ψ⁻¹ preserves convergence).
• Asymptotic Normality: The MoM estimator is also asymptotically normal. This can be shown using the multivariate Central Limit Theorem (which provides a limit for the vector of empirical moments) combined with the multivariate Delta Method (which extends the Delta Method to vector-valued functions ψ⁻¹).
• Performance Comparison: While often simpler to compute, the Method of Moments estimator is typically less accurate (has a larger asymptotic variance) than the Maximum Likelihood Estimator. MLE is generally considered the superior method when it is computationally feasible. However, MoM can be a useful starting point for iterative algorithms designed to find the MLE, especially when the likelihood function is not concave.

L7: Parametric Hypothesis Testing

The lecture begins by revisiting the key properties used to evaluate a statistical estimator θ̂.

• Bias: Bias(θ̂) = E[θ̂] - θ. An unbiased estimator is correct on average (Bias = 0).
• Variance: Var(θ̂). Measures the estimator's fluctuation or stability.
• Quadratic Risk (Mean Squared Error - MSE): MSE(θ̂) = E[(θ̂ - θ)²]. The primary measure of an estimator's overall quality.
• Bias-Variance Decomposition: MSE(θ̂) = (Bias(θ̂))² + Var(θ̂). This highlights the fundamental trade-off in estimator design. An estimator with low MSE is consistent (converges to the true parameter as sample size increases).

The instructor uses the example of estimating a Bernoulli parameter p to illustrate these concepts, showing that while the sample mean X̄n is unbiased and its variance decreases with n, a naive guess like 0.5 has zero variance but is biased. The choice between estimators involves a trade-off, but X̄n is generally superior as n grows.
2. Introduction to Hypothesis Testing
The lecture transitions from estimation (finding a value for a parameter) to hypothesis testing (making a yes/no decision about a parameter).

• The Goal: Hypothesis testing is a formal procedure for making decisions based on data in the face of uncertainty. It is the foundation of data-driven scientific inquiry.
• Example: Cherry Blossom Run: The central example is introduced: are runners in the 2012 Cherry Blossom Run faster than they were in 2009?
  • In 2009, the average time was 103.5 minutes, and the distribution was approximately Gaussian.
  • The scientific question "are runners faster?" is translated into a precise statistical hypothesis about the unknown mean running time μ for 2012: Is μ < 103.5?
• The Role of Modeling Assumptions: To answer this question, we must make simplifying assumptions. We assume the 2012 running times are also Gaussian and that the variance is the same as in 2009. This allows us to focus the statistical test solely on the parameter of interest, μ.

3. The Heuristic Approach and Its Flaw
A naive heuristic for testing the hypothesis would be to take a sample of n runners from 2012, compute their average time X̄n, and conclude that runners are faster if X̄n < 103.5.

• The Flaw: This approach is flawed because the sample mean X̄n is a random variable and is almost never exactly equal to the true mean μ. Due to random sampling fluctuations, X̄n could be greater than 103.5 even if the true μ is less than 103.5 (and vice-versa).
• The Solution: We need a "buffer zone." We should only conclude that μ < 103.5 if X̄n is significantly smaller than 103.5, by an amount that accounts for the statistical noise, which is typically on the order of 1/√n. The CLT will be used to quantify this buffer.

4. The Formal Framework of Hypothesis Testing
The lecture formalizes the language and structure of hypothesis testing.

• Hypotheses: A testing problem is defined by two disjoint sets of parameters, the null hypothesis ( and the alternative hypothesis (.
  • (The Null): Represents the "status quo" or a default position of no effect. In the running example, H₀: μ = 103.5.
  • (The Alternative): Represents the scientific discovery or the effect we are looking for. In the example, H₁: μ < 103.5.
• The Asymmetric Roles: H₀ and H₁ are not treated symmetrically. The framework is designed to be conservative, siding with the null hypothesis unless there is strong evidence against it. The analogy of "innocent until proven guilty" is used: H₀ is assumed true unless the data provides overwhelming evidence to reject in favor of H₁. One never "accepts" H₀, but rather "fails to reject" it.

5. A Test: The Decision Rule
A test is a formal rule for making a decision based on the data.

• Definition: A test ψ is a statistic (a function of the data X₁, ..., Xn) that outputs either 0 (fail to reject H₀) or 1 (reject H₀).
• Test Statistic: A test is typically constructed by first computing a test statistic T(X₁, ..., Xn), which is a single number that summarizes the evidence in the data. Then, a decision is made based on whether this statistic falls into a pre-defined rejection region.
• Example (Fair Coin Test):
  • To test if a coin is fair (H₀: p = 0.5 vs. H₁: p ≠ 0.5), we use the standardized sample mean as the test statistic: T = √n * (X̄n - 0.5) / √(0.5 * 0.5).
  • Under H₀, the CLT tells us that T is approximately N(0, 1).
  • We reject H₀ if the observed value of T is "unusual" for a standard normal variable, i.e., if |T| is greater than some critical value c (e.g., c = 1.96).

6. Types of Errors and the Level of a Test
Since the decision is based on random data, it can be wrong. There are two types of errors.

• Type I Error: Rejecting H₀ when it is actually true (e.g., concluding a fair coin is biased). The probability of a Type I error is denoted by α.
• Type II Error: Failing to reject H₀ when it is false (e.g., failing to detect that a biased coin is biased). The probability of a Type II error is denoted by β.
• Power: The probability of correctly rejecting H₀ when it is false (1 - β).

7. The Neyman-Pearson Paradigm: Controlling the Type I Error
There is an inherent trade-off between Type I and Type II errors. A test that never rejects H₀ has α = 0 but maximum β. A test that always rejects H₀ has β = 0 but maximum α.

• The Neyman-Pearson Paradigm: The standard approach in statistics is to fix the probability of the more "costly" error and then optimize for the other.
• Controlling the Type I Error: By convention, H₀ (the status quo) is given the benefit of the doubt. The Type I error (wrongly rejecting the status quo) is considered more severe.
• The Procedure:
  1. Fix a maximum allowable probability for the Type I error, called the level of the test, denoted α (commonly α = 0.05).
  2. Among all possible tests that have a Type I error probability of at most α, find the one that has the highest power (or, equivalently, the lowest Type II error probability).

This paradigm provides a principled way to choose the critical value c for a test. We choose c such that the probability of the test statistic exceeding c when is exactly α. For the coin example, c = 1.96 is chosen because the probability of a standard normal variable exceeding 1.96 in absolute value is 0.05.

L8: Parametric Hypothesis Testing cont.

The lecture begins with a review of the concepts from the previous session.

• Estimator Evaluation: The quality of an estimator θ̂ is assessed by its bias (accuracy on average) and variance (stability). The overall performance is captured by the quadratic risk (Mean Squared Error), which decomposes into (Bias)² + Variance. An ideal estimator has low bias and low variance.
• Hypothesis Testing Framework: The core setup of a hypothesis testing problem is revisited:
  • The parameter space Θ is partitioned into two disjoint subsets, Θ₀ and Θ₁.
  • This defines the null hypothesis ( and the alternative hypothesis (.
  • H₀ represents the "status quo," while H₁ represents a scientific discovery or effect.
• Tests and Errors:
  • A test ( is a function of the data that outputs a binary decision: 0 (fail to reject H₀) or 1 (reject H₀).
  • Type I Error: Rejecting H₀ when it is true.
  • Type II Error: Failing to reject H₀ when H₁ is true.
  • The probabilities of these errors, α(θ) and β(θ), are functions of the true, unknown parameter θ.

2. The Power Function, Level, and Neyman-Pearson Paradigm
The lecture further develops the framework for evaluating and designing tests.

• Power Function: The power function of a test ψ is the probability of rejecting H₀, viewed as a function of the true parameter θ.
  1. P_θ(ψ=1)
• Interpretation:
  1. When θ ∈ Θ₀, the power function gives the Type I error probability.
  2. When θ ∈ Θ₁, the power function gives the probability of making a correct decision (1 - β(θ)).
• Level and Power of a Test:
  1. Level ( The maximum Type I error probability over all parameters in the null hypothesis (α = sup_{θ ∈ Θ₀} P_θ(ψ=1)). This is the worst-case probability of a false discovery.
  2. Power ( The minimum probability of a correct discovery over all parameters in the alternative hypothesis (π = inf_{θ ∈ Θ₁} P_θ(ψ=1)). This is the worst-case power.
• The Neyman-Pearson Paradigm: This is the standard principle for designing statistical tests. It addresses the inherent trade-off between Type I and Type II errors.
  1. Fix the Level: First, an acceptable maximum rate for the Type I error is chosen (e.g., α = 0.05). This constrains the test.
  2. Maximize the Power: Among all tests that satisfy the level constraint, choose the one that has the highest possible power.

3. Constructing Tests and the p-value
The lecture details the practical construction of a test and introduces the concept of the p-value.

• Test Statistic: A test is typically built around a test statistic, T(X₁, ..., Xn), a single value computed from the data that summarizes the evidence.
• Rejection Region: The test rejects H₀ if the test statistic T falls into a pre-defined rejection region. For a simple test, this region might be |T| > c for some critical value c.
• Finding the Critical Value ( The critical value c is chosen to satisfy the level α. Specifically, we find c such that P_{H₀}(|T| > c) = α. To do this, we need to know the distribution of the test statistic T under the null hypothesis. Often, this is an asymptotic distribution given by the Central Limit Theorem.
• The p-value: Instead of just a "reject" or "fail to reject" decision, one can report a p-value.
  • Definition: The p-value is the smallest level for the given data.
  • Calculation: It is calculated as the probability, under the null hypothesis, of observing a test statistic value at least as extreme as the one actually observed. For the example test, p-value = P_{H₀}(|Z| > |T_obs|), where Z ~ N(0,1) and T_obs is the observed value of the test statistic.
  • Interpretation: A small p-value (e.g., < 0.05) indicates that the observed data is very unlikely if the null hypothesis were true, providing strong evidence to reject H₀. It provides a more nuanced measure of evidence than a simple binary decision.

4. The Chi-Squared (χ²) Distribution
To move beyond tests about means and proportions, which rely on the Gaussian distribution, the lecture introduces a new pivotal distribution: the chi-squared distribution. This is essential for making inferences about variances.

• Definition: The χ² distribution with d degrees of freedom (χ²_d) is defined as the distribution of the sum of d independent, squared standard normal random variables:
  Z₁² + Z₂² + ... + Z_d² ~ χ²_d, where Z_i ~ N(0,1).
• Properties:
  • It is a distribution for a non-negative random variable.
  • Its mean is d and its variance is 2d.
  • Its shape changes with d, gradually shifting to the right and becoming more symmetric as d increases.

5. Cochran's Theorem and Inference on Variance
The connection between the χ² distribution and statistical inference is made through Cochran's Theorem.

• The Sample Variance: The sample variance is defined as S_n² = (1/(n-1)) * Σ(X_i - X̄n)².
• Cochran's Theorem: If the observations X₁, ..., Xn are IID N(μ, σ²), then the scaled sample variance has an exact chi-squared distribution:
  (n-1) * S_n² / σ² ~ χ²_{n-1}
• Significance:
  • This provides an exact, non-asymptotic pivotal quantity for the variance σ².
  • The degrees of freedom are n-1, not n. Intuitively, one degree of freedom is "used up" by estimating the mean μ with the sample mean X̄n.
  • This result allows us to construct exact confidence intervals and hypothesis tests for the variance σ² when the underlying data is known to be Gaussian.

The lecture concludes by noting that other important distributions, like the Student's t-distribution, are built from the Gaussian and chi-squared distributions, which will be covered in the next session.

L9: Parametric Hypothesis Testing cont.

The lecture begins by reviewing the core concepts and vocabulary of hypothesis testing from the previous session.

• Hypotheses: A testing problem is defined by partitioning the parameter space Θ into two disjoint subsets, Θ₀ and Θ₁, corresponding to the null hypothesis ( and the alternative hypothesis (. H₀ represents the "status quo," while H₁ represents the "discovery."
• A Test ( A test is a statistic that maps the data to a binary decision, 0 (fail to reject H₀) or 1 (reject H₀). The set of data points that lead to rejection is called the rejection region.
• Errors:
  • Type I Error: Rejecting H₀ when it is true.
  • Type II Error: Failing to reject H₀ when H₁ is true.
• Performance Metrics:
  • Level ( The maximum probability of a Type I error over all parameters in H₀.
  • Power ( The minimum probability of correctly rejecting H₀ over all parameters in H₁.

2. The Student's t-distribution
A new pivotal distribution, the Student's t-distribution, is introduced. It is essential for conducting inference on the mean when the variance is unknown.

• Motivation: When performing a test or building a confidence interval for a Gaussian mean μ, the standardized quantity √n(X̄n - μ)/σ follows a standard normal distribution. However, this requires knowing the true standard deviation σ. In practice, σ is usually unknown and must be estimated from the data using the sample standard deviation S_n.
• The Problem: Replacing the true σ with its estimate S_n introduces additional randomness. The resulting statistic, √n(X̄n - μ)/S_n, no longer follows an exact normal distribution, especially for small sample sizes.
• Definition: The t-distribution with is defined as the distribution of a standard normal random variable Z divided by the square root of an independent χ²_d random variable scaled by its degrees of freedom: Z / √(V/d).
• Application (Cochran's Theorem): Cochran's Theorem states that for IID Gaussian data, the sample mean X̄n is independent of the sample variance S_n², and that (n-1)S_n² / σ² follows a χ²_{n-1} distribution. Combining these facts shows that the statistic √n(X̄n - μ)/S_n (using the unbiased version of S_n² with n-1 in the denominator) follows a t-distribution with .
• Properties: The t-distribution resembles the standard normal distribution but has heavier tails. This means it accounts for the extra uncertainty introduced by estimating the variance. As the degrees of freedom (n-1) increase, the t-distribution converges to the standard normal distribution. For n > 30 or 50, the difference is often negligible.

3. Wald's Test: A General Test Based on the MLE
The lecture introduces a general and powerful hypothesis test based on the Maximum Likelihood Estimator (MLE), known as Wald's Test.

• The Goal: To test a simple hypothesis H₀: θ = θ₀ versus H₁: θ ≠ θ₀ for a d-dimensional parameter θ.
• The Idea: The test is based on measuring the distance between the MLE θ̂ and the hypothesized value θ₀. If this distance is "large," we reject H₀.
• The Test Statistic: A naive Euclidean distance ||θ̂ - θ₀||² is not ideal because it doesn't account for the fact that the MLE has different variances and covariances in different directions. The Wald test statistic correctly normalizes this distance using the Fisher Information matrix, I(θ), which acts as a "metric" that captures the geometry of the parameter space. The test statistic is:
  T_n = n * (θ̂ - θ₀)ᵀ * I(θ₀) * (θ̂ - θ₀)
  (In practice, the unknown I(θ₀) is often replaced by I(θ̂) using Slutsky's Theorem).
• Asymptotic Distribution: A key theorem states that under H₀, the Wald test statistic T_n converges in distribution to a distribution with (χ²_d), where d is the dimension of the parameter θ.
• The Test: The Wald test rejects H₀ at level α if T_n is greater than the (1-α)-quantile of the χ²_d distribution.

4. The Likelihood Ratio Test (LRT)
Another general and fundamental hypothesis test is introduced: the Likelihood Ratio Test (LRT).

• The Idea: The LRT is based on a very intuitive principle: compare the maximized likelihood under the null hypothesis to the maximized likelihood under the alternative (or the entire space). If the likelihood under the null is significantly smaller, it suggests the null is a poor fit for the data.
• Simple vs. Composite Hypotheses:
  • For a simple hypothesis vs. a simple alternative (H₀: θ = θ₀ vs. H₁: θ = θ₁), the LRT rejects H₀ if the ratio L(θ₁)/L(θ₀) is large. The Neyman-Pearson Lemma proves this is the most powerful test possible.
  • For more general composite hypotheses (where Θ₀ or Θ₁ contain more than one point), we compare the maximum likelihood achievable within Θ₀ to the maximum likelihood achievable over the entire parameter space Θ.
• The Test Statistic:
  • Let θ̂_MLE be the standard MLE, maximizing L(θ) over all of Θ.
  • Let θ̂_c be the constrained MLE, maximizing L(θ) only over θ ∈ Θ₀.
  • The LRT statistic is T_n = 2 * (log L(θ̂_MLE) - log L(θ̂_c)).
• Asymptotic Distribution (Wilks' Theorem): Under H₀, the LRT statistic T_n converges in distribution to a distribution. The degrees of freedom is the difference in the number of free parameters between the full model (Θ) and the null model (Θ₀). For instance, if Θ is d-dimensional and H₀ constrains d-r of those parameters, the degrees of freedom is d-r.
• Relationship to Wald's Test: For many standard models (like the Gaussian model), the Wald test and the LRT are asymptotically equivalent.

5. The Delta Method for Hypothesis Testing
The lecture concludes by showing how the Delta Method can be used to construct tests for functions of parameters.

• The Goal: Test a hypothesis of the form H₀: g(θ) = 0, where g is a function mapping the d-dimensional parameter θ to a k-dimensional space.
• The Procedure:
  1. Start with the asymptotic normality of the MLE: √n(θ̂ - θ) → N(0, I(θ)⁻¹).
  2. Apply the multivariate Delta Method to find the asymptotic distribution of √n(g(θ̂) - g(θ)). The new asymptotic covariance matrix will involve the gradient of the function g.
  3. Under H₀, g(θ) = 0. The resulting quantity √n * g(θ̂) will be asymptotically normal with mean 0.
  4. A Wald-type test statistic can be formed by taking the squared Mahalanobis distance of √n * g(θ̂). This statistic will have an asymptotic distribution with , where k is the dimension of the output of g.

This provides a general recipe for constructing tests for a wide variety of complex hypotheses involving functions of the underlying model parameters.

L10: Testing Goodness of Fit

The lecture begins with a review of the hypothesis testing framework and the Chi-Squared (χ²) test, which was introduced in the previous session.

• Hypothesis Testing Jargon: The key terms are revisited: null hypothesis (H₀), alternative hypothesis (H₁), a test (ψ), rejection region, Type I error (rejecting a true H₀), and Type II error (failing to reject a false H₀).
• The Power Function: The concept of the power function, P_θ(ψ=1), is reintroduced. It represents the probability of rejecting H₀ as a function of the true parameter θ.
  • When θ is in the null set Θ₀, this function gives the Type I error probability.
  • When θ is in the alternative set Θ₁, this function gives the power of the test (the probability of a correct discovery).
• Level and Power:
  • Level ( The maximum Type I error probability over Θ₀.
  • Power ( The minimum power over Θ₁.
• Neyman-Pearson Paradigm: The standard approach is to fix an acceptable level α (e.g., 0.05) and then find a test that maximizes the power.
• Chi-Squared Test Review: The χ² test is a goodness-of-fit test for discrete (multinomial) data. It tests whether the observed frequencies of categories match a hypothesized set of probabilities (p₀). Its test statistic is shown to have an asymptotic χ² distribution, which is used to determine the rejection region.

2. The Student's t-test: Inference on Means with Unknown Variance
The lecture introduces a new, fundamental test: the Student's t-test.

• Motivation: When testing a hypothesis about the mean μ of a Gaussian distribution (H₀: μ = μ₀), the standard pivotal quantity √n(X̄n - μ₀)/σ requires knowledge of the true standard deviation σ. In most real-world scenarios, σ is unknown.
• The Problem: If we replace the true σ with its sample estimate S_n, the resulting statistic no longer follows a normal distribution, especially for small sample sizes (n). The uncertainty in estimating σ adds extra variability.
• The Solution: The t-distribution: The Student's t-distribution is the exact distribution of the statistic √n(X̄n - μ)/S_n (using the unbiased S_n with n-1 in the denominator) when the underlying data is Gaussian.
  • Definition: A t-distribution with d degrees of freedom is defined as the ratio of a standard normal random variable to the square root of an independent χ²_d variable scaled by d.
  • Connection: Cochran's theorem establishes that for Gaussian data, X̄n is independent of S_n² and that a scaled version of S_n² follows a χ²_{n-1} distribution, proving that the test statistic follows a t_{n-1} distribution.
• Properties: The t-distribution resembles a standard normal but has heavier tails, accounting for the extra uncertainty from estimating the variance. As the sample size n (and thus the degrees of freedom n-1) increases, the t-distribution converges to the normal distribution. For n > 30 or 50, the difference is often negligible.
• The t-test: The procedure is analogous to the Gaussian-based test (z-test), but the critical values are taken from the appropriate t-distribution table instead of the normal distribution table.

3. Goodness-of-Fit Tests for Continuous Distributions
The lecture then moves to the more general problem of testing whether a sample comes from a specific continuous distribution, a problem known as goodness-of-fit.

• The Challenge: Unlike the discrete multinomial case, we cannot simply compare frequencies in bins. We need a way to compare the entire shape of the observed data distribution to a hypothesized one.
• The Empirical CDF: The key tool is the Empirical Cumulative Distribution Function (ECDF), denoted F_n(t).
  • Definition: F_n(t) is the proportion of data points in the sample that are less than or equal to t.
  • F_n(t) = (1/n) * Σ I(X_i ≤ t)
• Consistency (Glivenko-Cantelli Theorem): The ECDF is a very good estimator of the true CDF, F(t). The Glivenko-Cantelli Theorem, a powerful result, states that the ECDF converges uniformly to the true CDF as n → ∞. That is, the maximum difference between the two functions over all t goes to zero.

4. The Kolmogorov-Smirnov (KS) Test
The Glivenko-Cantelli theorem motivates a classic goodness-of-fit test.

• The Idea: If the null hypothesis H₀: F = F₀ is true, then the empirical CDF F_n should be close to the hypothesized CDF F₀ everywhere. The test is based on measuring the maximum vertical distance between these two functions.
• The Test Statistic: The Kolmogorov-Smirnov (KS) statistic is T_n = sup_t |√n (F_n(t) - F₀(t))|.
• Asymptotic Distribution: A remarkable theorem by Donsker shows that, under H₀, the asymptotic distribution of this statistic is the supremum of a Brownian bridge.
  • A Brownian bridge is a specific stochastic process (a random function) that starts at 0, ends at 0, and behaves like a random walk in between.
• The Pivotal Property: The crucial property is that this limiting distribution is pivotal—it does not depend on the specific hypothesized distribution F₀ (as long as F₀ is continuous). This means a single table of critical values can be used to test for any continuous distribution.
• The KS Test: The test rejects H₀ if the observed KS statistic T_n is larger than the appropriate quantile (q_α) read from the table for the supremum of a Brownian bridge.

5. Practical Considerations for Goodness-of-Fit

• Computing the KS Statistic: The supremum over all real numbers t seems computationally intractable. However, the ECDF is a step function that only changes at the observed data points. It can be shown that the maximum difference between F_n and F₀ must occur at one of the data points X_i. Therefore, the supremum can be computed efficiently by checking only 2n values.
• Testing for Normality with Unknown Parameters: The standard KS test requires the null distribution F₀ to be fully specified (e.g., N(0, 5)). If we want to test for normality in general (i.e., is the data from any normal distribution?), we must first estimate the mean μ and variance σ² from the data and plug them into F₀.
  • The Problem: This "plug-in" procedure makes the null hypothesis "fit" the data better, which violates the assumptions of the standard KS test. Using the standard KS critical values in this case will lead to a test that is too conservative (it will fail to reject H₀ too often).
  • The Solution (Lilliefors Test): A different set of critical values, derived from the Lilliefors test, must be used in this situation. These values are smaller, making the test more powerful and appropriately accounting for the estimation of the parameters.

The lecture concludes by mentioning other goodness-of-fit tests, like the Cramér-von Mises and Anderson-Darling tests, which use different norms (e.g., L2-norm) to measure the distance between the ECDF and the hypothesized CDF.

L11: Testing Goodness of Fit cont.

The lecture begins by introducing goodness-of-fit tests as a way to answer the question: "Does my data come from a particular distribution?" This is a crucial step in many statistical analyses. For example, one might test if the zodiac signs of Fortune 500 CEOs are uniformly distributed, or more practically, if a dataset is normally distributed, which is a key assumption for applying a t-test.
Several standard tests exist for this purpose:

• Kolmogorov-Smirnov (K-S) Test: The most widely used test, it measures the greatest distance between the empirical cumulative distribution function (CDF) of the data and the CDF of the hypothesized distribution. Its popularity stems from its intuitive nature.
• Cramér-von Mises Test: This test uses an L2 distance (related to the average squared difference) between the empirical and theoretical CDFs.
• Anderson-Darling Test: A more complex test that gives more weight to the tails of the distribution.

The choice of test can be subjective, as each has advantages and may be more powerful against different types of deviations from the hypothesized distribution.
A Special Case: Testing for Normality
A primary application of goodness-of-fit tests is to check if data is normally distributed before applying methods like the Student's t-test. However, "normal" isn't a single distribution but a family of distributions defined by a mean (µ) and a variance (σ²).
When testing for normality, the specific µ and σ² are usually unknown. The natural approach is to estimate them from the data itself using the sample mean (μ̂) and sample variance (σ̂²). The test then checks if the data fits a normal distribution with these estimated parameters.
This estimation process introduces a complication:

• The standard theory behind the K-S test (specifically, Donsker's theorem) no longer applies.
• The distribution of the K-S test statistic changes because μ̂ and σ̂² are random variables, not fixed constants. The quantiles, or critical values for the test, are distorted as a result.

The Lilliefors Test: A Modification for Normality
To address the issue of testing for normality with estimated parameters, a modification of the K-S test known as the Lilliefors test is used.
The key insight is that while the test statistic's distribution is no longer the standard K-S distribution, it is pivotal. This means that under the null hypothesis (that the data is indeed Gaussian), its distribution does not depend on the specific true values of µ and σ.
This property arises from standardizing the data. The test statistic, which depends on terms like (μ̂ - μ) / σ and σ̂² / σ², relies on quantities whose distributions (Normal and chi-squared, respectively) do not depend on the unknown parameters. Because these quantities are also independent (a result of Cochran's theorem for Gaussian data), their joint distribution is also independent of μ and σ.
Interpreting Lilliefors Test Results:

• Since the test is pivotal, new tables of critical values can be calculated specifically for this scenario.
• These critical values are significantly smaller than those for the standard K-S test. For instance, a value of 0.4 might be used for the K-S test, whereas the Lilliefors test might use a value closer to 0.2.
• This means the Lilliefors test is more powerful or "harder not to reject." For the same data, one might reject the null hypothesis of normality with the Lilliefors test while failing to reject it with the standard K-S test.

Visual Diagnostics: Quantile-Quantile (QQ) Plots
Beyond numerical tests, visual tools are essential for diagnosing how a distribution deviates from a theoretical one. Quantile-Quantile (QQ) plots are a powerful method for this, especially for examining the tails of a distribution.
How QQ Plots Work:

• A QQ plot is a scatter plot comparing the empirical quantiles of the data against the theoretical quantiles of the distribution being tested (e.g., a normal distribution).
• The empirical quantiles are simply the sorted data points.
• If the data perfectly follows the theoretical distribution, the points on the QQ plot will form a straight line.

Interpreting Deviations:

• Heavy Tails: If the data has heavier tails than the normal distribution, the points at both ends of the plot will stray from the reference line. The points on the right (high values) will curve upwards, and the points on the left (low values) will dip downwards. This creates a characteristic elongated "S" shape. Student's t-distribution is an example of a heavy-tailed distribution.
• Light Tails: If the data has lighter tails, the pattern is reversed. The points on the right will dip downwards, and the points on the left will curve upwards, creating a compressed "S" shape. A uniform distribution, which has no tails, exhibits this pattern.
• Skewness: If one tail is heavier or lighter than the other, the plot will show an asymmetric curve. An exponential distribution, for example, has no left tail but a heavier right tail than a Gaussian, resulting in a "U" shape on the plot.

The Chi-Square Test for Parametric Families
The lecture revisits the chi-square goodness-of-fit test, extending it from testing against fixed probabilities (like a fair die) to testing against a parametric family of distributions, such as the Binomial or Poisson distributions.
The Process:

1. Binning: The data, whether discrete or continuous, is grouped into a finite number of categories or bins.
2. Parameter Estimation: The unknown parameters of the distribution (e.g., the probability p in a binomial distribution) are estimated from the data, typically using the maximum likelihood estimator.
3. Calculating Probabilities: The estimated parameters are used to calculate the expected probability for each bin according to the theoretical distribution.
4. Test Statistic: The chi-square statistic is computed by summing the squared differences between the observed proportions in each bin and the expected probabilities, scaled by the expected probabilities.

Degrees of Freedom:

• A crucial change is that the degrees of freedom for the test are reduced.
• The formula becomes K - 1 - d, where K is the number of bins and d is the number of parameters estimated from the data.
• Estimating parameters from the data makes the theoretical model fit the observed data more closely, thus "using up" degrees of freedom.

This method allows the chi-square test to be a flexible tool for testing goodness-of-fit for a wide range of distributions, including continuous ones like the normal distribution, by first discretizing the data through binning.

L12: Regression

Regression analysis aims to predict or explain a variable, called the explained variable (denoted as y), as a function of one or more other variables, known as explanatory variables (denoted as x). The core task is to estimate the underlying function f in the relationship y = f(x) using a set of data points (xi, yi).
2. The Linear Regression Model
The simplest and most common form of regression is linear regression, which assumes the relationship between x and y is a straight line.

• Model Equation: Since data points rarely fall perfectly on a line, the model includes a noise term, ε (epsilon), to account for deviations and unobserved factors. The model is expressed as:
  y = a + bx + ε
  • a is the intercept (the value of y when x is zero).
  • b is the slope (the change in y for a one-unit change in x).
  • ε is the random noise or error term.
• Univariate vs. Multivariate:
  • Univariate Regression: Involves one explanatory variable (x is a single number).
  • Multivariate Regression: Involves multiple explanatory variables (x is a vector (x1, x2, ..., xp)), and the model takes the form y = a + b1*x1 + b2*x2 + ... + bp*xp + ε.
• Transformations: Even if a relationship isn't inherently linear, it can often be transformed into one. For example, the Ideal Gas Law PV = nRT can be made linear by taking the logarithm of both sides, resulting in a linear relationship between log(P), log(V), and log(T).

3. Finding the "Best Fit" Line: The Method of Least Squares
To find the best possible line that explains the data, a criterion is needed to measure how well a given line fits the points. The standard method is least squares.

• The Principle: The method of least squares aims to find the values of the parameters (a and b) that minimize the sum of the squared vertical distances between each data point yi and the line a + bxi.
  • Mathematically, this means minimizing the expression: Σ(yi - (a + bxi))^2
• Theoretical Solution (Population Level): If the true underlying distribution of x and y were known, the optimal a and b that minimize the expected squared error E[(y - (a + bx))^2] are:
  • Slope b = Cov(x, y) / Var(x)
  • Intercept a = E[y] - b * E[x]
    This shows that the slope is fundamentally determined by the covariance between x and y, scaled by the variance of x.
• Practical Solution (Sample Level): In practice, we don't know the true expectations and covariances. We estimate them from our data sample. The estimators â and b̂ are found by replacing the theoretical values with their sample averages in the formulas above.

4. Properties of the Noise Term (ε)
The least squares method implicitly makes two important assumptions about the noise term ε:

1. Zero Mean: The expected value of the noise is zero (E[ε] = 0). This means the errors are not systematically positive or negative.
2. Uncorrelated with x: The noise is uncorrelated with the explanatory variable (Cov(x, ε) = 0). This ensures that the errors are random and not related to the value of x.
   More advanced statistical inference often requires stronger assumptions, such as the noise being normally distributed.

5. Multivariate Regression and Matrix Notation
When dealing with multiple explanatory variables, solving for each coefficient individually becomes complicated. The problem is simplified by using matrix notation.

• The Model in Matrix Form: All n observations can be represented in a single, compact equation:
  y = Xβ + ε
  • y is a vector of the observed outcomes.
  • X is the design matrix, where each row represents an observation and each column represents an explanatory variable (with an added column of ones for the intercept).
  • β (beta) is the vector of coefficients we want to estimate.
  • ε is the vector of noise terms.
• Least Squares in Matrix Form: The goal of finding the best-fit "plane" (in higher dimensions) is equivalent to finding the vector β̂ that minimizes the squared length (norm) of the residual vector:
  Minimize ||y - Xβ||^2

This matrix formulation is the standard approach for solving multivariate linear regression problems.

L13: Regression cont.

The lecture begins by formalizing the multivariate linear regression model using matrix notation, which simplifies the calculations.

• Model Setup: The n observations are stacked into vectors and a matrix:
  • y: An n x 1 vector of the observed outcomes (y₁ to yₙ).
  • X: An n x p design matrix, where each row is an observation (xᵢᵀ) and each column represents an explanatory variable. A column of ones is typically included to account for the intercept.
  • β: A p x 1 vector of the unknown coefficients to be estimated.
  • ε: An n x 1 vector of the random noise terms (ε₁ to εₙ).
    The model for all observations is concisely written as: y = Xβ + ε
• The Least Squares Criterion: The goal is to find the estimator β̂ that minimizes the sum of squared residuals. In matrix form, this is equivalent to minimizing the squared Euclidean norm of the residual vector:
  Minimize
• Derivation of the Estimator: By taking the gradient of this expression with respect to β, setting it to zero, and solving the resulting system of linear equations (known as the score equations), we arrive at the formula for the least squares estimator:
  β̂ = (XᵀX)⁻¹Xᵀy

2. Conditions for the Estimator to Exist
The formula for β̂ requires the matrix (XᵀX) to be invertible.

• The Invertibility Condition: For (XᵀX) to be invertible, the design matrix X must have full column rank, which means rank(X) = p. This implies that none of the explanatory variables can be a perfect linear combination of the others.
• Sample Size vs. Parameters ( A direct consequence is that the number of observations (n) must be greater than or equal to the number of parameters (p). If p > n, there are more parameters to estimate than data points, leading to an infinite number of solutions. This is the central challenge of high-dimensional statistics, where modern datasets (like in genomics) often have p >> n.

3. The Geometric Interpretation: Orthogonal Projection
Linear regression has a powerful geometric interpretation.

• Predicted Values ( The vector of fitted or predicted values is given by ŷ = Xβ̂.
• Projection: Substituting the formula for β̂, we get ŷ = X(XᵀX)⁻¹Xᵀy. The matrix P = X(XᵀX)⁻¹Xᵀ is a projection matrix. It performs an orthogonal projection of the data vector y onto the vector space spanned by the columns of the design matrix X (the column space).
• Pythagoras' Theorem: Because the projection is orthogonal, the total variation in y can be perfectly partitioned into the variation explained by the model and the residual variation:
  ||y||² = ||ŷ||² + ||y - ŷ||²
  (Total Sum of Squares = Explained Sum of Squares + Residual Sum of Squares)

4. Statistical Assumptions for Inference
To move from merely fitting a line to performing statistical inference (e.g., creating confidence intervals or hypothesis tests), we need a formal statistical model with three key assumptions:

1. Deterministic Design Matrix The values of the explanatory variables are treated as fixed, not random.
2. Homoscedasticity: The noise terms εᵢ are independent and have the same constant variance, σ².
3. Gaussian Noise: The noise vector ε follows a multivariate normal distribution with a mean of zero and a covariance matrix of σ²I (where I is the identity matrix).

5. The Distribution of the Estimator and Residuals
Under these assumptions, we can derive the exact distributions of our estimators.

• Distribution of The least squares estimator β̂ is a linear transformation of the Gaussian vector y, so it is also normally distributed:
  β̂ ~ N(β, σ²(XᵀX)⁻¹)
  This result is fundamental. It shows that β̂ is an unbiased estimator of β (its mean is the true β), and its variance depends on the noise level σ² and the geometry of the design matrix X.
• Least Squares as Maximum Likelihood: Under the Gaussian noise assumption, the least squares estimator is identical to the Maximum Likelihood Estimator (MLE). This is because maximizing the Gaussian log-likelihood function is equivalent to minimizing the sum of squared errors.
• Estimating the Noise Variance The unknown noise variance σ² can be estimated from the residuals (y - ŷ). The normalized sum of squared residuals provides an unbiased estimator for σ²:
  σ̂² = ||y - Xβ̂||² / (n - p)
  Furthermore, this statistic follows a scaled chi-squared distribution:
  (n - p)σ̂² / σ² ~ χ²(n-p)

6. Towards Hypothesis Testing
The lecture concludes by setting the stage for hypothesis testing. The key insight, derived from Cochran's Theorem, is that the estimator β̂ (which is normally distributed) is independent of the variance estimator σ̂² (which follows a chi-squared distribution). This independence is exactly the condition required to construct a Student's t-statistic to test hypotheses about individual coefficients, such as whether a particular βⱼ is equal to zero.

L14: Regression cont.

This lecture focuses on the core principles of statistical inference in linear regression, building upon the derivation of the least squares estimator. It covers how to perform hypothesis tests on the model's coefficients, address the challenges of multiple testing, and understand the critical assumptions that underpin the entire framework.
1. Recap: The Foundation for Inference
The lecture begins with a summary of the key results for the least squares estimator, β̂, under the Gaussian noise model (y = Xβ + ε, where ε ~ N(0, σ²I)):

• Distribution of the Estimator: The estimator β̂ is itself a random variable with a normal distribution. It is unbiased, meaning its expected value is the true parameter vector β.
  • β̂ ~ N(β, σ²(XᵀX)⁻¹)
• Estimating Noise Variance: The unknown noise variance σ² is estimated using the residuals (ε̂ = y - Xβ̂). The unbiased estimator is:
  • σ̂² = ||y - Xβ̂||² / (n - p)
  • This estimator, when properly scaled, follows a chi-squared distribution: (n - p)σ̂² / σ² ~ χ²(n-p).
• Independence (Cochran's Theorem): Crucially, the estimator β̂ is statistically independent of the noise variance estimator σ̂². This independence is a consequence of their geometric relationship: β̂ relates to the projection of y onto the column space of X, while σ̂² relates to the projection onto its orthogonal complement. This independence is the key that allows for the construction of t-tests.

2. Testing Individual Coefficients: The t-test
A primary goal in regression is to determine which explanatory variables are significantly related to the outcome variable. This is done by testing if their corresponding coefficients are non-zero.

• The Hypothesis: For each coefficient βⱼ, we test:
  1. Null Hypothesis (H₀): βⱼ = 0 (Variable xⱼ has no linear relationship with y).
  2. Alternative Hypothesis (H₁): βⱼ ≠ 0 (Variable xⱼ is a significant predictor).
• Constructing the t-statistic: The process mirrors the construction of a one-sample t-test:
  1. Isolate The distribution of a single estimated coefficient is β̂ⱼ ~ N(βⱼ, σ²γⱼ), where γⱼ is the j-th diagonal element of the (XᵀX)⁻¹ matrix.
  2. Standardize: Under H₀ (βⱼ = 0), the standardized value β̂ⱼ / (σ√γⱼ) follows a standard normal distribution N(0, 1).
  3. Replace Since σ is unknown, it is replaced with its estimate σ̂. This changes the distribution from normal to a Student's t-distribution.
  4. The t-statistic: The final test statistic is:
     T = β̂ⱼ / (σ̂√γⱼ)
     This statistic follows a t-distribution with .
• Decision Rule: We reject H₀ if the absolute value of the calculated T-statistic is greater than the critical value from the t-distribution for a given significance level α (e.g., q_(α/2)). Statistical software automates this by providing a p-value for each coefficient.

3. The Problem of Multiple Testing and the Bonferroni Correction
When a regression model has many variables, analysts typically examine the t-test for every single coefficient. This practice of conducting multiple simultaneous tests creates a problem.

• Inflated Type I Error: If you perform p tests, each at a 5% significance level, the probability of incorrectly rejecting at least one true null hypothesis (a "false discovery") is much higher than 5%.
• The Union Bound: The probability of at least one of several events occurring is less than or equal to the sum of their individual probabilities. This means the overall error rate can be as high as p * α.
• Bonferroni Correction: To control the overall "family-wise" error rate at α, a simple and common (though conservative) solution is the Bonferroni correction. It involves performing each individual test at a much stricter significance level:
  α' = α / p
  By dividing the desired error rate by the number of tests, the union bound guarantees that the overall probability of making any Type I error is at most α.

4. General Linear Hypotheses: The F-test
The t-test is limited to testing a single coefficient at a time. The F-test is a more general tool for testing multiple linear constraints on the coefficients simultaneously.

• The General Hypothesis: The hypothesis is formulated in matrix form:
  1. H₀: Gβ = λ, where G is a matrix defining k linear combinations of the coefficients, and λ is a vector of hypothesized values.
• Examples:
  1. Testing a subset of coefficients: Testing if β₂ = 0 and β₃ = 0 simultaneously.
  2. Testing if the entire model is significant: Testing if all coefficients (except the intercept) are zero. This is equivalent to setting G to an identity matrix and λ to zero.
• Constructing the F-statistic:
  1. The quantity Gβ̂ - λ has a multivariate normal distribution.
  2. Similar to the Wald test, a quadratic form is created using the inverse of this quantity's covariance matrix. This form follows a chi-squared distribution with .
  3. The F-statistic is the ratio of two independent chi-squared variables, each divided by its degrees of freedom:
     F = (Numerator χ²(k) / k) / (Denominator χ²(n-p) / (n-p))
     The numerator comes from the hypothesis being tested, and the denominator is the familiar scaled residual sum of squares (σ̂²).
• The F-distribution: This statistic follows a Fisher (or F) distribution with k and n-p degrees of freedom. The test is one-sided: a large F-value provides evidence against the null hypothesis.

5. Key Assumptions and Practical Considerations
The lecture concludes by highlighting the critical assumptions and potential misinterpretations in regression analysis.

• Correlation is Not Causation: A significant coefficient for a variable xⱼ only indicates a statistical association with y. It does not prove that x causes y. The classic example is finding a positive correlation between police expenditure and crime rates—it's more likely that higher crime rates cause increased police funding, not the other way around.
• The Design Matrix The structure of X is crucial. The variance of the estimators (γⱼ) depends on it. If explanatory variables are highly correlated (not orthogonal), the (XᵀX) matrix is nearly singular, leading to large variances and unstable estimates.
• Normality of the Noise: All the inference procedures discussed (t-tests, F-tests) depend heavily on the assumption that the noise term ε is normally distributed. It is essential to check this assumption, for example, by performing a Kolmogorov-Lilliefors test on the model's residuals.

L15: Bayesian Statistics

The fundamental difference between the two statistical paradigms lies in how they treat the unknown parameter (e.g., θ).

• Frequentist Statistics (What we've done so far):
  1. Data is considered a random sample from a generating process.
  2. The unknown parameter θ is a fixed, non-random constant.
  3. The goal is to create procedures (estimators, confidence intervals) that have good long-run properties over many hypothetical repetitions of the experiment.
• Bayesian Statistics:
  1. Data is also considered random, conditional on a parameter.
  2. The key difference: The parameter θ is treated as a random variable. We don't necessarily believe it's truly random in nature, but modeling it this way provides a mathematical framework for incorporating and updating our beliefs.
  3. The process involves a two-layer model:
     • First, the parameter θ is "drawn" from a prior distribution, π(θ), which represents our initial knowledge or belief about θ before seeing any data.
     • Second, the data X₁,...,Xₙ is drawn from a distribution conditional on this value of θ.
  4. The goal is to use the observed data to update our initial belief into a posterior distribution, π(θ | data), which represents our refined knowledge about θ after accounting for the evidence.

2. The Core Mechanism: Bayes' Rule
The updating of belief from prior to posterior is governed by Bayes' Rule. For a parameter θ and data X₁,...,Xₙ, the rule states:
Posterior ∝ Likelihood × Prior
Mathematically, this is written as:
π(θ | X₁,...,Xₙ) ∝ P(X₁,...,Xₙ | θ) * π(θ)

• (Posterior): The distribution of the parameter θ after observing the data. This is the main output of a Bayesian analysis.
• (Likelihood): The same likelihood function used in frequentist statistics. It is the probability (or density) of observing the data given a specific value of the parameter θ.
• (Prior): The distribution assigned to the parameter θ before observing the data. It quantifies our initial beliefs.
• The proportionality ∝ means we can ignore the normalizing constant P(data), which can be difficult to compute but does not depend on θ.

3. Priors in Practice: Conjugate and Uninformative Priors
Choosing a prior is a central part of the Bayesian workflow.
A. Conjugate Priors
A prior distribution is called conjugate for a given likelihood if the resulting posterior distribution belongs to the same family of distributions as the prior.

• Example: Bernoulli/Beta Model
  • Likelihood: For Bernoulli data (1 for success, 0 for failure), the likelihood for the parameter p is p^(Σxᵢ) * (1-p)^(n-Σxᵢ).
  • Prior: The Beta distribution is a flexible distribution on the interval, controlled by two parameters, a and b. Its density is proportional to p^(a-1) * (1-p)^(b-1).
  • Conjugacy: The Beta distribution is a conjugate prior for the Bernoulli likelihood. If you start with a Beta(a, b) prior, the posterior will also be a Beta distribution.
  • Posterior: The posterior distribution for p is Beta(a + Σxᵢ, b + n - Σxᵢ). The parameters are simply updated by adding the number of successes (Σxᵢ) and failures (n - Σxᵢ) to the initial prior parameters. This makes computation remarkably simple.

B. Uninformative Priors
What if you have no prior knowledge and want the data to speak for itself? You can use an uninformative prior.

• Principle: An uninformative prior attempts to assign equal belief to all possible parameter values.
• Example: Uniform Prior: For a parameter on a bounded interval like the Bernoulli p in, the uninformative prior is the uniform distribution.
• Improper Priors: For parameters on unbounded intervals (like the mean θ of a Gaussian on the real line), a uniform prior would be a constant function that does not integrate to 1. This is called an improper prior. It can still be used as long as the resulting posterior distribution is a proper (integrable) distribution.

4. Example: Bayesian Inference for a Gaussian Mean
Let's assume our data X₁,...,Xₙ comes from a N(θ, 1) distribution, and we want to find the posterior distribution for θ.

• Likelihood: The likelihood is proportional to exp(-½ Σ(xᵢ - θ)²).
• Prior: We use an uninformative (improper) prior for θ, so π(θ) ∝ 1.
• Posterior: The posterior is proportional to the likelihood: π(θ | data) ∝ exp(-½ Σ(xᵢ - θ)²).

To understand this distribution, we complete the square in the exponent with respect to θ. After some algebra, the exponent can be rewritten as:
exp(-n/2 * (θ - X̄ₙ)²), plus terms that do not depend on θ.
This is the kernel of a Gaussian density for θ. We can therefore identify the posterior distribution:
θ | data ~ N(X̄ₙ, 1/n)
This is a powerful result. It states that our updated belief about the mean θ is a normal distribution centered at the sample mean X̄ₙ, and our certainty about this mean increases (the variance 1/n decreases) as we collect more data.
5. Jeffrey's Priors
A more principled, generic type of uninformative prior is Jeffrey's prior, which is defined as being proportional to the square root of the determinant of the Fisher information matrix.

• Intuition: It gives more prior weight to parameter values where the data is expected to provide more information.
• Key Property: It is invariant to reparameterization. If you change your parameterization (e.g., from p to p²), the Jeffrey's prior for the new model is consistent with the old one, which is a highly desirable property.

L16: Bayesian Statistics cont.

This lecture concludes the short chapter on Bayesian inference, delving into principled ways of choosing priors, the interpretation and construction of Bayesian confidence regions, and how Bayesian estimators behave.
1. Principled Priors: The Jeffreys Prior
While priors can encode subjective, domain-specific knowledge, there are also principled, "objective" ways to select them. The lecture focuses on the Jeffreys prior.

• Definition: The Jeffreys prior is defined as being proportional to the square root of the determinant of the Fisher information matrix, I(θ).
  • π(θ) ∝ √det(I(θ))
• Intuition: The Fisher information measures how much information a random variable carries about an unknown parameter. The Jeffreys prior, therefore, gives more weight to parameter values where the data is expected to be more informative.
• Examples:
  • Bernoulli Model: The Fisher information is 1/(p(1-p)). The Jeffreys prior is π(p) ∝ 1/√p(1-p), which is a Beta(1/2, 1/2) distribution. This prior places more weight on values of p near 0 and 1, where observations are more informative.
  • Gaussian Model (Mean): The Fisher information for the mean θ is a constant. Therefore, the Jeffreys prior is π(θ) ∝ 1, which is the same flat, uninformative (and improper) prior discussed previously.
• Key Property: Reparameterization Invariance: The most important theoretical property of the Jeffreys prior is its invariance to reparameterization. If you decide to parameterize your model using a new parameter η = φ(θ) (e.g., θ³ instead of θ), the prior induced on η will be the correct Jeffreys prior for the model parameterized by η. This ensures that the physical meaning of your inference doesn't change just because you changed your mathematical description of the parameter.

2. Bayesian Confidence Regions
The main output of a Bayesian analysis is the entire posterior distribution, π(θ | data). A Bayesian confidence region is a way to summarize this distribution.

• Definition: A Bayesian confidence region R is a subset of the parameter space such that the posterior probability of the true parameter lying within it is 1 - α.
  • P(θ ∈ R | data) = 1 - α
• Construction: A common method is to use level sets of the posterior density. One finds a threshold τ such that the region where π(θ | data) > τ contains exactly 1 - α of the total probability mass. This method naturally produces the smallest possible confidence regions and can handle multi-modal posteriors, resulting in a region that is a union of disjoint intervals.

3. The Crucial Difference: Bayesian vs. Frequentist Interpretation
The interpretation of a confidence region is a major philosophical dividing line between the two paradigms.

• Frequentist Confidence Interval: A 95% confidence interval is a random interval calculated from the data. The interpretation is that if one were to repeat the entire experiment many times, 95% of the calculated intervals would contain the fixed, true parameter . For any single, realized interval, the true parameter is either in it or not—the probability is either 1 or 0. The 95% refers to the long-run reliability of the procedure, not the probability for a specific interval.
• Bayesian Confidence Region: A 95% Bayesian confidence region is a fixed region calculated from the specific data you have. The interpretation is that, given the observed data and the chosen prior, there is a 95% probability that the random parameter lies within this region. It is a direct statement of belief or certainty about the parameter's location, conditional on the evidence at hand.
• Illustrative Example:
  Consider an experiment where two observations are made: x₁ = θ + ε₁ and x₂ = θ + ε₂, with ε being either +1 or -1 with 50% probability each, and θ being an unknown integer.
  • A frequentist confidence region can be constructed that has 75% coverage.
  • However, if you observe x₁ = 5 and x₂ = 7, you know with 100% certainty that θ must be 6 (since ε₁ must have been -1 and ε₂ must have been +1).
  • The frequentist framework cannot update its confidence based on this specific data; the procedure still only has a 75% "guarantee."
  • The Bayesian framework, by conditioning on the observed data (5, 7), calculates a posterior distribution that puts 100% of its probability mass on . This aligns with our logical certainty and highlights the Bayesian ability to make probability statements conditional on the specific evidence.

4. From Posterior to Point Estimate
While the full posterior is the ideal output, a single-point estimate is often needed. Common choices include:

• Posterior Mean: The expected value of θ with respect to the posterior distribution (E[θ | data]). This is the "center of mass" of our updated belief.
• Posterior Mode (Maximum a Posteriori or MAP): The value of θ that maximizes the posterior density. With a flat, uninformative prior, the MAP estimate is identical to the Maximum Likelihood Estimate (MLE).
• Posterior Median: The median of the posterior distribution.

5. The Asymptotic Behavior of Bayesian Methods
As the sample size n grows, the influence of the prior is "washed away" by the overwhelming evidence from the data.

• Example: Bernoulli/Beta Model: The posterior mean for p is (a + Σxᵢ) / (2a + n). As n → ∞, this expression converges to Σxᵢ / n = X̄ₙ, which is the frequentist MLE.
• Example: Gaussian Model: With an uninformative prior, the posterior mean for θ is X̄ₙ for any sample size. The posterior variance is 1/n. The Bayesian posterior distribution N(X̄ₙ, 1/n) mirrors the frequentist sampling distribution of the estimator X̄ₙ.
• General Principle (Bernstein-von Mises Theorem): For large sample sizes, the posterior distribution often approximates a normal distribution centered at the MLE with a variance equal to the inverse of the Fisher information. In essence, with enough data, Bayesian and frequentist results tend to converge. This provides a frequentist justification for Bayesian procedures.

L17: PCA

This lecture introduces Principal Component Analysis (PCA), a fundamental technique for dimension reduction in unsupervised learning. In unsupervised learning, we are given a dataset of multi-dimensional observations without any corresponding outcome or "y" variable. The goal is to understand the inherent structure, patterns, and relationships within the data itself.
1. The Setup: High-Dimensional Data

• Data: We start with n independent observations, X₁, X₂, ..., Xₙ, where each observation Xᵢ is a d-dimensional random vector (a point in ℝᵈ).
• The Challenge: We assume the dimension d is large (e.g., d ≥ 4), making direct visualization and interpretation of the data cloud impossible.
• Matrix Representation: For computational convenience, the data is organized into an n x d matrix, X, where each row represents an observation (Xᵢᵀ) and each column represents a feature or coordinate.

2. The Language of Covariance
PCA is fundamentally based on analyzing the variance and covariance of the data.

• Population Mean and Covariance: For the underlying random vector X, we define:
  • Mean Vector ( μ = E[X], a d-dimensional vector representing the center of the data's distribution.
  • Covariance Matrix ( A d x d symmetric matrix where the entry Σⱼₖ is the covariance between the j-th and k-th coordinates of X. The covariance matrix can be expressed compactly as Σ = E[(X - μ)(X - μ)ᵀ]. It describes the shape and orientation of the data cloud.
• Sample Mean and Covariance: In practice, we don't know the true population parameters. We estimate them from our data:
  • Sample Mean Vector ( x̄ = (1/n) ΣXᵢ, the average of our observed vectors.
  • Sample Covariance Matrix ( S = (1/n) Σ(Xᵢ - x̄)(Xᵢ - x̄)ᵀ. This is the empirical counterpart to Σ.

3. The Matrix Formulation of Sample Covariance
The sample covariance matrix S can be expressed elegantly using the data matrix X.

• The Centering Matrix ( A special n x n matrix H = Iₙ - (1/n)11ᵀ is introduced. This matrix H is an orthogonal projector. When applied to a vector, it subtracts the vector's mean from each of its components.
• Compact Formula: The sample covariance matrix can be written as:
  S = (1/n) XᵀHX
  Geometrically, this means we first center the data by applying the projection H to the columns of X, and then compute the covariance.

4. The Core Idea: Maximizing Variance
The central principle of PCA is to find the directions in the d-dimensional space along which the data has the largest spread or variance. Projecting the data onto these directions preserves the most "information."

• Variance in a Direction A key insight is the link between the covariance matrix and variance. The variance of the data projected onto a direction (a unit vector) u is given by the quadratic form:
  • Population Variance: Var(uᵀX) = uᵀΣu
  • Sample Variance: SampleVar(uᵀXᵢ) = uᵀSu
• The Optimization Problem: PCA seeks to find the direction u that maximizes this projected variance. To prevent an infinite solution (by making u infinitely long), we constrain u to be a unit vector (||u|| = 1). The problem is thus:
  Maximize

5. The Solution: Eigenvalues and Eigenvectors
This seemingly complex, non-convex optimization problem has an elegant and efficient solution from linear algebra, based on the eigen-decomposition (or spectral decomposition) of the sample covariance matrix S.

• Spectral Theorem: Since S is a real, symmetric matrix, it can be diagonalized as:
  S = PDPᵀ
  • D is a diagonal matrix whose entries λ₁, λ₂, ..., λᵈ are the eigenvalues of S. Since S is a covariance matrix, these eigenvalues are all non-negative.
  • P is an orthogonal matrix whose columns v₁, v₂, ..., vᵈ are the corresponding eigenvectors. These eigenvectors form an orthonormal basis for ℝᵈ.
• The Solution: The solution to the variance maximization problem is:
  • The maximum possible variance in any direction is the largest eigenvalue, λ₁.
  • The direction u that achieves this maximum variance is its corresponding eigenvector, v₁.

6. The PCA Algorithm
Principal Component Analysis is the procedure of finding these directions of maximal variance sequentially.

1. First Principal Component (PC1): Find the eigenvector v₁ of the sample covariance matrix S that corresponds to the largest eigenvalue λ₁. This is the direction of maximum variance in the data.
2. Second Principal Component (PC2): Find the eigenvector v₂ corresponding to the second-largest eigenvalue λ₂. This is the direction that has the largest variance among all directions orthogonal to the first principal component.
3. Subsequent Components: Continue this process. The k-th principal component is the eigenvector vₖ corresponding to the k-th largest eigenvalue λₖ.

To perform dimension reduction from d dimensions to a smaller number d', you project the original data onto the subspace spanned by the first d' principal components (v₁, ..., v_d').

L18: PCA cont.

This lecture continues the exploration of Principal Component Analysis (PCA), detailing the algorithm, its geometric interpretation, the process of dimension reduction, and the statistical challenges associated with high-dimensional data.
1. Recap: The Core Idea of PCA
PCA is a technique for transforming a high-dimensional dataset into a new, lower-dimensional coordinate system. The goal is to retain as much of the original data's "information" as possible. In PCA, "information" is defined as variance.

• Mathematical Tool: The foundation of PCA is the spectral (or eigen) decomposition of the sample covariance matrix S.
  • S = PDPᵀ, where P is an orthogonal matrix of eigenvectors (vⱼ) and D is a diagonal matrix of non-negative eigenvalues (λⱼ).
• Interpretation: The eigenvalues λⱼ represent the amount of sample variance along the directions defined by their corresponding eigenvectors vⱼ. By convention, the eigenvalues are ordered from largest to smallest (λ₁ ≥ λ₂ ≥ ... ≥ λᵈ).

2. The PCA Algorithm: A Sequential Approach
PCA finds the optimal low-dimensional subspace by identifying the directions of maximal variance in a step-by-step manner.

1. First Principal Component (PC1): The first principal component is the eigenvector v₁ corresponding to the largest eigenvalue λ₁. This is the single direction in the data space that captures the most variance.
2. Second Principal Component (PC2): The second principal component is the eigenvector v₂ corresponding to the second-largest eigenvalue λ₂. This direction captures the most remaining variance, subject to being orthogonal to the first principal component.
3. Subsequent Components: This process continues. The k-th principal component, vₖ, is the direction that maximizes variance while being orthogonal to all the previously found components (v₁ through vₖ₋₁).

The complete set of eigenvectors v₁, ..., vᵈ forms a new, rotated orthonormal basis for the data space, where the axes are aligned with the directions of decreasing variance.
3. Dimension Reduction in Practice
To reduce data from d dimensions to a smaller dimension k (e.g., k=2 or 3 for visualization):

1. Compute the Covariance Matrix: Calculate the sample covariance matrix S from the data.
2. Find Eigenvectors: Perform an eigen-decomposition of S to find the eigenvectors v₁, ..., vᵈ and their corresponding eigenvalues λ₁, ..., λᵈ.
3. Select Choose the first k eigenvectors (v₁ through vₖ) corresponding to the k largest eigenvalues. These k vectors form the basis of the new, lower-dimensional subspace.
4. Project the Data: Project each original data point Xᵢ onto this subspace. The new k-dimensional coordinates for each point are calculated by Yᵢ = PₖᵀXᵢ, where Pₖ is the d x k matrix whose columns are the first k eigenvectors.

4. How to Choose
Choosing the right number of components to keep is a critical, and often subjective, part of PCA.

• For Visualization: If the goal is visualization, k=2 or k=3 is the obvious choice.
• Proportion of Explained Variance: A more principled approach is to look at the proportion of variance explained. The total variance in the data is the sum of all eigenvalues (the trace of S). The proportion of variance captured by the first k components is:
  (λ₁ + ... + λₖ) / (λ₁ + ... + λᵈ)
  A common heuristic is to choose k such that this ratio is high (e.g., 80% or 90%).
• The Scree Plot: This is a plot of the eigenvalues λⱼ versus their index j. Since the eigenvalues are sorted in decreasing order, the plot will be a downward-sloping curve. A common rule of thumb is to look for an "elbow" or inflection point in the plot. The number of components to keep is the index where the curve flattens out, suggesting that subsequent components contribute only marginal amounts of variance.

5. Interpretation and the Challenge of Sparsity
While PCA provides the mathematically optimal directions of variance, interpreting what these directions mean can be challenging.

• The Interpretation Problem: Each principal component vⱼ is a linear combination of all the original d variables. In high dimensions (e.g., genomics, where d can be hundreds of thousands), a direction that is a dense combination of all genes is not practically interpretable.
• Sparse PCA: To address this, a modern variant called Sparse PCA is introduced. It modifies the optimization problem to find directions that are not only high-variance but also sparse—meaning they have only a few non-zero coefficients. The problem becomes:
  Maximize
  (where ||u||₀ is the number of non-zero elements, and s is a small number like 10).
  While this problem is computationally much harder, its solutions are interpretable. For example, a sparse principal component in a medical dataset might be defined by only 10 key symptoms, making it easy for a doctor to understand.

6. Statistical Theory: The High-Dimensional Setting
The lecture concludes by touching on the statistical theory behind PCA, especially in the challenging modern setting where the dimension d can be comparable to or even larger than the sample size n.

• Classical Asymptotics ( In this regime, the sample covariance S converges to the true covariance Σ, and the sample eigenvectors converge to the true population eigenvectors.
• High-Dimensional Asymptotics ( Results from Random Matrix Theory show that when d is not negligible compared to n, the sample eigenvectors (vₖ of S) do not converge to the true eigenvectors (vₖ of Σ). There is a systematic angle between them that depends on the aspect ratio γ. When d >> n, this angle approaches 90 degrees, meaning the sample components are essentially noise and uncorrelated with the true underlying structure. This highlights the need for techniques like sparse PCA, which impose additional structural assumptions to make learning possible in the d > n regime.

L19: Generalized Linear Models

This lecture introduces Generalized Linear Models (GLMs), a powerful and flexible framework that extends the principles of linear regression to a much broader range of problems.
1. Limitations of the Linear Model
The standard linear model, y = Xβ + ε, is built on several key assumptions:

1. Gaussian Response: The response variable y is assumed to be normally distributed.
2. Linearity: The mean of y is a linear combination of the predictors, E[y|x] = xᵀβ.
3. Constant Variance: The variance of y is constant and does not depend on the predictors.

These assumptions are often violated in practice:

• The response variable might be binary (0/1), a count, or strictly positive.
• The relationship between the predictors and the mean response might be non-linear.

2. The Three Components of a GLM
GLMs generalize the linear model by relaxing these assumptions. Every GLM is defined by three components:

1. Random Component: Specifies the conditional distribution of the response variable y given the predictors x. Instead of being restricted to the Gaussian distribution, this distribution can be any member of the exponential family. This is a broad class of distributions that includes the Gaussian, Bernoulli, Poisson, Binomial, and Gamma distributions, among others.
2. Systematic Component: This is the linear predictor, η = xᵀβ. This part remains the same as in the standard linear model. It specifies that the predictors are combined linearly.
3. Link Function: This component connects the random and systematic parts. The link function g specifies how the mean of the response, μ = E[y|x], is related to the linear predictor η.
   • g(μ) = η = xᵀβ
     The link function's role is to map the range of the mean μ (which might be constrained, e.g., to the interval for a Bernoulli) to the entire real line (the range of xᵀβ).

3. The Exponential Family of Distributions
A distribution belongs to the exponential family if its probability density (or mass) function can be written in the form:
f(y; θ) = exp[ (η(θ)ᵀT(y) - b(θ)) ] * h(y)
In the simpler canonical form, this becomes:
f(y; θ) = exp[ (yθ - b(θ))/φ + c(y, φ) ]

• θ is the canonical parameter.
• φ is the dispersion parameter (often known, like σ² in a Gaussian with known variance).
• b(θ) is the cumulant generating function, which uniquely defines a member of the family.

Why is this family important?

• Breadth: It includes most of the commonly used distributions.
• Properties: Distributions in this family share convenient statistical properties, particularly related to their mean and variance, which can be derived directly from the b(θ) function.

4. Examples of GLMs in Practice

• Example 1: Disease Outbreak (Poisson Regression)
  • Response Number of new cases (a count).
  • Random Component: The Poisson distribution is a natural choice for count data.
  • Systematic Component: The predictor x is time. We hypothesize an exponential growth model, μ = γe^(δx).
  • Link Function: To make this linear, we take the logarithm: log(μ) = log(γ) + δx. This is a linear function of x where β₀ = log(γ) and β₁ = δ. The link function g(μ) = log(μ) is called the log link.
• Example 2: Binary Response (Logistic Regression)
  • Response Presence/absence of a medical condition (0 or 1).
  • Random Component: The Bernoulli distribution is the only choice for a binary outcome. The mean μ is the probability of a "1", so 0 ≤ μ ≤ 1.
  • Systematic Component: A linear combination of predictors xᵀβ, which can take any real value.
  • Link Function: We need a function g that maps the interval (0, 1) to the entire real line. The standard choice is the logit function: g(μ) = log(μ / (1-μ)). This model is known as logistic regression. Other choices for the link function, such as the inverse of the standard normal CDF (probit link), are also possible.

5. Properties of the Canonical Exponential Family
For distributions in the canonical form, the mean and variance of the response variable y have a simple and direct relationship with the cumulant generating function b(θ).

• Mean: The expected value of y is the first derivative of b(θ) with respect to θ.
  • E[y] = μ = b'(θ)
• Variance: The variance of y is the second derivative of b(θ) times the dispersion parameter φ.
  • Var(y) = b''(θ) * φ

This relationship is proven using the identity E[∂log(L)/∂θ] = 0, a property used in the derivation of the Cramér-Rao lower bound. A key consequence is that b''(θ) = Var(y)/φ, which must be positive. This implies that the function b(θ) is always convex, a property that is crucial for the optimization of the likelihood function in GLMs.

L20: Generalized Linear Models cont.

This lecture dives deeper into the mechanics of Generalized Linear Models (GLMs), focusing on the properties derived from the exponential family and the role of the link function in connecting the model components.
1. Properties of the Canonical Exponential Family
The lecture begins by revisiting the canonical exponential family, a simpler, one-parameter subclass of the exponential family. Its probability density/mass function has the form:
f(y; θ) = exp[ (yθ - b(θ))/φ + c(y, φ) ]

• θ is the canonical parameter.
• φ is a known dispersion parameter.
• b(θ) is the cumulant generating function, a key function that encapsulates the properties of the distribution.

A crucial insight is that the mean and variance of the response variable y can be derived directly from b(θ) using calculus, without needing to perform explicit integration or summation.

• Derivation: This is achieved using two fundamental identities related to the log-likelihood function l(θ; y) = log(f(y; θ)).
  • E[∂l/∂θ] = 0 (The expected value of the score function is zero).
  • E[∂²l/∂θ²] = -E[(∂l/∂θ)²] (The expected value of the Hessian is the negative of the Fisher information).
• Resulting Properties: By applying these identities to the log-likelihood of the canonical exponential family, we find:
  • Mean: E[y] = b'(θ) (The mean is the first derivative of b(θ)).
  • Variance: Var(y) = b''(θ)φ (The variance is the second derivative of b(θ) times the dispersion parameter).
• Convexity of Since Var(y) and φ are always positive, it follows that b''(θ) must be positive. This means that the cumulant generating function b(θ) is always a strictly convex function. This convexity is a cornerstone property that ensures the likelihood function for a GLM is well-behaved and can be efficiently optimized.

2. The Role of the Link Function
The link function g is the bridge between the mean of the response, μ = E[y|x], and the linear predictor, η = xᵀβ.

• The Link Equation: g(μ) = xᵀβ
• Purpose: The link function's primary job is to resolve the mismatch between the range of μ and the range of xᵀβ.
  • xᵀβ can take any value on the real line (-∞, +∞).
  • μ is often constrained. For a Poisson count, μ > 0. For a Bernoulli probability, μ ∈ (0, 1).
• The link function g must be a strictly increasing, invertible function that maps the domain of μ to the entire real line.

3. Common Link Functions

• Identity Link ( Used for the standard linear model where y is Gaussian. μ can take any real value, so no transformation is needed.
• Log Link ( Used for models where the mean must be positive, such as Poisson regression. It maps (0, ∞) to (-∞, ∞).
• Logit Link ( The standard link for binary data (logistic regression). It maps (0, 1) to (-∞, ∞).
• Probit Link ( An alternative to the logit link for binary data, where Φ⁻¹ is the inverse of the standard normal CDF.
• Reciprocal Link ( Used for models like the Gamma distribution.

4. The Canonical Link Function
For each member of the canonical exponential family, there is a "natural" or canonical link function.

• Definition: The canonical link is the function g that sets the linear predictor xᵀβ equal to the canonical parameter θ.
  • g(μ) = θ
• Derivation: Since we know μ = b'(θ), it follows that θ = (b')⁻¹(μ). Therefore, the canonical link function is the inverse of the first derivative of the cumulant generating function:
  • g(μ) = (b')⁻¹(μ)
• Why is it "Canonical"? It provides the most direct and mathematically simplest connection between the systematic component (xᵀβ) and the parameter (θ) that governs the probability distribution. For example:
  • Gaussian: b(θ) = θ²/2, b'(θ) = θ. The inverse is the identity function. The canonical link is the identity link.
  • Poisson: b(θ) = e^θ, b'(θ) = e^θ. The inverse is the natural logarithm. The canonical link is the log link.
  • Bernoulli: b(θ) = log(1+e^θ), b'(θ) = e^θ/(1+e^θ). The inverse is log(μ/(1-μ)). The canonical link is the logit link.

5. Estimation in GLMs
With the model fully specified, the goal is to estimate the coefficient vector β from the data.

• The Log-Likelihood Function: We have n independent (but not identically distributed) observations (xᵢ, yᵢ). The total log-likelihood is the sum of the individual log-likelihoods:
  l(β) = Σ [ (yᵢθᵢ - b(θᵢ))/φ + c(yᵢ, φ) ]
• Connecting to The canonical parameter θᵢ for each observation depends on β through the link function:
  θᵢ = (b')⁻¹(μᵢ) = (b')⁻¹(g⁻¹(xᵢᵀβ))
• Canonical Link Simplification: If we use the canonical link, g = (b')⁻¹, then these functions cancel, and the relationship simplifies dramatically:
  θᵢ = xᵢᵀβ
  The log-likelihood becomes: l(β) = Σ [ (yᵢ(xᵢᵀβ) - b(xᵢᵀβ))/φ + c(yᵢ, φ) ]
• Optimization: This log-likelihood function l(β) is a concave function of β (due to the convexity of b(θ)). This is a crucial property, as it guarantees that a unique maximum exists and can be found efficiently using numerical optimization algorithms. The lecture concludes by setting up the next topic: the Newton-Raphson method and its specialization for GLMs, known as Iteratively Reweighted Least Squares (IRLS), which is the standard algorithm for fitting these models.

L21: Generalized Linear Models cont.

This lecture delves into the core mechanics of fitting Generalized Linear Models (GLMs), focusing on maximum likelihood estimation. It bridges the theoretical properties of the exponential family with the practical algorithms used to find the model parameters.
1. Setting Up the Likelihood Function
The goal is to estimate the parameter vector β using the principle of maximum likelihood. This requires constructing the log-likelihood function for all n observations.

• Model Recap: For each observation i:
  1. The response Yᵢ given predictors xᵢ follows a distribution from the canonical exponential family with canonical parameter θᵢ.
  2. The mean μᵢ = E[Yᵢ|xᵢ] is related to θᵢ by μᵢ = b'(θᵢ).
  3. The mean is linked to the linear predictor by a link function g: g(μᵢ) = xᵢᵀβ.
• From To write the likelihood in terms of β, we need to express θᵢ as a function of β. This is done by composing the chain of relationships:
  θᵢ = (b')⁻¹(μᵢ) = (b')⁻¹(g⁻¹(xᵢᵀβ))
  This composite function (b')⁻¹ ∘ g⁻¹ is denoted as h, so θᵢ = h(xᵢᵀβ).
• The Log-Likelihood: The total log-likelihood for n independent observations is the sum of the individual log-likelihoods:
  l(β) = Σ [ (Yᵢθᵢ - b(θᵢ))/φ + c(Yᵢ, φ) ]
  Substituting θᵢ = h(xᵢᵀβ), we get the function l(β) that we need to maximize.
• Canonical Link Simplification: The process simplifies significantly if we use the canonical link function, where g = (b')⁻¹. In this case, h becomes the identity function, and the relationship is simply θᵢ = xᵢᵀβ. The log-likelihood becomes:
  l(β) = Σ [ (Yᵢ(xᵢᵀβ) - b(xᵢᵀβ))/φ ] (ignoring constant terms).

2. The Convexity of the Log-Likelihood
A key property of GLMs with a canonical link is that the log-likelihood function l(β) is concave.

• Intuition: The term Yᵢ(xᵢᵀβ) is linear in β. The term -b(xᵢᵀβ) is the negative of a convex function composed with a linear function, which remains convex. Therefore, -b is concave. The sum of a linear function and a concave function is concave.
• Significance: Because l(β) is concave, it has a unique maximum. This guarantees that we can find a single best estimate for β and that optimization algorithms will converge to this global maximum.

3. Numerical Optimization: The Newton-Raphson Method
Except for the simple Gaussian case (which is just linear regression), the equation ∇l(β) = 0 does not have a closed-form solution. We must use an iterative numerical method to find the maximum. The standard algorithm is the Newton-Raphson method.

• The Idea: Newton-Raphson works by iteratively creating a local quadratic approximation of the function and finding the maximum of that approximation.
• Taylor Expansion: At a current estimate β⁽ᵏ⁾, the log-likelihood is approximated by its second-order Taylor expansion:
  l(β) ≈ l(β⁽ᵏ⁾) + ∇l(β⁽ᵏ⁾)ᵀ(β - β⁽ᵏ⁾) + ½(β - β⁽ᵏ⁾)ᵀH(β⁽ᵏ⁾)(β - β⁽ᵏ⁾)
  where ∇l is the gradient and H is the Hessian matrix (matrix of second derivatives).
• The Update Step: Finding the β that maximizes this quadratic approximation leads to the Newton-Raphson update rule:
  β⁽ᵏ⁺¹⁾ = β⁽ᵏ⁾ - [H(β⁽ᵏ⁾)]⁻¹∇l(β⁽ᵏ⁾)
  This process is repeated until the estimates converge.

4. Fisher Scoring: A Statistical Simplification
The Newton-Raphson method requires computing and inverting the Hessian matrix H at each step. The Fisher Scoring algorithm simplifies this by replacing the observed Hessian with its expected value.

• The Identity: We know from likelihood theory that E[H] = -I(β), where I(β) is the Fisher Information matrix.
• Fisher Scoring Update: The update rule becomes:
  β⁽ᵏ⁺¹⁾ = β⁽ᵏ⁾ + [I(β⁽ᵏ⁾)]⁻¹∇l(β⁽ᵏ⁾)
  This is often computationally simpler because the Fisher Information can sometimes have a cleaner form than the Hessian, especially when the Hessian depends on the random variable Yᵢ. (Note: For GLMs with a canonical link, the Hessian does not depend on Yᵢ, and thus Newton-Raphson and Fisher Scoring are identical).

5. Iteratively Reweighted Least Squares (IRLS)
The final step is to show that the Fisher Scoring update for a GLM can be cleverly rewritten as a weighted least squares problem. This algorithm is known as Iteratively Reweighted Least Squares (IRLS) and is the standard method for fitting GLMs.

• Derivation: Through careful application of the chain rule and using the properties μᵢ = b'(θᵢ) and Var(Yᵢ) = b''(θᵢ)φ, the Fisher Scoring update step can be shown to be equivalent to solving the following for β⁽ᵏ⁺¹⁾:
  β⁽ᵏ⁺¹⁾ = (XᵀWX)⁻¹XᵀWZ
• Components of IRLS:
  1. X is the usual design matrix.
  2. W is a diagonal weights matrix where the weights wᵢᵢ change at each iteration. These weights depend on the variance of Yᵢ and the derivative of the link function.
  3. Z is a "working response" vector, also updated at each iteration. It is constructed from the linear predictor and the residuals.
• The Algorithm:
  1. Initialize β⁽⁰⁾.
  2. Iterate:
     a. Using β⁽ᵏ⁾, calculate the current fitted values μᵢ, linear predictors ηᵢ, weights wᵢᵢ, and the working response zᵢ.
     b. Solve the weighted least squares problem to get the next estimate β⁽ᵏ⁺¹⁾.
  3. Repeat until β converges.

This algorithm is elegant because it recasts a complex non-linear optimization problem into a sequence of familiar linear regression problems, which can be solved efficiently.

L22: Generalized Linear Models cont.

This lecture details the numerical algorithms used to perform maximum likelihood estimation for Generalized Linear Models (GLMs). Since a closed-form solution for the parameter vector β generally does not exist, iterative methods are required.
1. Recap: The GLM Log-Likelihood
The foundation for estimation is the log-likelihood function. For n independent observations (xᵢ, yᵢ), where each yᵢ follows a distribution from the canonical exponential family, the log-likelihood is:
l(β) = Σ [ (yᵢθᵢ - b(θᵢ))/φ + c(yᵢ, φ) ]
The canonical parameter θᵢ is linked to the linear predictor xᵢᵀβ via the composite function h = (b')⁻¹ ∘ g⁻¹, such that θᵢ = h(xᵢᵀβ). When using the canonical link g = (b')⁻¹, this simplifies to θᵢ = xᵢᵀβ.
Crucially, this log-likelihood function is concave, which guarantees that a unique maximum exists and can be found by numerical optimization algorithms.
2. The Newton-Raphson Method
Newton-Raphson is a general-purpose, second-order optimization algorithm for finding the maximum (or minimum) of a function.

• Core Idea: It works by iteratively approximating the function with a quadratic (its second-order Taylor expansion) around the current estimate and then jumping to the maximum of that quadratic approximation.
• The Update Rule: Given a current estimate β⁽ᵏ⁾, the next estimate β⁽ᵏ⁺¹⁾ is calculated as:
  β⁽ᵏ⁺¹⁾ = β⁽ᵏ⁾ - [H(β⁽ᵏ⁾)]⁻¹∇l(β⁽ᵏ⁾)
  • ∇l(β⁽ᵏ⁾) is the gradient (vector of first partial derivatives) of the log-likelihood evaluated at β⁽ᵏ⁾.
  • H(β⁽ᵏ⁾) is the Hessian matrix (matrix of second partial derivatives) evaluated at β⁽ᵏ⁾.
• Properties: Newton-Raphson typically converges very quickly (quadratically) when close to the solution but requires computing and inverting the p x p Hessian matrix at each step, which can be computationally expensive.

3. Fisher Scoring: A Statistical Adaptation
The Fisher Scoring algorithm is a modification of Newton-Raphson specifically for statistical estimation problems. It simplifies the update step by replacing the observed Hessian matrix with its expected value.

• The Identity: From likelihood theory, we know that the expected value of the Hessian is the negative of the Fisher Information matrix, I(β).
  • E[H(β)] = -I(β)
• The Update Rule: The Fisher Scoring update rule is:
  β⁽ᵏ⁺¹⁾ = β⁽ᵏ⁾ + [I(β⁽ᵏ⁾)]⁻¹∇l(β⁽ᵏ⁾)
• Advantages: The Fisher Information matrix can sometimes have a simpler mathematical form than the Hessian, especially if the Hessian depends on the random response Yᵢ. Taking the expectation can average out this randomness, leading to a cleaner and sometimes more stable algorithm. For GLMs with a canonical link, the Hessian does not depend on Yᵢ, and in this important case, Fisher Scoring and Newton-Raphson are identical.

4. Iteratively Reweighted Least Squares (IRLS)
The most common algorithm for fitting GLMs is Iteratively Reweighted Least Squares (IRLS). It is a brilliant statistical trick that shows the Fisher Scoring update step is equivalent to solving a weighted least squares problem at each iteration.

• Weighted Least Squares (WLS): WLS is a modification of ordinary least squares (OLS) used when observations have unequal variance (heteroscedasticity). It minimizes the weighted sum of squared residuals, Σwᵢ(yᵢ - xᵢᵀβ)². The solution is:
  β̂_WLS = (XᵀWX)⁻¹XᵀWy
  where W is a diagonal matrix of weights.
• The IRLS Connection: Through careful derivation involving the chain rule, it can be shown that the Fisher Scoring update β⁽ᵏ⁺¹⁾ = β⁽ᵏ⁾ + [I(β⁽ᵏ⁾)]⁻¹∇l(β⁽ᵏ⁾) is mathematically identical to the solution of a specific WLS problem.
  The solution β⁽ᵏ⁺¹⁾ is found by performing a weighted least squares regression of a "working response" Z on the predictors X with a specific set of weights W.
• The IRLS Algorithm:
  1. Initialize: Obtain starting values for β⁽⁰⁾ and the corresponding means μ⁽⁰⁾.
  2. Iterate (at step k):
     a. Calculate the weights matrix W. W is a diagonal matrix with weights wᵢᵢ that depend on the current estimate β⁽ᵏ⁾ (via μᵢ and the variance and link functions).
     b. Calculate the working response vector Z. Each element zᵢ also depends on the current estimate β⁽ᵏ⁾ (via the linear predictor ηᵢ, the observed yᵢ, and the mean μᵢ).
     c. Solve the WLS problem: Update the estimate by computing:
     β⁽ᵏ⁺¹⁾ = (XᵀWX)⁻¹XᵀWZ
  3. Repeat step 2 until the estimates for β converge.

This algorithm is powerful because it recasts a potentially complex non-linear optimization into a sequence of well-understood and efficiently solvable linear algebra problems (weighted least squares).