Math for Data Science · MDS-6
Turn data into defensible conclusions — estimate, test, and quantify uncertainty with confidence intervals, the bootstrap, A/B tests, regression, and Bayes, all in code.
Statistical inference is how data scientists move from a sample to a claim about the world while honestly reporting uncertainty. This course builds the toolkit end to end — descriptive summaries and sampling, estimation and maximum likelihood, confidence intervals and the bootstrap, hypothesis testing, experiment design and A/B testing, regression and GLMs, and a first look at Bayesian inference. Every idea is paired with a hands-on analysis in Python (SciPy/statsmodels) or R.
Course Outline
Every lesson is self-contained: a full explanation, worked step-by-step examples (in Python and R/Julia), common mistakes, a try-it-yourself, and an interactive quiz. Jump to a unit:
Inference starts by separating what we want to know (the population) from what we actually observe (the sample).
A population is the complete set of units you care about; a parameter is a fixed, usually unknown number describing it (e.g., the true mean income mu). Because we rarely measure everyone, we take a sample and compute a statistic (e.g., the sample mean x-bar) to estimate the parameter. The whole of statistical inference is using a known statistic to make a defensible, uncertainty-aware claim about an unknown parameter. Notation matters: Greek letters (mu, sigma, p) denote parameters, Latin letters (x-bar, s, p-hat) denote statistics. A statistic is random because it changes from sample to sample; a parameter is not. Confusing the two is the root cause of most beginner errors, because every confidence interval and hypothesis test is fundamentally a statement that connects a random statistic to a fixed parameter.
Worked Example 1
Problem. A city has 50,000 households. You survey 200 of them and find an average of 1.8 cars per household. Identify the population, sample, parameter, and statistic.
Answer. Population = 50,000 households; sample = 200 surveyed; parameter = mu (unknown true mean); statistic = x-bar = 1.8 cars.
Worked Example 2
Problem. From a sample, p-hat = 0.42 favor a policy. A report says 'the population proportion is 0.42.' What is wrong, and how should it be phrased?
Answer. The report conflates a statistic with the parameter; p-hat = 0.42 estimates the unknown p and should be reported with uncertainty.
Worked Example 3
Problem. Two analysts each draw a different sample of size 30 from the same population and get x-bar = 51.2 and x-bar = 48.7. Does this mean one is wrong?
Answer. Both are valid estimates of the same fixed mu; differing values reflect normal sampling variability, not error.
Problem. A factory produces millions of bolts. An inspector measures 150 bolts and finds a mean length of 20.04 mm with the goal of judging the true mean length. Label the population, sample, parameter, and statistic, and state which value is known.
Solution. Population = all bolts produced. Sample = the 150 measured bolts. Parameter = mu, the true mean length of all bolts (unknown). Statistic = x-bar = 20.04 mm (known, computed from the sample). The known value is the statistic x-bar = 20.04 mm; it is our estimate of the unknown parameter mu.
Summaries compress a dataset into a few numbers describing where it sits, how spread out it is, and whether it is symmetric.
Center locates a typical value: the mean (sum divided by n) is the balance point but is pulled by outliers; the median (middle rank) and mode (most frequent) resist them. Spread measures variability: variance is the mean squared deviation from the mean, and standard deviation (its square root) is in original units. Shape describes asymmetry via skewness: right-skewed data have a long upper tail and mean > median; left-skewed have mean < median. In data science these summaries drive feature understanding, outlier detection, and the choice of model. A key subtlety is the divisor in variance: dividing by n - 1 (Bessel's correction) gives an unbiased estimate of the population variance from a sample, whereas dividing by n underestimates it. Always match center to spread: report mean with standard deviation, or median with IQR.
Worked Example 1
Problem. For the data 2, 4, 4, 6, 9, compute the mean, median, and mode.
Answer. Mean = 5, median = 4, mode = 4.
Worked Example 2
Problem. Compute the sample standard deviation of 2, 4, 4, 6, 9.
Answer. Sample variance = 7, sample standard deviation = 2.646 (using n - 1).
Worked Example 3
Problem. Incomes (in thousands) are 30, 32, 35, 38, 200. Compare mean vs. median and describe the skew.
Answer. Mean = 67, median = 35; mean > median signals strong right skew driven by the outlier, so the median better represents a typical income.
Problem. For the values 10, 12, 12, 13, 18, compute the mean, median, sample standard deviation, and state the likely skew direction.
Solution. Mean = (10+12+12+13+18)/5 = 65/5 = 13. Median = middle of sorted = 12. Deviations from 13: -3,-1,-1,0,5; squares 9,1,1,0,25 sum to 36; s^2 = 36/(5-1) = 9; s = 3. Mean (13) > median (12), so mild right skew. In Python: np.mean, np.median, np.std(d, ddof=1) confirm 13, 12, 3.
Quantiles slice a distribution by rank, giving robust summaries that ignore extreme values.
A quantile (or percentile) is the value below which a given fraction of data falls: the 0.25 quantile (Q1, 25th percentile) has a quarter of the data below it, and the median is the 0.5 quantile. The five-number summary (min, Q1, median, Q3, max) and the interquartile range IQR = Q3 - Q1 describe a distribution without assuming a shape. These are robust: a single huge outlier barely moves a quantile, unlike the mean. In data science, quantiles power percentile-based SLAs (e.g., p95 latency), winsorizing, and outlier rules. A common outlier fence flags points beyond Q1 - 1.5*IQR or Q3 + 1.5*IQR. Because quantile estimation depends on an interpolation convention, different tools can report slightly different values, so know which method your library uses.
Worked Example 1
Problem. For 1, 3, 5, 7, 9, 11, 13, find Q1, the median (Q2), Q3, and the IQR.
Answer. Q1 = 3, median = 7, Q3 = 11, IQR = 8 (by the simple split method; library defaults may interpolate to Q1=4, Q3=10).
Worked Example 2
Problem. Latency data has p95 = 420 ms. Explain what that means for an SLA.
Answer. p95 = 420 ms means 95% of requests are at or under 420 ms; the SLA target of 500 ms is satisfied.
Worked Example 3
Problem. Using the 1.5*IQR rule, identify outliers in 4, 5, 6, 7, 8, 9, 40.
Answer. 40 is the only outlier, since it exceeds the upper fence of 15.
Problem. For 2, 4, 6, 8, 10, 12, 100, compute Q1, Q3, the IQR, and use the 1.5*IQR rule to check whether 100 is an outlier.
Solution. Median = 8 (4th of 7). Lower half 2,4,6 -> Q1 = 4; upper half 10,12,100 -> Q3 = 12. IQR = 12 - 4 = 8. Upper fence = 12 + 1.5*8 = 24. Since 100 > 24, 100 is an outlier. (np.percentile may interpolate to Q1=5, Q3=11, IQR=6, fence=20; 100 is still an outlier.)
Three core plots reveal a distribution's shape, spread, outliers, and exact percentiles.
A histogram bins values and shows counts as bars, revealing modality and skew; the bin width is a tuning choice that can hide or invent structure. A boxplot draws the five-number summary as a box (Q1 to Q3) with a median line and whiskers, marking 1.5*IQR outliers as points; it compresses spread and outliers into one compact glyph, excellent for comparing groups side by side. An empirical cumulative distribution function (ECDF) plots, for each x, the fraction of data <= x as a rising step function from 0 to 1; unlike a histogram it uses no binning, so it shows every data point and lets you read any percentile directly off the curve. Together these let a data scientist judge shape, detect outliers, and compare distributions before modeling.
Worked Example 1
Problem. Sketch what the ECDF value is at x = 5 for the data 1, 3, 5, 7, 9.
Answer. ECDF(5) = 0.6 — 60% of the data is at or below 5.
Worked Example 2
Problem. A boxplot has Q1 = 20, median = 25, Q3 = 40 with a point plotted at 75. Is 75 an outlier and why is it shown separately?
Answer. Yes, 75 exceeds the upper fence of 70, so the boxplot displays it as an individual outlier point.
Worked Example 3
Problem. A histogram of response times has a tall bar near 0 and a long thin tail to the right. Name the shape and which plot best reads the median.
Answer. The distribution is right-skewed; the ECDF best pinpoints the median by reading the x where the curve reaches 0.5.
Problem. For the data 2, 2, 4, 6, 6, 6, 9, give the ECDF value at x = 6 and describe what a histogram with bin width 2 starting at 1 would show as the tallest bar.
Solution. ECDF(6): values <= 6 are 2,2,4,6,6,6 = 6 of 7, so ECDF(6) = 6/7 ~ 0.857. Bins [1,3),[3,5),[5,7),[7,9),[9,11): counts are 2 (the two 2s), 1 (the 4), 3 (the three 6s), 0, 1 (the 9). The tallest bar is the [5,7) bin with 3 observations, reflecting the mode at 6.
How you select a sample determines whether your estimates are trustworthy or systematically wrong.
A sampling design is the rule for choosing units. Simple random sampling gives every unit an equal chance and is the gold standard; stratified sampling draws within homogeneous subgroups to improve precision; cluster sampling samples whole groups for convenience. Bias is a systematic error that pushes estimates away from the truth regardless of sample size. Selection bias arises when the sampling method over- or under-represents parts of the population (e.g., convenience or volunteer samples); nonresponse bias arises when non-responders differ from responders; survivorship bias arises when only 'surviving' units are observable. Crucially, bias is not cured by a bigger sample — a large biased sample is confidently wrong. For data scientists, recognizing how data was collected is often more important than the analysis applied to it.
Worked Example 1
Problem. An app surveys only users who opened it this week to estimate satisfaction of all registered users. Name the bias.
Answer. Selection bias: sampling only active users excludes dissatisfied or churned users, inflating the satisfaction estimate.
Worked Example 2
Problem. A population is 70% region A and 30% region B. To estimate mean spending precisely, design a stratified sample of size 100.
Answer. Draw 70 random units from region A and 30 from region B (proportional stratified sampling), reducing variance versus a simple random sample.
Worked Example 3
Problem. A WWII analysis examined returning planes' bullet holes to decide where to add armor and proposed reinforcing the most-hit areas. What bias is at play?
Answer. Survivorship bias: the sample omits downed planes, so armor belongs on the unhit areas of survivors, not the hit areas.
Problem. A radio station asks listeners to call in to vote on a city policy, and 5,000 people respond. Why might this estimate of citywide support be biased, and what design would be better?
Solution. Call-in voting is a self-selected (volunteer) sample, so it suffers selection and nonresponse bias: only listeners with strong opinions and access to the station respond, who may differ systematically from the general population. The large n = 5,000 does not fix this. A better design is a simple random or stratified random sample of all residents (e.g., random-digit dialing across the city with follow-ups to reduce nonresponse), giving every resident a known, roughly equal chance of inclusion.
A statistic computed from random samples has its own distribution; its spread is the standard error.
If you could draw every possible sample of size n and compute a statistic (say x-bar) for each, the distribution of those values is the sampling distribution. Its standard deviation is the standard error (SE), which measures how much the statistic bounces around the parameter from sample to sample. For the sample mean, SE = sigma / sqrt(n), so precision improves with the square root of sample size: quadrupling n halves the SE. The Central Limit Theorem says that for large n the sampling distribution of the mean is approximately normal regardless of the population's shape, which justifies normal-based confidence intervals and tests. When sigma is unknown we estimate SE with s / sqrt(n). The SE is the bridge from a single sample to a statement of uncertainty, making it the backbone of all the inference that follows.
Worked Example 1
Problem. A population has sigma = 12. For a sample of n = 36, find the standard error of the mean.
Answer. The standard error of the mean is 2.
Worked Example 2
Problem. With sigma = 12, how large must n be to cut the standard error to 1?
Answer. n = 144 yields SE = 1; halving the SE requires quadrupling the sample size.
Worked Example 3
Problem. A sample of n = 50 has s = 8 (sigma unknown). Estimate the SE of the mean and explain why CLT matters here.
Answer. Estimated SE ~ 1.131; the CLT lets us treat x-bar as approximately normal at n = 50, enabling normal/t-based inference.
Problem. A population has sigma = 20. (a) Find the SE of the mean for n = 100. (b) What n gives an SE of 1? (c) If only s = 18 is known with n = 81, estimate the SE.
Solution. (a) SE = 20/sqrt(100) = 20/10 = 2. (b) Set 20/sqrt(n) = 1 -> sqrt(n) = 20 -> n = 400. (c) With sigma unknown, SE ~ s/sqrt(n) = 18/sqrt(81) = 18/9 = 2. In Python: 20/np.sqrt(100)=2.0, (20)**2=400, 18/np.sqrt(81)=2.0.
Load a public dataset (e.g., a CSV of city temperatures or housing prices) in Python or R and compute mean, median, standard deviation, IQR, and skewness. Produce a histogram, boxplot, and ECDF and describe the distribution's shape and any outliers. Then draw 1000 random samples of size n=30 from the data, record each sample mean, and plot the resulting sampling distribution alongside its standard error.
Deliverable · A notebook or script with the summary table, the three distribution plots, and the sampling-distribution plot, plus a paragraph interpreting the shape and the observed standard error.
1. Which statement correctly distinguishes a parameter from a statistic?
Answer B. A parameter is a fixed (often unknown) population quantity; a statistic is calculated from observed sample data and estimates the parameter.
2. A dataset of incomes has a few very large values. Which measure of center is most resistant to those outliers?
Answer B. The median depends only on rank order, so extreme values barely move it, whereas the mean is pulled toward large outliers.
3. As the sample size n increases, the standard error of the sample mean:
Answer C. The standard error of the mean equals sigma over the square root of n, so it falls as n grows, just more slowly than 1 over n.
4. A convenience sample of volunteers is used to estimate average exercise time in a city. The main concern is:
Answer B. Convenience samples are prone to selection bias: people who volunteer may exercise more (or less) than the general population, biasing the estimate.
5. The interquartile range (IQR) measures:
Answer C. The IQR is Q3 minus Q1, capturing the range of the central half of the data and ignoring the extreme tails.
I can summarize a dataset with appropriate measures of center, spread, and shape and choose robust alternatives when needed.
I can evaluate a sampling design for bias and explain how the sampling distribution of a statistic gives rise to standard error.
I can produce and read histograms, boxplots, and ECDFs to describe a distribution.
An estimator is judged by whether it is centered on the truth, how much it varies, and whether it improves with data.
A point estimator theta-hat is a function of the sample used to guess a parameter theta. Three properties characterize it. Bias = E[theta-hat] - theta measures systematic error: an unbiased estimator is centered on the truth on average. Variance = Var(theta-hat) measures how much it scatters across samples. Consistency means theta-hat converges in probability to theta as n grows, so with enough data it lands on the truth. These can trade off: a slightly biased estimator can have much smaller variance and thus be preferable. For data scientists this framework explains why we choose s^2 with n - 1 (unbiased for variance) and why the sample mean is both unbiased and consistent for mu. Crucially, bias and variance are properties of the estimator's sampling distribution, not of any one estimate.
Worked Example 1
Problem. Show the sample mean x-bar is an unbiased estimator of the population mean mu.
Answer. E[x-bar] = mu, so the sample mean is unbiased for the population mean.
Worked Example 2
Problem. Why does dividing the sum of squared deviations by n (instead of n - 1) give a biased variance estimator?
Answer. Dividing by n underestimates sigma^2 by a factor (n-1)/n; using n - 1 removes the bias.
Worked Example 3
Problem. Estimator A has bias 0 and variance 10; estimator B has bias 1 and variance 2. Which has smaller MSE, and is the sample mean consistent for mu?
Answer. Estimator B (MSE = 3) beats A (MSE = 10); and yes, x-bar is consistent because its variance vanishes as n grows.
Problem. Estimator P of theta has bias = 2 and variance = 4; estimator Q has bias = 0 and variance = 9. Compute each MSE and say which you would prefer, and state whether either being biased rules it out.
Solution. MSE = variance + bias^2. MSE_P = 4 + 2^2 = 8. MSE_Q = 9 + 0^2 = 9. P has the smaller MSE (8 < 9), so prefer P even though it is biased; bias alone does not disqualify an estimator, because lower variance can more than compensate.
MSE bundles bias and variance into one number, exposing the central tradeoff of estimation and modeling.
Mean squared error is MSE(theta-hat) = E[(theta-hat - theta)^2], the expected squared distance from the truth. It decomposes exactly as MSE = Var(theta-hat) + Bias(theta-hat)^2. This identity is the bias-variance tradeoff: reducing one often raises the other, and the best estimator (or model) minimizes their sum, not either alone. In machine learning the same decomposition governs over- and under-fitting: a flexible model has low bias but high variance (overfits), while a rigid model has high bias but low variance (underfits). Regularization deliberately adds bias to cut variance and lower total error. Understanding MSE lets a data scientist reason about why a more complex model can predict worse on new data even though it fits the training data better.
Worked Example 1
Problem. An estimator has bias = 3 and variance = 16. Compute its MSE.
Answer. MSE = 25.
Worked Example 2
Problem. Consider estimating mu with the shrinkage estimator c * x-bar for constant c in [0,1], where Var(x-bar) = sigma^2/n. Express its bias and variance.
Answer. Bias = (c - 1)*mu, Variance = c^2 * sigma^2/n, MSE = c^2 sigma^2/n + (c-1)^2 mu^2; c < 1 trades added bias for reduced variance.
Worked Example 3
Problem. A flexible model has training MSE 0.5 but test MSE 8; a simpler model has training MSE 3 and test MSE 4. Interpret via bias-variance.
Answer. The flexible model overfits (low bias, high variance); the simpler model wins on test MSE (4 vs 8), illustrating the bias-variance tradeoff.
Problem. Estimator R has variance 5 and bias 4; estimator S has variance 18 and bias 1. Compute both MSEs and explain which generalizes better.
Solution. MSE_R = 5 + 4^2 = 5 + 16 = 21. MSE_S = 18 + 1^2 = 18 + 1 = 19. S has the lower MSE (19 < 21), so it has smaller expected squared error and generalizes better, even though R has lower variance — its large bias dominates.
Match sample moments to their theoretical formulas and solve for the parameters — a quick, intuitive estimator.
The method of moments (MoM) estimates parameters by equating sample moments to their model-implied counterparts and solving. The k-th sample moment is (1/n) sum X_i^k; the first is the sample mean, the second relates to variance. If a distribution has p parameters, you write p moment equations and solve the system. MoM is computationally cheap and needs no optimization, making it a fast starting point or initializer for MLE. It is consistent but generally not as efficient (higher variance) as the MLE, and it can produce out-of-range estimates. In data science MoM is handy for quick distribution fits, sanity checks, and providing good starting values for numerical likelihood maximization.
Worked Example 1
Problem. Data are exponential with rate lambda; the mean is 1/lambda. Given x-bar = 4, find the MoM estimate of lambda.
Answer. lambda-hat = 1/x-bar = 0.25.
Worked Example 2
Problem. For a Poisson(lambda), E[X] = lambda. A sample has counts with mean 3.2. Give the MoM estimate and note its relation to the MLE.
Answer. lambda-hat = 3.2; for the Poisson, MoM and MLE give the same estimate, x-bar.
Worked Example 3
Problem. A Uniform(0, theta) has E[X] = theta/2. With x-bar = 5, find the MoM estimate, and note a flaw if the sample max exceeds it.
Answer. theta-hat = 2*x-bar = 10, but MoM can yield an estimate below the observed maximum, an impossibility the MLE (sample max) avoids.
Problem. A Geometric distribution (number of trials to first success) has mean 1/p. A sample of waiting times has mean x-bar = 5. Find the MoM estimate of p, and state whether it lies in a valid range.
Solution. The first moment is E[X] = 1/p. Setting x-bar = 1/p gives 5 = 1/p, so p-hat = 1/x-bar = 1/5 = 0.2. Since 0 < 0.2 < 1, the estimate is a valid probability. In Python: 1/np.mean(data) = 0.2.
The likelihood reinterprets the data's joint density as a function of the parameter, scoring how well each parameter explains the data.
For i.i.d. data x_1,...,x_n from a model with density f(x; theta), the likelihood is L(theta) = product f(x_i; theta) — the same joint density, but viewed as a function of theta with the data fixed. It answers: how plausible is each theta given what we observed? Because products of many small numbers underflow and are awkward to differentiate, we take the log-likelihood ell(theta) = sum log f(x_i; theta), turning the product into a sum. The log is monotonic, so it has the same maximizer as L. The log-likelihood is the central object linking estimation, hypothesis testing (likelihood-ratio tests), and Bayesian inference (where it multiplies the prior). Mastering how to write down a log-likelihood for a model is the key skill that unlocks MLE.
Worked Example 1
Problem. Write the log-likelihood for n i.i.d. Bernoulli(p) observations with s successes.
Answer. ell(p) = s*log(p) + (n - s)*log(1 - p).
Worked Example 2
Problem. For n i.i.d. Poisson(lambda) counts with sum x_i = T, write the log-likelihood (dropping constants).
Answer. ell(lambda) = -n*lambda + T*log(lambda) - const, where T = sum x_i.
Worked Example 3
Problem. Two values of p give L(0.4) = 0.012 and L(0.6) = 0.030 for some Bernoulli data. Which p is more likely, and why compare log-likelihoods instead?
Answer. p = 0.6 is more likely (higher L and higher log-L); logs are preferred for numerical stability and easier differentiation.
Problem. For n i.i.d. Exponential(lambda) observations with density f(x;lambda)=lambda*e^{-lambda x}, write the log-likelihood in terms of n and S = sum x_i.
Solution. L(lambda) = product lambda*e^{-lambda x_i} = lambda^n * e^{-lambda * sum x_i} = lambda^n e^{-lambda S}. Taking the log: ell(lambda) = n*log(lambda) - lambda*S. In Python: ell = lambda lam: n*np.log(lam) - lam*S.
The MLE is the parameter that maximizes the (log-)likelihood, found by setting its derivative to zero.
Maximum likelihood estimation chooses theta-hat to maximize the log-likelihood ell(theta) — the value under which the observed data are most probable. Analytically, you differentiate ell with respect to theta, set the score equation ell'(theta) = 0, solve, and confirm it is a maximum (second derivative negative). MLEs have excellent large-sample properties: they are consistent, asymptotically unbiased, asymptotically normal, and efficient (achieving the smallest possible variance, the Cramer-Rao bound). The curvature of the log-likelihood at the maximum (the observed information) yields standard errors. These properties are why MLE underlies logistic regression, GLMs, and most modern statistical models. The recipe is always the same: write ell, differentiate, solve the score equation, verify a maximum.
Worked Example 1
Problem. Derive the MLE of p for Bernoulli data with s successes in n trials.
Answer. p-hat = s/n, the sample proportion of successes.
Worked Example 2
Problem. Derive the MLE of lambda for n i.i.d. Exponential(lambda) data with S = sum x_i.
Answer. lambda-hat = n/S = 1/x-bar, the reciprocal of the sample mean.
Worked Example 3
Problem. Derive the MLE of mu for n i.i.d. Normal(mu, sigma^2) with sigma^2 known.
Answer. mu-hat = x-bar; the MLE of a normal mean is the sample mean.
Problem. Derive the MLE of lambda for n i.i.d. Poisson(lambda) counts with T = sum x_i.
Solution. From ell(lambda) = -n*lambda + T*log(lambda) - const, differentiate: ell'(lambda) = -n + T/lambda. Set to 0: T/lambda = n -> lambda-hat = T/n = x-bar. The second derivative -T/lambda^2 < 0 confirms a maximum. So the Poisson MLE is the sample mean: lambda-hat = x-bar. In Python: np.mean(counts).
When no closed form exists, maximize the log-likelihood numerically with an optimizer.
Many models (logistic regression, custom mixtures) have no closed-form MLE, so we maximize ell(theta) numerically. The standard trick is to minimize the negative log-likelihood (NLL), since optimizers minimize by default. Tools like scipy.optimize.minimize or R's optim use gradient-based methods (e.g., BFGS) that iterate from a starting guess toward the optimum. Good practice includes: providing sensible starting values (method-of-moments estimates work well), optimizing on an unconstrained scale (e.g., optimize log(lambda) to keep lambda > 0), and checking convergence flags. The Hessian (second derivatives) at the optimum estimates the information matrix, whose inverse gives the parameter covariance and hence standard errors. This numerical approach is how virtually all real statistical and ML models are actually fit.
Worked Example 1
Problem. Write Python to numerically find the MLE rate of exponential data and confirm it matches 1/x-bar.
Answer. The optimizer's lambda-hat equals 1/x-bar, confirming the closed-form MLE numerically.
Worked Example 2
Problem. Why optimize log(lambda) instead of lambda directly, and how do you recover lambda?
Answer. Optimizing on the log scale keeps lambda > 0; recover lambda-hat = exp(theta-hat).
Worked Example 3
Problem. After fitting by MLE, how do you get a standard error for the estimate using the Hessian?
Answer. Invert the Hessian of the NLL at the optimum; the square roots of its diagonal are the parameter standard errors.
Problem. Outline Python code to find the MLE of the Poisson rate by numerical optimization and verify it equals x-bar.
Solution. import numpy as np; from scipy.optimize import minimize. Define the negative log-likelihood (dropping the constant log factorial term): nll = lambda lam, x: -(-len(x)*lam[0] + x.sum()*np.log(lam[0])). Then res = minimize(nll, x0=[1.0], args=(counts,), method='Nelder-Mead'). Finally print(res.x[0], np.mean(counts)); the optimizer's estimate equals the sample mean x-bar, matching the closed-form Poisson MLE.
Simulate (or load) data you believe is exponentially or Poisson distributed. Derive the MLE of the rate parameter by hand, then implement it in Python or R and compare it to the method-of-moments estimate. Finally, use a numerical optimizer (scipy.optimize or R's optim) to maximize the log-likelihood directly and confirm it matches your closed-form MLE.
Deliverable · A short report showing the hand-derived MLE, the method-of-moments estimate, and the optimizer result side by side, with a note on whether the estimates agree and why.
1. An estimator is unbiased if:
Answer B. Unbiasedness means E[estimator] equals the parameter. Convergence as n grows is consistency, a different property.
2. Mean squared error of an estimator decomposes into:
Answer C. MSE = Var(estimator) + (Bias)^2, which is the basis of the bias-variance tradeoff.
3. Why is the log-likelihood maximized instead of the likelihood itself?
Answer B. Because log is monotonic, the maximizer is unchanged; taking logs converts a product of densities into a sum, avoiding underflow and simplifying derivatives.
4. For n i.i.d. observations from a Bernoulli(p) distribution, the MLE of p is:
Answer B. Maximizing the Bernoulli log-likelihood gives p-hat equal to the number of successes divided by n, i.e., the sample proportion.
5. The method of moments estimates parameters by:
Answer B. Method-of-moments matches sample moments (e.g., the sample mean) to their model expressions and solves for the parameters.
I can describe an estimator's bias, variance, and consistency and use mean squared error to compare estimators.
I can write down a likelihood function for a model and find the maximum likelihood estimate analytically or numerically.
I can apply the method of moments as an alternative estimation strategy.
A confidence interval is a procedure with a long-run coverage guarantee, not a probability about one fixed interval.
A confidence interval (CI) is a range computed from data designed so that, if the sampling and computation were repeated many times, a stated proportion (e.g., 95%) of the resulting intervals would contain the true parameter. The confidence level describes the reliability of the method, not the probability that the one interval you computed contains the parameter — once computed, that interval either does or does not contain the fixed parameter, so the probability is 0 or 1. The correct frequentist statement is about the procedure across repetitions. This distinction matters in data science because stakeholders routinely misread '95% CI' as 'a 95% chance the truth is in this range,' which is a Bayesian-style claim that requires a prior. Coverage is the property we actually want: the realized capture rate should match the nominal level.
Worked Example 1
Problem. A 95% CI for a mean is [12.1, 15.3]. Is the statement 'there is a 95% probability mu is between 12.1 and 15.3' correct?
Answer. No; the 95% describes the method's long-run coverage, not the probability that this specific interval contains mu.
Worked Example 2
Problem. You simulate 1000 datasets and build a 95% CI for the known true mean each time. Roughly how many should contain the truth, and what does that demonstrate?
Answer. About 950 of 1000 intervals should contain the true mean, demonstrating the 95% coverage of the procedure.
Worked Example 3
Problem. A team reports 'we are 95% confident the click rate is between 4.8% and 5.6%.' Rewrite this so it is technically accurate but still readable.
Answer. 'Using a procedure that captures the true rate in 95% of repeated samples, this sample's interval is 4.8% to 5.6%' — accurate and clear.
Problem. An analyst computes a 90% CI of [0.30, 0.38] for a proportion and writes 'p has a 90% chance of being in [0.30, 0.38].' State why this is technically wrong and give a correct one-sentence interpretation.
Solution. It is wrong because p is a fixed unknown constant and the computed interval is also fixed, so p is either in it or not (probability 0 or 1); the 90% is a property of the procedure, not this interval. Correct interpretation: 'If we repeated the sampling and interval construction many times, about 90% of the intervals produced would contain the true proportion p, and this particular interval is [0.30, 0.38].'
Build a mean interval as estimate plus-or-minus a critical value times the standard error, using t when sigma is unknown.
A confidence interval for a mean is x-bar +/- (critical value) * SE. When the population sigma is known (rare), the critical value is a z from the normal (1.96 for 95%) and SE = sigma/sqrt(n). When sigma is unknown (the usual case), we estimate it with s and use the t-distribution with n - 1 degrees of freedom, whose heavier tails widen the interval to account for estimating sigma. As n grows the t critical value approaches the z value, so the two methods converge. The interval requires either approximately normal data or a large enough n for the CLT to make x-bar approximately normal. This is the workhorse interval of applied statistics: reporting a mean estimate with a t-based interval is the default way to quantify uncertainty about a population average.
Worked Example 1
Problem. n = 25, x-bar = 50, s = 10. Build a 95% t-interval for the mean.
Answer. 95% CI = [45.87, 54.13].
Worked Example 2
Problem. With sigma known = 10, n = 25, x-bar = 50, build the 95% z-interval and compare to the t-interval above.
Answer. z-interval = [46.08, 53.92], slightly narrower than the t-interval because the known sigma needs no extra-uncertainty correction.
Worked Example 3
Problem. A sample of n = 100 has x-bar = 8.4, s = 3. Build a 99% confidence interval for the mean.
Answer. 99% CI = [7.61, 9.19].
Problem. A sample of n = 16 measurements has x-bar = 100 and s = 8 (sigma unknown). Construct a 95% confidence interval for the mean.
Solution. SE = s/sqrt(n) = 8/sqrt(16) = 8/4 = 2. df = 16 - 1 = 15; t* for 95% two-sided = 2.131 (scipy: stats.t.ppf(0.975, 15) -> 2.1314). Margin = 2.131 * 2 = 4.262. Interval = 100 +/- 4.262 = [95.74, 104.26].
A proportion interval is p-hat plus-or-minus z times the standard error of a proportion.
For a binary outcome, the sample proportion p-hat = s/n estimates the population proportion p. Its standard error is SE = sqrt(p-hat*(1 - p-hat)/n), and the classic (Wald) 95% interval is p-hat +/- 1.96*SE. This relies on a normal approximation that works well when both n*p-hat and n*(1 - p-hat) are at least about 10; for small samples or proportions near 0 or 1 it behaves poorly (intervals can fall outside [0,1] or undercover). Better alternatives are the Wilson score interval and the Agresti-Coull adjustment, which add pseudo-counts for more reliable coverage. Proportion intervals are everywhere in data science: conversion rates, click-through rates, defect rates, and poll results all need them.
Worked Example 1
Problem. Out of n = 400 visitors, 80 converted. Build a 95% Wald interval for the conversion rate.
Answer. 95% CI ~ [0.161, 0.239], i.e., about 16.1% to 23.9%.
Worked Example 2
Problem. Check the normal-approximation condition for the previous example.
Answer. Both success and failure counts (80 and 320) exceed 10, so the normal approximation is valid.
Worked Example 3
Problem. Only 1 of 20 trials succeeded. Why is the Wald interval unreliable, and what should you use?
Answer. The Wald interval gives an impossible negative bound; use the Wilson interval (about [0.009, 0.236]) for small-sample reliability.
Problem. A poll of n = 1000 finds 540 in favor. Build a 95% Wald confidence interval for the true proportion and verify the normal-approximation condition.
Solution. p-hat = 540/1000 = 0.54. SE = sqrt(0.54*0.46/1000) = sqrt(0.2484/1000) = sqrt(0.0002484) = 0.01576. Margin = 1.96 * 0.01576 = 0.0309. Interval = 0.54 +/- 0.0309 = [0.509, 0.571]. Condition: n*p-hat = 540 and n*(1 - p-hat) = 460, both well above 10, so the approximation is valid. In statsmodels: proportion_confint(540, 1000, method='normal').
The margin of error links confidence, variability, and sample size; solving for n lets you plan a study.
The margin of error (ME) is the half-width of a confidence interval: ME = critical value * SE. Three levers control it. Higher confidence raises the critical value (z grows from 1.645 at 90% to 1.96 at 95% to 2.576 at 99%), widening the interval. Greater variability (larger sigma or p near 0.5) widens it. Larger sample size shrinks SE like 1/sqrt(n), narrowing it. Inverting the formula lets you plan: to achieve a target ME, solve for n. For a mean, n = (z*sigma/ME)^2; for a proportion, n = z^2 * p*(1 - p)/ME^2, using p = 0.5 for the most conservative (largest) n when p is unknown. This is the core calculation behind survey design and experiment sizing.
Worked Example 1
Problem. Estimate a mean within ME = 1 at 95% confidence, with sigma = 5. Find the required n.
Answer. n = 97 (rounding 96.04 up).
Worked Example 2
Problem. Estimate a proportion within ME = 0.03 at 95% with no prior guess for p. Find the conservative n.
Answer. n = 1068 (the familiar '~1067 for +/-3%' poll size).
Worked Example 3
Problem. For a fixed sample, how does the margin of error change going from 95% to 99% confidence? Use z = 1.96 vs 2.576.
Answer. The margin grows by a factor of ~1.31 (about 31% wider) when moving from 95% to 99% confidence.
Problem. You want to estimate a population mean within a margin of error of 2 at 99% confidence, and a pilot suggests sigma = 12. How large must the sample be?
Solution. Use n = (z*sigma/ME)^2 with z = 2.576 for 99%. z*sigma/ME = 2.576*12/2 = 30.912/2 = 15.456. n = 15.456^2 = 238.9. Round up to n = 239. In Python: import math; math.ceil((2.576*12/2)**2) -> 239.
The bootstrap treats your sample as the population and resamples it to approximate a statistic's sampling distribution.
The bootstrap is a resampling method that estimates the sampling distribution of almost any statistic without formulas or distributional assumptions. The idea: your observed sample is the best available stand-in for the population, so repeatedly draw new samples of the same size n from it with replacement, compute the statistic on each resample, and collect those values. The spread of this bootstrap distribution estimates the statistic's standard error, and its quantiles can form confidence intervals. The bootstrap shines for statistics with no neat formula (the median, a correlation, a ratio, a trimmed mean) and for small or non-normal data. Its main assumptions are that the sample is representative and observations are independent. In data science it is a go-to tool for honest uncertainty when classical theory is hard or unavailable.
Worked Example 1
Problem. Describe one bootstrap resample of the data [3, 7, 7, 10] and why duplicates appear.
Answer. One resample might be [7, 3, 7, 7] with mean 6.0; with-replacement sampling naturally produces repeats and omissions.
Worked Example 2
Problem. Write Python to bootstrap the standard error of the median for a dataset.
Answer. The standard deviation of 10,000 bootstrap medians estimates the SE of the median directly from data.
Worked Example 3
Problem. Roughly what fraction of the original observations appear in a typical large-n bootstrap resample?
Answer. About 63.2% of the distinct original observations appear in a typical resample (the rest are omitted; some appear multiple times).
Problem. Outline a bootstrap to estimate the standard error of the mean for a sample of 200 values, and state how many resamples you would use.
Solution. Treat the 200 values as the population. For B = 10,000 iterations, draw a resample of size 200 with replacement and record its mean. The standard deviation of those 10,000 bootstrap means estimates the SE of the mean. In Python: import numpy as np; rng=np.random.default_rng(0); boot=[rng.choice(data,200,replace=True).mean() for _ in range(10000)]; se=np.std(boot,ddof=1). Using B = 10,000 keeps Monte Carlo noise small.
Read interval bounds directly off the bootstrap distribution's quantiles, with refinements when needed.
Once you have a bootstrap distribution of a statistic, the percentile method forms a (1 - alpha) interval by taking its alpha/2 and 1 - alpha/2 quantiles — for 95%, the 2.5th and 97.5th percentiles. It is simple and assumption-light, but it can have biased coverage when the statistic's distribution is skewed or the estimator is biased. Refinements address this: the basic (pivotal) bootstrap reflects the percentiles around the observed estimate, and the bias-corrected and accelerated (BCa) interval adjusts for both bias and skewness, generally giving the most accurate coverage. A normal-approximation bootstrap interval uses estimate +/- z * (bootstrap SE). Data scientists default to percentile or BCa intervals for statistics where classical formulas are unavailable or unreliable.
Worked Example 1
Problem. From 10,000 bootstrap means, the 2.5th percentile is 48.2 and the 97.5th is 53.6. Give the 95% percentile interval.
Answer. 95% percentile bootstrap CI = [48.2, 53.6].
Worked Example 2
Problem. The observed mean is 50 with bootstrap SE 1.4. Give a normal-approximation bootstrap 95% interval and compare to the percentile one above.
Answer. Normal-approx CI = [47.26, 52.74]; its forced symmetry differs from the skewed percentile interval, which better reflects the bootstrap shape.
Worked Example 3
Problem. When the bootstrap distribution is skewed and the estimator biased, which interval is preferred and why?
Answer. Use the BCa interval; it corrects for bias and skewness and gives more accurate coverage than the plain percentile method.
Problem. You bootstrap a correlation coefficient 5000 times; the 2.5th percentile of the bootstrap values is 0.31 and the 97.5th is 0.67. State the 95% percentile CI and note one situation where you would switch to BCa.
Solution. The 95% percentile bootstrap CI is just the [2.5th, 97.5th] percentiles: [0.31, 0.67]. In Python: np.percentile(boot_corrs, [2.5, 97.5]) -> [0.31, 0.67]. You would switch to BCa if the bootstrap distribution of the correlation is noticeably skewed (common for correlations near +/-1) or the estimator is biased, because BCa adjusts the endpoints for bias and skewness to restore proper coverage.
Take a sample of continuous data (real or simulated) and compute a 95% confidence interval for the mean using the t-distribution. Then write a bootstrap routine in Python or R that resamples the data thousands of times and builds a percentile confidence interval for the same mean. Compare the two intervals and discuss when the bootstrap is preferable.
Deliverable · A script that prints both intervals plus a histogram of the bootstrap distribution with the interval bounds marked, and a short comparison of the two approaches.
1. Which interpretation of a 95% confidence interval is correct?
Answer C. Confidence refers to the long-run procedure: across many repetitions, 95% of such intervals capture the parameter. A single computed interval either does or does not contain it.
2. When constructing a confidence interval for a mean with an unknown population standard deviation, you use:
Answer A. With sigma unknown and estimated from the sample, the t-distribution accounts for the extra uncertainty via its heavier tails.
3. Holding everything else fixed, increasing the confidence level from 90% to 99% will:
Answer B. Higher confidence requires a larger critical value, increasing the margin of error and widening the interval.
4. The bootstrap approximates a statistic's sampling distribution by:
Answer B. The bootstrap treats the sample as a stand-in for the population and resamples it with replacement many times to mimic sampling variability.
5. A percentile bootstrap 95% interval is formed by:
Answer B. The percentile method reads the central 95% of the bootstrap distribution, i.e., its 2.5th and 97.5th percentiles, as the interval bounds.
I can construct and correctly interpret a confidence interval for a mean or proportion.
I can use the bootstrap to estimate the sampling distribution of a statistic and build a confidence interval without strong assumptions.
I can reason about how confidence level, variability, and sample size affect the margin of error.
Testing frames a question as a default 'no effect' claim that the data must argue against.
A hypothesis test starts with a null hypothesis H0 stating no effect or no difference (e.g., mu = 0, or two means equal), and an alternative H1 stating the effect we suspect. The logic is proof by contradiction under uncertainty: we assume H0 is true, compute a test statistic measuring how far the data deviate from H0 in standard-error units, and ask how surprising that deviation would be if H0 held. If the data are sufficiently unlikely under H0, we reject H0 in favor of H1; otherwise we fail to reject (we never 'accept' H0). The alternative can be two-sided (different) or one-sided (greater or less). Framing the right H0/H1 — and choosing one- vs. two-sided before seeing data — is the foundation that keeps every downstream test honest.
Worked Example 1
Problem. A coin is suspected of being unfair. State H0 and H1 and whether the test is one- or two-sided.
Answer. H0: p = 0.5; H1: p != 0.5; two-sided test.
Worked Example 2
Problem. A new drug is claimed to lower blood pressure versus placebo. Write directional hypotheses for the mean difference d = mu_drug - mu_placebo.
Answer. H0: d = 0; H1: d < 0; a one-sided (left-tailed) test.
Worked Example 3
Problem. A test yields a result that does not reach significance. Is it correct to conclude 'H0 is true'?
Answer. No; we 'fail to reject H0' — the data simply do not provide enough evidence against it, which is not proof H0 is true.
Problem. A company believes a redesign changed (in either direction) the average time users spend on a page, formerly mu = 120 seconds. State the hypotheses and the test sidedness.
Solution. Let mu be the new mean time on page. H0: mu = 120 (no change from the old average). H1: mu != 120 (the redesign changed the time, up or down). Because the claim is 'changed in either direction,' this is a two-sided test. We would assume H0, compute a t-statistic from the new sample, and reject H0 only if the data are sufficiently unlikely under mu = 120.
A p-value measures surprise under the null; alpha sets the threshold and equals the Type I error rate.
The p-value is the probability, computed assuming H0 is true, of observing a test statistic at least as extreme as the one obtained. A small p-value means the data would be surprising under H0, casting doubt on it. The significance level alpha (commonly 0.05) is a pre-chosen threshold: reject H0 if p <= alpha. Two errors are possible. A Type I error rejects a true null (false positive); its probability is exactly alpha. A Type II error fails to reject a false null (false negative); its probability is beta, and power = 1 - beta. There is a tradeoff: lowering alpha reduces false positives but raises beta (lowers power) unless n increases. A p-value is not the probability that H0 is true and not a measure of effect size — always report an effect estimate alongside it.
Worked Example 1
Problem. A two-sided test gives a z-statistic of 2.1. Find the p-value and decide at alpha = 0.05.
Answer. p ~ 0.036 <= 0.05, so reject H0.
Worked Example 2
Problem. Interpret p = 0.20 correctly and state which error you risk if you fail to reject.
Answer. p = 0.20 is the chance of equally/more extreme data under H0; failing to reject risks a Type II error if H0 is really false.
Worked Example 3
Problem. Explain how lowering alpha from 0.05 to 0.01 affects the two error types, holding n fixed.
Answer. Lower alpha cuts Type I errors but raises Type II errors (less power), unless you increase n to compensate.
Problem. A one-sided (upper-tail) test produces a z-statistic of 1.5. Compute the p-value and decide at alpha = 0.05, then name the error you could be committing.
Solution. One-sided upper-tail p = P(Z >= 1.5) = 1 - 0.9332 = 0.0668 (scipy: 1 - stats.norm.cdf(1.5) = 0.0668). Since 0.0668 > 0.05, we fail to reject H0. If the true effect actually exists (H0 false), failing to reject it is a Type II error (a false negative).
t-tests compare means using t = (estimate - null) / standard error, judged against the t-distribution.
A t-test assesses claims about means when sigma is unknown. The one-sample t-test checks whether a mean equals a hypothesized value: t = (x-bar - mu0)/(s/sqrt(n)) with n - 1 degrees of freedom. The paired t-test applies the one-sample test to within-pair differences (e.g., before/after on the same subjects). The two-sample (independent) t-test compares two group means; Welch's version, which does not assume equal variances, is the safe default and uses t = (x1-bar - x2-bar)/sqrt(s1^2/n1 + s2^2/n2). You compare the t-statistic to a t-distribution to get a p-value. Assumptions are approximately normal data (or large n via the CLT) and independence; for paired designs, the pairs must be matched. t-tests are the most common tool for comparing averages in experiments and A/B analyses.
Worked Example 1
Problem. n = 16, x-bar = 52, s = 8, test H0: mu = 50 vs H1: mu != 50.
Answer. t = 1.0 on 15 df, p ~ 0.33; fail to reject H0 — no evidence mu differs from 50.
Worked Example 2
Problem. Paired before/after differences have mean d-bar = 3, s_d = 5, n = 25. Test whether the change is nonzero.
Answer. t = 3.0 on 24 df, p ~ 0.006; reject H0 — there is a significant mean change.
Worked Example 3
Problem. Two groups: group 1 (n1=30, x1-bar=10.0, s1=2.0), group 2 (n2=30, x2-bar=11.5, s2=2.5). Run a Welch two-sample test.
Answer. Welch t ~ -2.57, p ~ 0.013; reject H0 — group 2's mean is significantly higher.
Problem. A sample of n = 9 has x-bar = 104 and s = 6. Test H0: mu = 100 against H1: mu != 100 at alpha = 0.05.
Solution. SE = s/sqrt(n) = 6/sqrt(9) = 6/3 = 2. t = (104 - 100)/2 = 2.0 on df = 8. Two-sided p = 2*P(T8 >= 2.0); scipy: 2*(1-stats.t.cdf(2.0, 8)) = 0.0805. Since 0.0805 > 0.05, we fail to reject H0 — there is not enough evidence that mu differs from 100 (note the critical t for 8 df is 2.306, and 2.0 < 2.306).
The chi-squared test compares observed counts to counts expected under a hypothesis about categories.
The chi-squared (chi^2) test works with categorical counts. The goodness-of-fit test checks whether one categorical variable's observed frequencies match a hypothesized distribution; the test of independence checks whether two categorical variables in a contingency table are associated. Both use the statistic chi^2 = sum (O - E)^2 / E, where O are observed counts and E are counts expected under H0. For independence, E for a cell = (row total * column total)/grand total. The statistic follows a chi-squared distribution; degrees of freedom are (categories - 1) for goodness of fit and (rows - 1)*(cols - 1) for independence. A large chi^2 (small p) means observed counts deviate from expectation more than chance allows. A rule of thumb requires expected counts of at least about 5 per cell for the approximation to hold.
Worked Example 1
Problem. A die is rolled 60 times with counts 8,12,9,11,10,10. Test goodness of fit to a fair die.
Answer. chi^2 = 1.0 on 5 df, p ~ 0.96; no evidence the die is unfair.
Worked Example 2
Problem. Compute the expected count for a contingency cell whose row total is 80, column total is 50, and grand total is 200.
Answer. Expected count E = 20.
Worked Example 3
Problem. A 2x2 table of (treatment vs control) by (success vs failure) is [[30,20],[20,30]]. Outline the independence test and its df.
Answer. chi^2 = 4.0 on 1 df, p ~ 0.046; borderline-significant association (use chi2_contingency, which by default applies Yates' continuity correction for 2x2 tables).
Problem. A survey expects responses split 50/30/20 across three options, and out of 100 people observes 40/35/25. Run a goodness-of-fit test.
Solution. Expected counts for n = 100: 50, 30, 20. chi^2 = sum (O-E)^2/E = (40-50)^2/50 + (35-30)^2/30 + (25-20)^2/20 = 100/50 + 25/30 + 25/20 = 2.0 + 0.833 + 1.25 = 4.083. df = 3 - 1 = 2; p = P(chi^2_2 >= 4.083); scipy: 1 - stats.chi2.cdf(4.083, 2) = 0.130. Since 0.130 > 0.05, fail to reject — observed counts are consistent with the 50/30/20 split.
ANOVA uses an F-test to ask whether three or more group means differ, comparing between- to within-group variation.
Analysis of variance (ANOVA) tests H0 that several group means are all equal against the alternative that at least one differs. Rather than running many pairwise t-tests (which inflates the error rate), one-way ANOVA partitions total variation into between-group (variation of group means around the grand mean) and within-group (variation inside each group). The F-statistic is the ratio of mean squares: F = MSB/MSW, where MSB = SSB/(k - 1) and MSW = SSW/(N - k), with k groups and N total observations. Under H0, F follows an F-distribution with (k - 1, N - k) degrees of freedom, centered near 1; a large F (small p) signals that between-group differences exceed what within-group noise explains. ANOVA assumes independent, roughly normal groups with similar variances. A significant F says 'some means differ' but not which — follow up with post-hoc comparisons.
Worked Example 1
Problem. Three groups give SSB = 60 (k = 3) and SSW = 90 (N = 30). Compute the F-statistic.
Answer. F = 9.0 on (2, 27) df, p ~ 0.001; reject H0 — at least one group mean differs.
Worked Example 2
Problem. Why does a large F lead to rejecting equal means? Interpret F near 1 vs F much greater than 1.
Answer. F near 1 is consistent with equal means; F much greater than 1 means between-group spread exceeds within-group noise, so we reject H0.
Worked Example 3
Problem. Outline running a one-way ANOVA in Python on three samples and what a significant result does not tell you.
Answer. stats.f_oneway gives F and p; a significant result flags that some means differ but requires a post-hoc test (e.g., Tukey HSD) to identify which.
Problem. Four groups (k = 4) with N = 40 observations give SSB = 120 and SSW = 144. Compute MSB, MSW, and the F-statistic, then state the df.
Solution. MSB = SSB/(k - 1) = 120/(4 - 1) = 120/3 = 40. MSW = SSW/(N - k) = 144/(40 - 4) = 144/36 = 4. F = MSB/MSW = 40/4 = 10.0. Degrees of freedom are (k - 1, N - k) = (3, 36). The p-value is P(F_{3,36} >= 10.0); scipy: 1 - stats.f.cdf(10.0, 3, 36) ~ 0.00006, so we strongly reject H0 that all four means are equal.
Running many tests inflates false positives; corrections and pre-registration keep conclusions honest.
Each test at level alpha has an alpha chance of a false positive. Run many and the family-wise error rate (FWER), the chance of at least one false positive, balloons: with m independent true nulls it is 1 - (1 - alpha)^m. p-hacking — trying many analyses, subgroups, or metrics and reporting only the significant ones — exploits this and produces irreproducible findings. Corrections restore control. The Bonferroni method tests each hypothesis at alpha/m, controlling FWER but conservatively. The Benjamini-Hochberg procedure controls the false discovery rate (FDR), the expected proportion of false positives among rejections, offering more power for large m (common in genomics or many-metric experiments). The cure is also procedural: pre-register hypotheses and the primary metric before looking at data so the analysis is not contorted to find significance.
Worked Example 1
Problem. You run 10 independent tests at alpha = 0.05 with all nulls true. What is the chance of at least one false positive?
Answer. About 40% — there is a 0.40 chance of at least one false positive across 10 tests.
Worked Example 2
Problem. Apply a Bonferroni correction for 10 tests at family alpha = 0.05. What is the per-test threshold, and which of p = 0.004, 0.02, 0.06 survive?
Answer. Per-test threshold = 0.005; only p = 0.004 stays significant.
Worked Example 3
Problem. Explain why reporting only the 'best' of 20 metrics tested is p-hacking, and one fix.
Answer. It is p-hacking because selecting the best of 20 tests hides the multiplicity; fix it by pre-registering the primary metric or correcting across all tests.
Problem. An experiment tests 5 independent metrics at alpha = 0.05, all nulls true. (a) Find the family-wise false-positive probability. (b) Give the Bonferroni per-test threshold. (c) Decide which of p = 0.008, 0.03, 0.04 are significant after Bonferroni.
Solution. (a) P(at least one false positive) = 1 - (1 - 0.05)^5 = 1 - 0.95^5 = 1 - 0.7738 = 0.2262, about 23%. (b) Bonferroni threshold = 0.05/5 = 0.01. (c) Compare each p to 0.01: 0.008 <= 0.01 is significant; 0.03 and 0.04 exceed 0.01 and are not. Only p = 0.008 survives the correction.
You are given three mini-datasets: two groups of continuous measurements, a contingency table of categorical counts, and several groups for a one-way comparison. Using Python (SciPy/statsmodels) or R, pick and run the correct test for each — a two-sample t-test, a chi-squared test, and an ANOVA — and report the test statistic, p-value, and a plain-language conclusion. Where you run several tests, apply a multiple-comparison correction and explain its effect.
Deliverable · A report with the chosen test for each dataset, the statistics and p-values, the corrected results, and a sentence each on what the conclusion does and does not justify.
1. A p-value of 0.03 means:
Answer B. A p-value is computed under the assumption that the null is true; it is the probability of equally or more extreme data, not the probability of any hypothesis.
2. A Type I error occurs when you:
Answer B. A Type I error is a false positive: rejecting a null that is actually true. Failing to reject a false null is a Type II error.
3. To test whether two categorical variables in a contingency table are independent, you use:
Answer B. The chi-squared test of independence compares observed cell counts to those expected under independence of the two categorical variables.
4. ANOVA is most appropriate when you want to:
Answer A. One-way ANOVA uses an F-test to compare several group means simultaneously, controlling the overall error rate better than many pairwise t-tests.
5. Running 20 independent tests at alpha = 0.05 with all nulls true, you would expect roughly how many false positives?
Answer B. At a 5% rate, about 1 in 20 true-null tests will be falsely significant by chance, which is why multiple-comparison corrections are needed.
I can state null and alternative hypotheses and interpret a p-value and significance level correctly.
I can choose and run the appropriate test (t-test, chi-squared, or ANOVA) for a given question and report the result responsibly.
I can explain Type I and Type II errors and why multiple comparisons inflate false positives.
Random assignment is what turns a correlation into a credible causal claim.
Observational data can show that two things move together, but correlation alone cannot establish causation because of confounders — variables that influence both the treatment and the outcome. A randomized controlled experiment breaks this by randomly assigning units to a treatment group and a control group. Randomization balances all confounders (known and unknown) across groups in expectation, so the two groups are comparable except for the treatment. Any systematic difference in outcomes can then be attributed to the treatment, estimating the average treatment effect (ATE) = mean(outcome | treated) - mean(outcome | control). The control group provides the counterfactual baseline of what would have happened without the change. This logic — randomize, compare, attribute — is the gold standard for causal claims and the backbone of A/B testing in data science.
Worked Example 1
Problem. An observational study finds users of a feature have higher retention. Why can't we conclude the feature causes retention?
Answer. Self-selection means confounders (like baseline engagement) are not balanced, so the association is not necessarily causal.
Worked Example 2
Problem. A randomized test gives treated mean conversion 0.062 and control mean 0.050. Estimate the average treatment effect.
Answer. Estimated ATE = +0.012 (1.2 percentage points absolute, a 24% relative lift).
Worked Example 3
Problem. How does randomization handle an unknown confounder that an analyst never even measures?
Answer. Randomization balances even unmeasured confounders in expectation, which observational adjustment cannot guarantee.
Problem. An online store randomly shows half its visitors a new checkout flow. Treated visitors convert at 8.0% and control at 7.0%. (a) Estimate the absolute and relative treatment effect. (b) Explain in one sentence why randomization makes this causal.
Solution. (a) Absolute ATE = 8.0% - 7.0% = 1.0 percentage point; relative lift = 1.0/7.0 = 0.143, about a 14.3% increase. (b) Because visitors were randomly assigned, the treatment and control groups are comparable in expectation on all other factors (device, intent, traffic source, etc.), so the 1.0-point difference can be attributed to the new checkout flow rather than to a confounder.
A sound A/B test fixes a single primary metric, clear hypotheses, and a consistent randomization unit before launch.
Designing an A/B test means specifying three things up front. The metric: one primary success metric (e.g., conversion rate, revenue per user) chosen before launch, ideally with a few guardrail metrics to catch harm. The hypothesis: H0 of no difference between control (A) and treatment (B), with a directional or two-sided H1 and a pre-stated minimum detectable effect worth acting on. The randomization unit: what gets assigned to A or B — usually the user (not the page view), so a person sees a consistent experience and units are independent. Mismatched units (e.g., randomizing by session while measuring per-user metrics) break independence and bias inference. Defining these before collecting data prevents p-hacking and ensures the test answers a single, decision-relevant question.
Worked Example 1
Problem. A team wants to test a new homepage. Why randomize by user rather than by page view?
Answer. Randomize by user so each person sees one consistent variant and units are independent, avoiding correlated within-user observations.
Worked Example 2
Problem. State a complete primary hypothesis for testing whether variant B raises conversion above a 5% baseline.
Answer. Metric = per-user conversion; H0: p_B = p_A; H1: p_B != p_A; MDE = +1 point (0.05 -> 0.06), all fixed before launch.
Worked Example 3
Problem. Why include a guardrail metric like page-load time alongside the primary conversion metric?
Answer. Guardrail metrics catch unintended harm, so you do not ship a primary-metric win that degrades performance or another key outcome.
Problem. You want to test whether a new onboarding email increases 7-day activation. Specify the primary metric, the randomization unit, H0/H1, and one guardrail metric.
Solution. Primary metric: 7-day activation rate (activated users / enrolled users), measured per user. Randomization unit: the user, randomly assigned to receive the new email (B) or the old one (A), so each person gets one consistent experience. Hypotheses: H0: activation_B = activation_A; H1: activation_B > activation_A (one-sided, since we only care about an increase) with a pre-set MDE such as +2 percentage points. Guardrail metric: unsubscribe rate, monitored to ensure the new email does not drive users away even if activation rises.
Power is the probability of detecting a real effect; the MDE is the smallest effect a test can reliably catch.
Statistical power is the probability that a test correctly rejects H0 when a true effect of a given size exists; it equals 1 - beta, where beta is the Type II error rate. Power rises with larger true effect size, larger sample size, larger alpha, and smaller variance. Studies are typically designed for 80% power (beta = 0.20). The minimum detectable effect (MDE) is the smallest effect the test is powered to find at the chosen alpha and power: smaller MDEs demand much larger samples. There is a direct link via the standardized effect size (Cohen's d for means, or a standardized proportion difference). Computing power before launch prevents underpowered tests that waste resources and miss real effects, and it forces an honest conversation about what effect size is worth detecting.
Worked Example 1
Problem. A two-sample mean test has standardized effect d = 0.5, n = 64 per group, alpha = 0.05 two-sided. Roughly estimate the power.
Answer. Power is about 0.80 — roughly the conventional 80% target.
Worked Example 2
Problem. Holding alpha and effect size fixed, what happens to power as you double the sample size, and why?
Answer. Power increases because a larger n reduces the standard error, making a fixed true effect easier to detect.
Worked Example 3
Problem. Explain the tradeoff between MDE and sample size at fixed power and alpha.
Answer. Smaller MDEs require far larger samples — n scales like 1/MDE^2, so halving the MDE roughly quadruples the needed n.
Problem. A study is powered at 80% to detect a standardized effect of d = 0.4. If the team now wants to detect d = 0.2 at the same 80% power and alpha, by what factor must the per-group sample size change?
Solution. Required n per group scales like 1/d^2 (since n ~ 2*(z_alpha/2 + z_beta)^2 / d^2). Going from d = 0.4 to d = 0.2 halves d, so 1/d^2 changes by a factor of (0.4/0.2)^2 = 2^2 = 4. The per-group sample size must increase about fourfold. In statsmodels: TTestIndPower().solve_power(effect_size=0.2, power=0.8, alpha=0.05) returns roughly 4x the n for d = 0.4.
Translate the MDE, baseline, power, and alpha into the visitors and days a test needs.
Before launching, you compute how many units the test needs and how long that takes. For comparing two proportions, a common formula for n per group is n = (z_{alpha/2} + z_beta)^2 * (p1*(1 - p1) + p2*(1 - p2)) / (p2 - p1)^2, where p1 is the baseline and p2 = p1 + MDE. The z-values come from the chosen alpha (e.g., 1.96 for two-sided 0.05) and power (e.g., 0.84 for 80%). Duration follows from traffic: days = (total units needed)/(daily eligible units), rounded up and ideally spanning full weeks to cover weekday/weekend cycles. Tools like statsmodels' NormalIndPower or proportion effect-size helpers automate this. Planning duration prevents two failure modes: stopping too early (underpowered) and running needlessly long (wasting traffic and delaying decisions).
Worked Example 1
Problem. Baseline conversion p1 = 0.10, MDE to p2 = 0.12, alpha = 0.05 two-sided, power = 0.80. Estimate n per group.
Answer. About 3,834 users per group (roughly 7,700 total).
Worked Example 2
Problem. If each group needs 3,834 users and the site gets 2,000 eligible users/day split 50/50, how long must the test run?
Answer. Mathematically ~4 days, but run a full 7 days to cover weekly cycles.
Worked Example 3
Problem. Sketch statsmodels code to get the per-group sample size for detecting a lift from 0.10 to 0.12.
Answer. proportion_effectsize then NormalIndPower().solve_power returns about 3,800 users per group, matching the formula.
Problem. Baseline conversion p1 = 0.05, you want to detect a lift to p2 = 0.06 at alpha = 0.05 (two-sided) and 80% power. Estimate the per-group sample size.
Solution. (z_{alpha/2} + z_beta)^2 = (1.96 + 0.84)^2 = 2.8^2 = 7.84. Variance term = p1(1-p1) + p2(1-p2) = 0.05*0.95 + 0.06*0.94 = 0.0475 + 0.0564 = 0.1039. Denominator = (0.06 - 0.05)^2 = 0.01^2 = 0.0001. n = 7.84 * 0.1039 / 0.0001 = 0.8146 / 0.0001 = 8146 per group (round up). So roughly 8,150 users per variant, about 16,300 total.
Run to the planned sample, then compare the variants with the appropriate two-sample test and report effect plus uncertainty.
Once the test reaches its pre-computed sample size, analyze it with the test that matches the metric. For a conversion (proportion) metric, use a two-proportion z-test (or chi-squared on the 2x2 table). For a continuous metric like revenue per user, use a two-sample (Welch) t-test. Report not just the p-value but the estimated effect and its confidence interval, since a tiny but significant lift may not be worth shipping. Also check for a sample-ratio mismatch (the actual split deviating from the intended 50/50), which signals a logging or assignment bug that can invalidate results. Good analysis pairs a clear decision rule (ship if significant and the lift exceeds the MDE, respecting guardrails) with honest reporting of uncertainty. Avoid re-deciding based on data you peeked at mid-run.
Worked Example 1
Problem. Control: 500/10000 convert; Treatment: 560/10000 convert. Run a two-proportion z-test.
Answer. z ~ 1.89, p ~ 0.058; not significant at 0.05 — the 0.6-point lift is not yet distinguishable from chance.
Worked Example 2
Problem. The lift above is +0.006 with SE 0.003169. Give a 95% CI for the difference and interpret.
Answer. 95% CI for the lift ~ [-0.0002, 0.0122]; it includes 0, so the effect is not statistically distinguishable from no change.
Worked Example 3
Problem. You intended a 50/50 split but observe 10,000 vs 11,200 users. Why investigate before trusting the result?
Answer. The split deviates far more than chance (chi^2 ~ 68, p ~ 0), a sample-ratio mismatch indicating a bug; results are untrustworthy until resolved.
Problem. Control converts 300/5000 and treatment 360/5000. Compute both proportions, the pooled proportion, and the two-proportion z-statistic, then state the decision at alpha = 0.05.
Solution. p1 = 300/5000 = 0.060; p2 = 360/5000 = 0.072; pooled p = 660/10000 = 0.066. SE = sqrt(0.066*0.934*(1/5000 + 1/5000)) = sqrt(0.066*0.934*0.0004) = sqrt(2.4658e-5) = 0.004966. z = (0.072 - 0.060)/0.004966 = 0.012/0.004966 = 2.416. Two-sided p = 2*(1 - Phi(2.416)) ~ 0.0157. Since 0.0157 < 0.05, reject H0 — treatment significantly increases conversion (statsmodels: proportions_ztest([360,300],[5000,5000]) gives z~2.42, p~0.016).
Common traps quietly invalidate A/B tests; recognizing them protects your conclusions.
Several pitfalls undermine experiments. Peeking is repeatedly checking results and stopping as soon as significance appears; each look is another chance to cross the threshold by luck, so naive peeking can inflate the Type I error rate from 5% to well over 20%. The fix is to fix the sample size in advance or use sequential methods with adjusted boundaries (e.g., alpha-spending or always-valid p-values). Novelty (and primacy) effects mean a new variant temporarily changes behavior simply because it is new, so short tests can overstate or understate the durable effect; running long enough lets behavior settle. Simpson's paradox is when an effect within every subgroup reverses after pooling, caused by a confounder unevenly distributed across groups; segmenting and using proper weighting or stratification reveals the true picture. Awareness of these traps separates trustworthy experimentation from misleading dashboards.
Worked Example 1
Problem. An analyst checks an A/B test 5 times and stops at the first p < 0.05. Why does this inflate false positives?
Answer. Multiple peeks give multiple chances to cross 0.05 by luck, inflating the false-positive rate to roughly 20%+ instead of 5%.
Worked Example 2
Problem. A new feature shows a huge lift in week 1 that fades to near zero by week 4. What is happening and how do you handle it?
Answer. It is a novelty effect; run the experiment long enough for behavior to stabilize and judge by the durable effect, not the early spike.
Worked Example 3
Problem. Variant B beats A within both mobile and desktop, yet A wins overall. Explain this Simpson's paradox.
Answer. An uneven segment mix (a confounder) reverses the pooled result; stratify or reweight by segment to recover B's true within-group advantage.
Problem. An A/B test is monitored daily for 10 days and the team will stop at the first day with p < 0.05, with the null actually true. Give a rough upper bound on the false-positive probability treating the looks as independent, and name one proper fix.
Solution. Treating the 10 daily looks as independent, P(at least one false positive) = 1 - (1 - 0.05)^10 = 1 - 0.95^10 = 1 - 0.5987 = 0.4013, roughly a 40% false-positive rate instead of 5%. (Real sequential looks are positively correlated, so the true inflation is somewhat lower but still far above 5%.) A proper fix is to pre-commit to a fixed sample size and analyze once, or use a sequential testing method with alpha-spending boundaries (or always-valid/Bayesian sequential tests) that control the error rate across repeated looks.
Pick a conversion-rate scenario (e.g., baseline 5%, hoped-for lift to 6%). Use Python (statsmodels) or R to compute the power for a fixed sample size and the sample size needed for 80% power at alpha = 0.05. Then simulate the experiment: generate treatment and control outcomes, run the appropriate test, and report whether you would have detected the effect. Repeat the simulation many times to estimate empirical power and compare it to the formula.
Deliverable · A notebook showing the power and sample-size calculations, the single-run test result, and the empirical power from repeated simulations, with a note on how peeking would distort the conclusion.
1. Why does randomly assigning units to treatment and control support causal claims?
Answer B. Randomization balances both known and unknown confounders in expectation, so a measured difference can be causally attributed to the treatment.
2. Statistical power is:
Answer B. Power is the probability of detecting a true effect, equal to one minus the Type II error rate; it rises with sample size and effect size.
3. Holding the effect size and alpha fixed, increasing the sample size will generally:
Answer B. Larger samples shrink the standard error, making true effects easier to detect, which raises power without changing alpha.
4. Repeatedly checking an A/B test and stopping as soon as it looks significant (peeking) primarily causes:
Answer B. Each peek is another chance to cross the threshold by luck, so naive repeated testing inflates the overall Type I error rate well above alpha.
5. Simpson's paradox refers to a situation where:
Answer B. Simpson's paradox occurs when aggregation hides or reverses subgroup trends, usually because of an unbalanced confounding variable across groups.
I can design a randomized A/B test with a clear metric, hypothesis, and randomization unit.
I can compute the power and required sample size for a target effect, and analyze the result correctly.
I can identify common experiment pitfalls such as peeking, novelty effects, and Simpson's paradox.
Fit a straight line to two variables by minimizing squared vertical distances, yielding a slope and intercept.
Simple linear regression models a continuous response y as y = beta0 + beta1*x + error, a straight line plus random noise. Ordinary least squares (OLS) picks the intercept beta0 and slope beta1 that minimize the sum of squared residuals, sum (y_i - y-hat_i)^2. The closed-form solutions are beta1-hat = Sxy/Sxx = (sum (x_i - x-bar)(y_i - y-bar)) / (sum (x_i - x-bar)^2) and beta0-hat = y-bar - beta1-hat*x-bar, so the fitted line always passes through (x-bar, y-bar). The slope is the estimated change in y per one-unit increase in x. Squaring residuals penalizes large misses and gives a unique, differentiable solution. This is the foundational model of statistics and machine learning: interpretable, fast, and the base case for everything that follows.
Worked Example 1
Problem. Fit a line to x = [1,2,3,4,5], y = [2,4,5,4,5]. Find the slope and intercept.
Answer. y-hat = 2.2 + 0.6*x (slope 0.6, intercept 2.2).
Worked Example 2
Problem. Using y-hat = 2.2 + 0.6x, predict y at x = 4 and compute that residual.
Answer. Predicted y = 4.6; residual = 4 - 4.6 = -0.6.
Worked Example 3
Problem. Interpret a fitted housing model price-hat = 50 + 0.2*sqft (price in $thousands). What does the slope mean and what is a 1500 sqft prediction?
Answer. Each extra sqft adds ~$200 to predicted price; a 1,500 sqft home is predicted at $350,000.
Problem. Fit a line by least squares to x = [0,1,2,3], y = [1,3,2,5]. Find the slope, intercept, and predict y at x = 5.
Solution. x-bar = 1.5, y-bar = 2.75. x deviations: -1.5,-0.5,0.5,1.5; y deviations: -1.75,0.25,-0.75,2.25. Sxy = (-1.5)(-1.75)+(-0.5)(0.25)+(0.5)(-0.75)+(1.5)(2.25) = 2.625 - 0.125 - 0.375 + 3.375 = 5.5. Sxx = 2.25+0.25+0.25+2.25 = 5. beta1-hat = 5.5/5 = 1.1. beta0-hat = 2.75 - 1.1*1.5 = 2.75 - 1.65 = 1.1. So y-hat = 1.1 + 1.1*x; at x = 5, y-hat = 1.1 + 5.5 = 6.6. (np.polyfit confirms slope 1.1, intercept 1.1.)
With several predictors, each coefficient is a partial effect holding the others constant.
Multiple linear regression extends the model to several predictors: y = beta0 + beta1*x1 + ... + betap*xp + error. OLS still minimizes the sum of squared residuals, now solving for a vector of coefficients (in matrix form beta-hat = (X'X)^{-1} X'y). The key interpretive shift: each beta_j is a partial effect — the expected change in y for a one-unit increase in x_j while holding all other predictors fixed. This 'controlling for' interpretation is what lets regression disentangle overlapping influences, but it weakens when predictors are highly correlated (multicollinearity), which inflates coefficient variances and makes individual effects hard to separate. Categorical predictors enter via dummy (indicator) variables interpreted relative to a reference level. Multiple regression is the everyday tool for prediction and for adjusted, ceteris paribus association in data science.
Worked Example 1
Problem. A model is salary-hat = 30 + 2*years + 5*degree (degree is 1 if graduate, else 0; salary in $thousands). Interpret each coefficient.
Answer. Base $30k; +$2k per year (holding degree); +$5k for a graduate degree (holding years); predicted $45k for 5 years with a degree.
Worked Example 2
Problem. Why might a predictor's coefficient change sign when a second correlated predictor is added?
Answer. Because the simple-regression coefficient was biased by an omitted correlated variable; controlling for it reveals the true partial effect, which can change sign.
Worked Example 3
Problem. Sketch statsmodels code to fit and read coefficients for y on x1 and x2.
Answer. sm.add_constant then sm.OLS(...).fit(); model.params gives the intercept and each predictor's partial slope.
Problem. A model predicts monthly spend (dollars): spend-hat = 20 + 0.5*age + 15*member, where member = 1 for members, 0 otherwise. Interpret the 15 coefficient and predict spend for a 40-year-old member.
Solution. The coefficient 15 on member means that, holding age fixed, members are predicted to spend $15 more per month than non-members. Prediction for a 40-year-old member: spend-hat = 20 + 0.5*40 + 15*1 = 20 + 20 + 15 = 55, so $55 per month. (A 40-year-old non-member would be predicted at 20 + 20 + 0 = $40.)
Each coefficient comes with a standard error, t-statistic, p-value, and confidence interval for inference.
Regression estimates are statistics, so they carry uncertainty. Each coefficient beta_j-hat has a standard error SE(beta_j-hat) reflecting its sampling variability. To test H0: beta_j = 0 (the predictor has no linear effect), compute t = beta_j-hat / SE(beta_j-hat), compared to a t-distribution with n - p - 1 degrees of freedom (n observations, p predictors plus intercept). A small p-value indicates the predictor is significantly associated with y. A confidence interval is beta_j-hat +/- t* * SE(beta_j-hat). The overall F-test asks whether the model as a whole explains more than the intercept-only model. Valid inference rests on the regression assumptions (linearity, independent errors, constant variance, approximately normal errors); violations distort the standard errors. Reporting coefficients with their intervals — not just point estimates — is essential honest practice.
Worked Example 1
Problem. A coefficient is beta-hat = 0.80 with SE = 0.20, n = 50, p = 3 predictors. Test H0: beta = 0.
Answer. t = 4.0 on 46 df, p ~ 0.0002; the coefficient is significantly different from 0.
Worked Example 2
Problem. For that coefficient (beta-hat = 0.80, SE = 0.20, df = 46), give an approximate 95% confidence interval.
Answer. 95% CI ~ [0.397, 1.203]; it excludes 0, agreeing with the significant t-test.
Worked Example 3
Problem. What does the model's overall F-test assess, and how does it differ from the per-coefficient t-tests?
Answer. The F-test jointly tests whether all slopes are zero (overall model usefulness); t-tests assess each coefficient one at a time.
Problem. A regression coefficient is beta-hat = -1.5 with SE = 0.5 from a model with n = 24 and p = 2 predictors. Test H0: beta = 0 at alpha = 0.05 and give the t-statistic and df.
Solution. t = beta-hat/SE = -1.5/0.5 = -3.0. df = n - p - 1 = 24 - 2 - 1 = 21. Two-sided p = 2*P(T21 >= 3.0); scipy: 2*(1 - stats.t.cdf(3.0, 21)) = 0.0069. Since 0.0069 < 0.05 (and |t| = 3.0 exceeds the critical t ~ 2.080 for 21 df), reject H0 — the predictor has a significant negative association with the response.
Residual plots reveal whether the regression assumptions hold and flag influential points.
OLS inference relies on four assumptions, all checkable through residuals (observed minus fitted). Linearity: a residuals-vs-fitted plot should show no pattern; curvature signals a missing nonlinear term. Constant variance (homoscedasticity): the residual spread should be roughly uniform across fitted values; a funnel shape indicates heteroscedasticity. Normality of errors: a Q-Q plot of residuals should be roughly straight, important for small-sample inference. Independence: residuals should not be autocorrelated (relevant for time series). Beyond assumptions, watch for influential points using leverage (extreme x) and Cook's distance (impact on the fitted coefficients); a high-leverage outlier can single-handedly tilt the line. Diagnostics turn regression from a black box into a model you can trust or fix — transform variables, add terms, or use robust methods as the plots dictate.
Worked Example 1
Problem. A residuals-vs-fitted plot shows a clear U-shaped (parabolic) pattern. What is violated and what is a fix?
Answer. The linearity assumption is violated; add a polynomial term (e.g., x^2) or transform a variable to capture the curvature.
Worked Example 2
Problem. Residual spread is small for low fitted values and large for high ones (a funnel). Name the issue and a remedy.
Answer. Heteroscedasticity (non-constant variance); fix via a log/variance-stabilizing transform, weighted least squares, or robust standard errors.
Worked Example 3
Problem. One data point has very high leverage and a large Cook's distance. Why does it matter and how do you respond?
Answer. It is an influential point that can dominate the fit; investigate it, refit with and without it, and report how the conclusions change.
Problem. After fitting a regression you make a Q-Q plot of the residuals and it bends sharply up at both ends (heavy tails). Which assumption is in question, when does it matter most, and what is one response?
Solution. An S-shaped or both-ends-bending Q-Q plot indicates the residuals are not normally distributed (here, heavier tails than normal). The normality assumption matters most for inference (t-tests, confidence intervals) in small samples; with large n the CLT makes coefficient estimates approximately normal regardless. One response is to transform the response (e.g., log y) to tame the tails, add a missing predictor, or use robust/bootstrap standard errors that do not rely on normal errors. In Python: import statsmodels.api as sm; sm.qqplot(model.resid, line='45').
Model the log-odds of a yes/no outcome as a linear function of predictors, fit by maximum likelihood.
When the response is binary (0/1), linear regression is inappropriate because it can predict probabilities outside [0,1]. Logistic regression instead models the log-odds (logit) of the positive class as linear: log(p/(1 - p)) = beta0 + beta1*x1 + ... + betap*xp. Solving for p gives the S-shaped logistic (sigmoid) curve p = 1/(1 + e^{-(linear predictor)}), always between 0 and 1. Coefficients are interpreted on the odds scale: a one-unit increase in x_j multiplies the odds by e^{beta_j} (the odds ratio). So e^{beta_j} > 1 raises the odds, < 1 lowers them. There is no closed-form solution, so coefficients are estimated by maximum likelihood via iterative optimization. Logistic regression is the standard interpretable classifier in data science for problems like churn, default, and click prediction.
Worked Example 1
Problem. A logistic model has logit(p) = -2 + 0.5*x. Find the predicted probability at x = 4.
Answer. Predicted probability p = 0.5 at x = 4.
Worked Example 2
Problem. Interpret a coefficient beta = 0.7 on a predictor as an odds ratio.
Answer. Odds ratio = e^{0.7} ~ 2.01; each one-unit increase roughly doubles the odds of the positive outcome.
Worked Example 3
Problem. At x = 6 with logit(p) = -2 + 0.5*x, find the predicted probability and the odds.
Answer. p ~ 0.731 and odds ~ 2.72 (e^1) at x = 6.
Problem. A logistic model is logit(p) = -1 + 0.8*x. (a) Find the predicted probability at x = 1. (b) Interpret the coefficient 0.8 as an odds ratio.
Solution. (a) Linear predictor at x = 1 is -1 + 0.8*1 = -0.2. p = 1/(1 + e^{0.2}) = 1/(1 + 1.2214) = 1/2.2214 = 0.450. So the predicted probability is about 0.45. (b) The odds ratio is e^{0.8} = 2.2255, meaning each one-unit increase in x multiplies the odds of the positive outcome by about 2.23 (a 123% increase in odds), holding other predictors fixed. In Python: 1/(1+np.exp(0.2)) = 0.4502; np.exp(0.8) = 2.2255.
GLMs unify regression by pairing a response distribution with a link function connecting predictors to the mean.
A generalized linear model (GLM) generalizes linear regression through three components: a random component (the response's distribution from the exponential family, e.g., normal, Bernoulli, Poisson), a systematic component (the linear predictor eta = beta0 + sum beta_j*x_j), and a link function g connecting them via g(mean) = eta. The link maps the constrained mean onto the unbounded linear scale: the identity link recovers ordinary linear regression (normal response), the logit link gives logistic regression (Bernoulli response), and the log link gives Poisson regression for count data (ensuring positive means). This framework lets one fitting machinery (iteratively reweighted least squares maximizing the likelihood) handle continuous, binary, and count outcomes with consistent inference. Recognizing which distribution and link suit a problem is a core modeling skill that subsumes the earlier regressions as special cases.
Worked Example 1
Problem. Name the response distribution and link for (a) ordinary linear regression and (b) logistic regression.
Answer. (a) Normal distribution, identity link; (b) Bernoulli distribution, logit link.
Worked Example 2
Problem. For count data (e.g., number of website visits), which GLM fits, and why a log link?
Answer. Poisson regression with a log link; the log link keeps the predicted mean count positive.
Worked Example 3
Problem. In a Poisson (log-link) model, log(mu) = 0.5 + 0.3*x. Interpret the 0.3 coefficient and find mu at x = 2.
Answer. Each unit of x multiplies the expected count by e^{0.3} ~ 1.35; at x = 2, mu = e^{1.1} ~ 3.0.
Problem. You model the number of support tickets per customer (a count) with predictors. (a) State the appropriate GLM distribution and link. (b) If the fitted model is log(mu) = 0.2 + 0.5*plan, where plan is 1 for premium and 0 otherwise, interpret the plan coefficient and predict the mean tickets for a premium customer.
Solution. (a) Counts call for Poisson regression with a log link. (b) The coefficient 0.5 on plan means premium customers' expected ticket count is multiplied by e^{0.5} = 1.6487, about a 65% increase versus non-premium, holding other predictors fixed. Prediction for a premium customer: log(mu) = 0.2 + 0.5*1 = 0.7, so mu = e^{0.7} = 2.0138, about 2.0 tickets. (A non-premium customer: mu = e^{0.2} = 1.221.) In Python: np.exp(0.5) = 1.6487; np.exp(0.7) = 2.0138.
Using a dataset with a continuous outcome, fit a multiple linear regression in Python (statsmodels) or R, interpret the coefficients and R-squared, and run residual diagnostics to check the assumptions. Then take a binary outcome from the same domain and fit a logistic regression, interpreting one coefficient as an odds ratio. Compare what each model framework is doing.
Deliverable · A report with the linear model summary and diagnostic plots, the logistic model summary with an odds-ratio interpretation, and a paragraph contrasting the two models.
1. Ordinary least squares fits a regression line by:
Answer A. OLS chooses coefficients that minimize the sum of squared vertical distances (residuals) between observed and fitted values.
2. In a multiple regression, a coefficient represents:
Answer B. Each slope is a partial effect: the change in the response for a unit change in that predictor while the other predictors are held constant.
3. A residual plot showing a clear curved pattern most likely indicates:
Answer B. Systematic curvature in residuals signals that a linear form is inadequate; the relationship may be nonlinear or need a transformation.
4. Logistic regression models which quantity as a linear function of the predictors?
Answer B. Logistic regression uses the logit link, so the log-odds of the positive class is linear in the predictors; probabilities are recovered via the inverse logit.
5. A generalized linear model is defined by three components: a response distribution, a linear predictor, and a:
Answer B. A GLM combines a distribution from the exponential family, a linear predictor, and a link function relating the predictor to the mean of the response.
I can fit and interpret simple and multiple linear regression models and test their coefficients.
I can check regression assumptions with residual diagnostics and identify influential points.
I can fit a logistic regression and explain how GLMs extend linear models through link functions.
The two schools differ in what 'probability' means and whether a parameter can have a distribution.
Frequentist inference treats a parameter theta as a fixed unknown constant; probability is the long-run frequency of events over repeated sampling, so a confidence interval is a statement about the procedure, not the parameter. Bayesian inference treats theta as a random variable with a probability distribution that encodes our belief, updated by data via Bayes' rule. The Bayesian workflow is: choose a prior p(theta), observe data with likelihood p(data | theta), and obtain the posterior p(theta | data) proportional to likelihood times prior. The big payoff is interpretation: a Bayesian 95% credible interval lets you say directly 'there is a 95% probability theta lies here,' a statement a frequentist CI cannot make. For data scientists, the Bayesian view shines when you have prior knowledge, small samples, or want probability statements about hypotheses themselves.
Worked Example 1
Problem. A frequentist reports a 95% CI of [0.30, 0.38] for a conversion rate; a Bayesian reports a 95% credible interval of [0.30, 0.38]. How do their interpretations differ?
Answer. Numerically identical here, but the Bayesian interval makes a direct probability statement about p, while the frequentist interval describes the procedure's long-run coverage.
Worked Example 2
Problem. You flip a coin 3 times and get 3 heads. Contrast the frequentist MLE of p(heads) with a Bayesian estimate using a uniform prior.
Answer. The frequentist MLE is an overconfident 1.0; the Bayesian posterior mean of 0.8 (Beta(4,1)) regularizes the estimate using the prior, more sensible for tiny n.
Worked Example 3
Problem. Why can a Bayesian directly answer 'what is the probability the new variant beats the old one?' while a frequentist test cannot?
Answer. Because Bayesian parameters have distributions, P(p_B > p_A | data) is a direct, computable quantity; the frequentist framework treats parameters as fixed and only bounds error rates.
Problem. Explain, in one or two sentences each, (a) what probability means to a frequentist vs. a Bayesian, and (b) which framework lets you say 'there is a 90% probability the parameter is between 2 and 5.'
Solution. (a) To a frequentist, probability is the long-run relative frequency of an event over many repetitions, and a parameter is a fixed unknown constant; to a Bayesian, probability is a degree of belief, and the parameter itself has a probability distribution updated by data. (b) Only the Bayesian framework can make the direct statement 'there is a 90% probability the parameter is between 2 and 5,' because it treats the parameter as a random variable with a posterior distribution; the frequentist can only say 90% of such intervals would capture the parameter over repeated sampling.
Bayes' rule combines what you believed before (prior) with what the data say (likelihood) into an updated belief (posterior).
Bayesian inference runs on Bayes' rule: p(theta | data) = p(data | theta) * p(theta) / p(data). In words, posterior proportional to likelihood times prior, since the denominator p(data) (the evidence) is just a normalizing constant that does not depend on theta. The prior p(theta) encodes belief before data — it can be informative (a sharp guess) or weak/flat (let the data speak). The likelihood p(data | theta) is the same object from MLE: how probable the observed data are for each theta. Multiplying them and renormalizing yields the posterior, the full updated distribution over theta. The shape of the posterior is a compromise: with little data the prior dominates, and as data accumulate the likelihood takes over and the posterior concentrates around the data-driven value. This update can be done repeatedly, yesterday's posterior becoming today's prior.
Worked Example 1
Problem. State Bayes' rule for a parameter theta and explain why we usually write 'posterior proportional to likelihood times prior.'
Answer. p(theta | data) proportional to p(data | theta) * p(theta); the omitted denominator is a theta-independent normalizing constant fixed by making the posterior integrate to 1.
Worked Example 2
Problem. A disease has prevalence 1% (prior). A test is 99% sensitive and 95% specific. Given a positive test, find the posterior probability of disease.
Answer. P(disease | positive) ~ 0.167; despite a positive test, the low prior keeps the posterior near 17%, illustrating how strongly the prior matters with rare conditions.
Worked Example 3
Problem. Show how the same final posterior arises whether you update with two data batches sequentially or all at once.
Answer. Because independent likelihoods multiply, sequential updating (yesterday's posterior becomes today's prior) gives exactly the same posterior as updating on all data at once.
Problem. Spam filtering: 20% of emails are spam (prior). The word 'free' appears in 60% of spam emails and 5% of non-spam. An email contains 'free.' What is the posterior probability it is spam?
Solution. Prior: P(spam) = 0.20, P(ham) = 0.80. Likelihoods: P('free' | spam) = 0.60, P('free' | ham) = 0.05. Evidence P('free') = 0.60*0.20 + 0.05*0.80 = 0.12 + 0.04 = 0.16. Posterior P(spam | 'free') = (0.60*0.20)/0.16 = 0.12/0.16 = 0.75. In Python: (0.6*0.2)/(0.6*0.2 + 0.05*0.8) = 0.75. So there is a 75% posterior probability the email is spam.
When the prior and likelihood pair so the posterior stays in the prior's family, updating becomes simple arithmetic.
A conjugate prior is chosen so the posterior belongs to the same distributional family as the prior, making Bayesian updating analytic with no integration. The flagship example is the Beta-Binomial model for a proportion p. The Beta(a, b) distribution lives on [0, 1] and is the natural prior for a probability; a and b act like prior pseudo-counts of successes and failures. With binomial data of s successes in n trials, the likelihood is proportional to p^s (1 - p)^{n - s}, and multiplying by the Beta(a, b) prior proportional to p^{a-1}(1 - p)^{b-1} gives a posterior proportional to p^{a+s-1}(1 - p)^{b+n-s-1} — exactly Beta(a + s, b + n - s). So you just add successes to a and failures to b. Beta(1,1) is the uniform (flat) prior. Conjugacy makes the Beta-Binomial the textbook entry point to Bayesian estimation of rates.
Worked Example 1
Problem. With a uniform Beta(1,1) prior, you observe 8 successes in 10 trials. Find the posterior distribution.
Answer. Posterior = Beta(9, 3).
Worked Example 2
Problem. Compute the posterior mean for Beta(9, 3) and compare it to the MLE p-hat = 8/10.
Answer. Posterior mean = 0.75, shrunk from the MLE of 0.80 toward the prior mean 0.5 because of the prior's two pseudo-observations.
Worked Example 3
Problem. Use an informative prior Beta(20, 20) (belief centered at 0.5) with the same data (8 of 10). How does the posterior compare to the flat-prior case?
Answer. Posterior = Beta(28, 22), mean ~ 0.56 — far closer to the prior's 0.5 than the flat-prior result of 0.75, showing a strong prior resists a small dataset.
Problem. Starting from a Beta(2, 2) prior, you run an A/B variant and see 30 conversions out of 50 visitors. Give the posterior distribution and its mean.
Solution. Posterior = Beta(a + s, b + n - s) = Beta(2 + 30, 2 + (50 - 30)) = Beta(32, 22). Posterior mean = a/(a + b) = 32/(32 + 22) = 32/54 = 0.5926, about 0.593. In Python: from scipy.stats import beta; beta(32, 22).mean() -> 0.5926. The raw conversion rate 30/50 = 0.60 is pulled slightly toward the prior mean 0.5.
The posterior is a full distribution; we report point summaries and a credible interval that carries a direct probability.
A posterior distribution holds everything we know about theta, but we usually summarize it. Common point estimates are the posterior mean (minimizes squared-error loss), the posterior median (robust, minimizes absolute-error loss), and the posterior mode or MAP (the peak). To express uncertainty we report a credible interval: a (1 - alpha) credible interval is any region containing (1 - alpha) of the posterior probability. The equal-tailed interval uses the alpha/2 and 1 - alpha/2 posterior quantiles; the highest posterior density (HPD) interval is the shortest interval with the required mass and is preferable for skewed posteriors. Unlike a confidence interval, a 95% credible interval supports the direct statement 'P(theta in interval | data) = 0.95.' For the Beta posterior these summaries are one-line computations in SciPy or R.
Worked Example 1
Problem. For a Beta(9, 3) posterior, compute the posterior mean and an equal-tailed 95% credible interval.
Answer. Posterior mean = 0.75; 95% equal-tailed credible interval ~ [0.482, 0.940], i.e., P(0.482 <= p <= 0.940 | data) = 0.95.
Worked Example 2
Problem. State the direct interpretation of the credible interval [0.482, 0.940] and contrast it with a confidence interval.
Answer. The credible interval says P(p in [0.482, 0.940] | data) = 0.95 directly about p, whereas a confidence interval makes a long-run statement about the procedure, not this interval.
Worked Example 3
Problem. When does the highest posterior density (HPD) interval differ noticeably from the equal-tailed interval, and how do you compute it?
Answer. They diverge for skewed posteriors, where the HPD is the shortest 95% interval (closer to the mode); compute it from posterior samples with a tool like ArviZ's hdi.
Problem. You have a Beta(32, 22) posterior for a conversion rate. Report the posterior mean and outline how to get a 90% equal-tailed credible interval in Python.
Solution. Posterior mean = a/(a + b) = 32/(32 + 22) = 32/54 = 0.5926. For a 90% equal-tailed credible interval, take the 5th and 95th posterior percentiles: from scipy.stats import beta; lo, hi = beta(32, 22).ppf([0.05, 0.95]), which returns approximately [0.483, 0.700]. Interpretation: given the data and prior, there is a 90% posterior probability the true conversion rate lies between about 0.483 and 0.700.
When the posterior has no closed form, approximate it on a grid or by drawing samples from it.
Conjugate models are convenient but rare; most real posteriors lack a closed form. Two beginner-friendly approximations help. Grid approximation discretizes the parameter space into a fine grid, evaluates the unnormalized posterior (likelihood times prior) at each grid point, and normalizes by dividing by the sum — giving a discrete approximation you can use for means, intervals, and plots. It is exact in the limit of a fine grid but scales badly beyond one or two parameters (the curse of dimensionality). Sampling methods instead draw many values from the posterior; the collection of samples represents the distribution, and any summary becomes a simple average or quantile of the draws. For conjugate models you can sample directly; for general models, Markov chain Monte Carlo (MCMC) tools like PyMC or Stan draw correlated samples. Both turn integration into counting.
Worked Example 1
Problem. Outline a grid approximation for a Beta(1,1) prior with 7 successes in 10 trials, and how to get the posterior mean.
Answer. Evaluate p^7(1-p)^3 on a fine grid, normalize by its sum, and the weighted average sum(p*post) ~ 0.667 reproduces the analytic Beta(8,4) mean.
Worked Example 2
Problem. Estimate the posterior mean and a 95% credible interval of a Beta(9,3) posterior by drawing samples.
Answer. Sampling 100k Beta(9,3) draws gives mean ~ 0.75 and 95% interval ~ [0.48, 0.94], matching the analytic results to Monte Carlo precision.
Worked Example 3
Problem. Why does grid approximation become impractical for models with many parameters?
Answer. Because a grid needs g^d points, cost explodes exponentially with the number of parameters (curse of dimensionality); MCMC sampling is used instead for high-dimensional posteriors.
Problem. Using a flat Beta(1,1) prior and data of 3 successes in 4 trials, outline a grid approximation in Python and state which Beta posterior it should reproduce.
Solution. import numpy as np; p = np.linspace(0, 1, 1001); post = p**3 * (1 - p)**1; post = post / post.sum(). The posterior mean is np.sum(p * post) and a credible interval comes from the cumulative sum of post. Analytically the posterior is Beta(1 + 3, 1 + 1) = Beta(4, 2) with mean 4/6 = 0.667; the grid approximation's mean should match ~0.667, confirming the method.
Model each variant's rate with a Beta posterior and answer the business question directly: P(B beats A).
Bayesian A/B testing reframes the comparison in terms most stakeholders actually want. Model each variant's conversion rate with a Beta posterior: with a Beta(a, b) prior and s successes in n trials, variant A's posterior is Beta(a + s_A, b + n_A - s_A), and likewise for B. The key quantity is P(p_B > p_A | data), computed by drawing many samples from each posterior and counting the fraction where B's draw exceeds A's. You can also report the posterior of the difference p_B - p_A or the expected uplift, and the expected loss of each choice for a decision-theoretic stop rule. Unlike fixed-horizon frequentist tests, Bayesian methods make optional stopping less hazardous (no p-value to inflate), though priors and decision thresholds should still be set in advance. The output — 'a 97% probability B is better, with an expected uplift of 1.2 points' — is directly actionable.
Worked Example 1
Problem. Variant A: 90 conversions in 1000. Variant B: 110 in 1000. With flat Beta(1,1) priors, write the two posteriors.
Answer. p_A ~ Beta(91, 911) (mean ~ 0.091) and p_B ~ Beta(111, 891) (mean ~ 0.111).
Worked Example 2
Problem. Estimate P(p_B > p_A) by Monte Carlo for the posteriors above.
Answer. P(p_B > p_A | data) ~ 0.93 from 200k posterior draws — strong but not overwhelming evidence that B is better.
Worked Example 3
Problem. Beyond P(B > A), report the posterior of the uplift p_B - p_A and the expected loss of wrongly shipping B.
Answer. Posterior mean uplift ~ 0.020 with 95% credible interval ~ [-0.006, 0.046] (it includes 0), and a small expected loss; at a 0.95 threshold you would keep testing rather than ship B yet.
Problem. Variant A converts 40 of 500 and variant B converts 55 of 500. Using flat Beta(1,1) priors, give both posteriors and outline how to estimate P(p_B > p_A) and the expected uplift in Python.
Solution. Posteriors: A = Beta(1 + 40, 1 + 460) = Beta(41, 461) (mean ~ 0.0817); B = Beta(1 + 55, 1 + 445) = Beta(56, 446) (mean ~ 0.1116). In Python: import numpy as np; from scipy.stats import beta; a = beta(41, 461).rvs(200000); b = beta(56, 446).rvs(200000); print((b > a).mean()) gives P(p_B > p_A) ~ 0.97, and print((b - a).mean()) gives the expected uplift ~ 0.030 (about 3 percentage points). With ~97% probability B beats A and a ~3-point uplift, shipping B is well supported.
Suppose you observe a number of successes out of trials for a conversion rate. Choose a Beta prior (try a uniform Beta(1,1) and a more informed prior) and derive the Beta posterior analytically. In Python or R, plot the prior, likelihood, and posterior, report the posterior mean and a 95% credible interval, and then confirm your result with a grid approximation. Finally, compare your Bayesian conclusion to a frequentist confidence interval for the same data.
Deliverable · A notebook plotting prior/likelihood/posterior for both priors, the posterior summaries and credible interval, the grid-approximation check, and a short comparison to the frequentist interval.
1. Bayes' rule for inference states that the posterior is proportional to:
Answer B. Posterior is proportional to likelihood times prior; the normalizing constant (the evidence) makes it a proper distribution.
2. A conjugate prior is one that:
Answer B. Conjugacy means the posterior stays in the prior's family (e.g., a Beta prior with binomial data gives a Beta posterior), making updates analytic.
3. With a Beta(a, b) prior and observing s successes in n trials (binomial), the posterior is:
Answer A. The Beta-Binomial update adds successes to the first parameter and failures (n - s) to the second: Beta(a + s, b + n - s).
4. How does a 95% credible interval differ from a 95% confidence interval?
Answer B. Under the Bayesian view the parameter has a distribution, so a credible interval carries a direct probability statement, unlike a frequentist confidence interval.
5. As you collect more data, the posterior distribution typically:
Answer B. With more data the likelihood dominates, so the posterior tightens around the data-driven estimate and the prior's influence fades.
I can explain the Bayesian view of probability and combine a prior with a likelihood to obtain a posterior.
I can use a conjugate Beta-Binomial model and summarize the posterior with a credible interval.
I can approximate a posterior by grid or sampling and apply Bayesian reasoning to an A/B test.
Where this leads
Part of the open-source Crunch Academy tech-education path. Pair it with the data-science and ML tracks when you're ready.
All Math for Data Science courses Crunch Academy home