Comprehensive R Code Guide for Statistical
Methods (Exam-oriented)
Abstract
This document organizes R code patterns by technique rather than by assignment. It
prioritizes IRLS and Bootstrapping (from your screenshot) while covering Newton–Raphson,
EM, GLMs, power simulation and quick R references that frequently appear in exams.
1 Part 1: Iterative Optimization (IRLS & Newton–Raphson)
Context: Screenshot (Q2) and Assignment 2. Concept: Find estimators (MLE, median)
when closed forms don’t exist.
1.1 1.1 IRLS for Least Absolute Deviations (Screenshot problem)
P
Theory: The median minimizes i |yi − θ|. IRLS implements this by iterating weighted means
with weights
1
wi = ,
|yi − θold |
(or a regularized version wi = 1/ max(|yi − θold |, ε) to avoid division by zero).
Generic exam template (R):
# IRLS for Median / LAD
[Link](123)
n <- 201; sigma <- 2; theta_true <- 1
delta <- 2 * rbinom(n, size=1, prob=0.5) - 1
epsilon <- rexp(n, rate = 1/sigma) * delta
y <- theta_true + epsilon
theta_old <- 1
theta_new <- 2
tol <- 1e-10
max_iter <- 100
iter <- 0
while(abs(theta_new - theta_old) > tol * abs(theta_old) && iter < max_iter) {
theta_old <- theta_new
weights <- 1 / pmax(abs(y - theta_old), 1e-12)
theta_new <- sum(y * weights) / sum(weights)
iter <- iter + 1
}
print(paste("Estimated Theta:", theta_new))
print(paste("Sample Median:", median(y)))
1
Exam tip: If weights are 1/|r| the criterion is L1 (absolute deviations); if weights are constant
the criterion is L2 (squared errors).
1.2 1.2 Steepest Ascent / Newton–Raphson
Context: Assignment 2 (Problem 6). Use when score (gradient) and optionally information
are available.
Steepest-ascent template (R):
# Steepest Ascent (gradient ascent)
theta <- 0.5 # initial
alpha <- 0.01 # step size
tol <- 1e-6
repeat {
grad <- gradient(theta, data) # user-supplied
theta_new <- theta + alpha * grad
theta_new <- min(max(theta_new, 0.01), 0.99) # enforce bounds if needed
if (abs(theta_new - theta) < tol) break
theta <- theta_new
}
Newton–Raphson variant: replace update by
U (θ)
θnew = θold −
I(θ)
where U is score and I is observed/expected information.
2 Part 2: The EM Algorithm (Expectation–Maximization)
Context: Assignment 2 (Poisson missing data, Gaussian mixtures) and Assignment 3.
2.1 2.1 EM for Gaussian Mixture Models
Model: data yi ∼ p N (µ1 , σ 2 ) + (1 − p) N (µ2 , σ 2 ). Goal: estimate p, µ1 , µ2 , σ 2 .
EM template (R):
EM_Algorithm <- function(y, p_init, mu1_init, mu2_init, s2_init, tol=1e-6) {
n <- length(y)
p <- p_init; mu1 <- mu1_init; mu2 <- mu2_init; s2 <- s2_init
converged <- FALSE; iter <- 0
while(!converged) {
p_old <- p; mu1_old <- mu1; mu2_old <- mu2; s2_old <- s2
d1 <- dnorm(y, mean = mu1, sd = sqrt(s2))
d2 <- dnorm(y, mean = mu2, sd = sqrt(s2))
w1 <- (p * d1) / (p * d1 + (1 - p) * d2)
w2 <- 1 - w1
p <- mean(w1)
mu1 <- sum(w1 * y) / sum(w1)
mu2 <- sum(w2 * y) / sum(w2)
2
rss1 <- sum(w1 * (y - mu1)^2)
rss2 <- sum(w2 * (y - mu2)^2)
s2 <- (rss1 + rss2) / n
param_change <- max(abs(c(p - p_old, mu1 - mu1_old, mu2 - mu2_old)))
if (param_change < tol) converged <- TRUE
iter <- iter + 1
}
list(p = p, mu1 = mu1, mu2 = mu2, sigma2 = s2, iterations = iter)
}
Notes:
• Initialize sensibly (k-means or random restarts).
• Monitor log-likelihood for monotone increase.
• Handle label-switching by ordering components or fixing constraints.
3 Part 3: Bootstrapping (Resampling)
Context: Assignment 5 and screenshot Q3. Estimate SE or bias without analytic formula.
3.1 3.1 Nonparametric Bootstrap
Template (R) for estimating SE of median:
original_data <- c(2.1, 2.9, 4.0, 4.8, 6.2)
n <- length(original_data)
B <- 1000
boot_estimates <- numeric(B)
[Link](123)
for (i in 1:B) {
idx <- sample(1:n, size = n, replace = TRUE)
boot_sample <- original_data[idx]
boot_estimates[i] <- median(boot_sample)
}
boot_se <- sd(boot_estimates)
boot_bias <- mean(boot_estimates) - median(original_data)
cat("Bootstrap SE:", boot_se, "\n")
3.2 3.2 Parametric (Model-based) Bootstrap
Template (R) for regression residual bootstrap:
fit <- lm(y ~ x)
fitted_vals <- fitted(fit)
resids <- residuals(fit) - mean(residuals(fit))
B <- 1000
boot_coefs <- matrix(NA, nrow = B, ncol = 2)
for (i in 1:B) {
boot_errors <- sample(resids, size = length(resids), replace = TRUE)
3
y_boot <- fitted_vals + boot_errors
fit_boot <- lm(y_boot ~ x)
boot_coefs[i, ] <- coef(fit_boot)
}
se_beta <- sd(boot_coefs[, 2])
4 Part 4: Regression (GLM)
Context: Assignment 4 (Challenger data) — logistic and probit.
GLM template:
# Logistic
model_logit <- glm(y ~ x1 + x2, family = binomial(link = "logit"))
summary(model_logit)
predict(model_logit, newdata = [Link](x1=30, x2=200), type = "response")
# Probit
model_probit <- glm(y ~ x1 + x2, family = binomial(link = "probit"))
Interpretation: type="response" returns probabilities; type="link" returns linear predictor
(log-odds for logit).
5 Part 5: Simulation (Hypothesis Testing & Power)
Context: Assignments 6 & 9. Simulate many datasets under H0 or H1 to estimate size/power.
Power/size template (R):
nsim <- 5000; n <- 100
lambda_0 <- 3; lambda_true <- 3.5; alpha <- 0.05
z_crit <- qnorm(1 - alpha)
gamma_crit <- qgamma(1 - alpha, shape = n, rate = lambda_0)
reject_large <- reject_exact <- 0
for (i in 1:nsim) {
x <- rexp(n, rate = lambda_true)
T_obs <- sum(x)
x_bar <- mean(x)
S <- sqrt(n) * (1 - lambda_0 * x_bar) # example score-like stat
if (S > z_crit) reject_large <- reject_large + 1
if (T_obs < gamma_crit) reject_exact <- reject_exact + 1
}
cat("Large-sample power:", reject_large / nsim, "\n")
cat("Exact power:", reject_exact / nsim, "\n")
4
6 Quick Reference: Common R Functions
Task Function / Syntax
Normal dist dnorm(), pnorm(), qnorm(), rnorm()
Exponential dexp(), pexp(), qexp(), rexp()
Binomial dbinom(), pbinom(), qbinom(), rbinom()
Gamma dgamma(), pgamma(), qgamma(), rgamma()
Optimization optim(par, fn); manual loops for NR/IRLS
Linear model lm(y x)
GLM (logistic/probit) glm(y x, family=binomial(link="logit"))
Matrix mult %*%
Summaries sum(), mean(), median(), var(), sd()
7 How to handle “similar but not exact” exam questions
1. Identify the distributional form implied by the code (e.g., rexp(...)*delta ⇒
Laplace/double exponential).
2. Identify the objective (what is updated in the loop: θ ⇒ estimation).
3. Identify the method (weights indicate IRLS; gradient indicates steepest-ascent or NR).
4. Map weights to loss: w = 1/|r| ⇒ L1 ; w = const ⇒ L2 .
5. Adapt code line-by-line (change density/function calls or weight formula to match new
model).
Final Notes
• Use sensible initialisation and monitor convergence (parameter change or log-likelihood).
• Prefer multiple random starts for EM to avoid local maxima.
• Regularize IRLS denominators to avoid numerical blow-ups.
• Translate any exam variant to these templates — changing distributional functions or
weight definitions usually suffices.