B Worked Exercise Solutions
This appendix provides complete, validated code implementations and analytical answers for all exercise questions across Chapters 1 to 5.
B.1 Chapter 1 Solutions
B.1.1 Subsection 1.1: Foundations of Research
B.1.1.1 Exercise 1: Research Methods vs. Research Methodology
- Research Methods: Refer to the actual techniques and instruments used to perform research operations (e.g., questionnaires, interview checklists, statistical tests).
- Research Methodology: Refers to the systematic, theoretical analysis of the methods applied to a field of study. It is the science of understanding how research is done scientifically and why specific methods are selected.
B.1.1.2 Exercise 2: Objectives of Exploratory Research
Exploratory research aims to gain background information, define terms, clarify problems, formulate hypotheses, and establish research priorities in areas where little prior study exists.
- Example: A researcher conducting interviews and focus groups to explore how users interact with a brand-new generative AI tool to identify user patterns before designing a large-scale quantitative survey.
B.1.2 Subsection 1.2: Research Design & Formats
B.1.2.1 Exercise 3: Why Double-Blind RCT is the Gold Standard
A double-blind Randomized Controlled Trial (RCT) is the gold standard for causal inference because:
- Randomization: Distributes both known and unknown confounding variables equally across treatment groups.
- Double-Blinding: Prevents observer bias and subject placebo effects, ensuring that differences in outcomes are solely due to the intervention.
B.1.2.2 Exercise 4: Research Protocol Outline for Student Stress
- Title Page: Study title, principal investigator, institution, version number, and date.
- Project Summary: Rationale, objectives, and summary of the method.
- Background/Literature Review: Gap in current literature regarding university stress.
- Study Objectives: Primary hypothesis (e.g., mindfulness program reduces salivary cortisol).
- Methodology:
- Study design (2-arm parallel RCT).
- Selection criteria (undergraduate students).
- Interventions and measurements (cortisol assays, PSS questionnaire).
- Data Analysis Plan: Student’s t-test and ANOVA.
- Ethical Considerations: Informed consent, data security, and counseling resources.
- Budget & Timeline: Timeline chart (Gantt chart) and itemized budget.
B.1.3 Subsection 1.3: Sampling & Ethics
B.1.3.1 Exercise 5: Probability vs. Non-Probability Sampling Designs
- Probability Sampling: Every member of the population has a known, non-zero chance of being selected (e.g., Simple Random Sampling, Stratified Sampling). It allows for unbiased generalization and statistical inference.
- Non-Probability Sampling: Selection is based on non-random criteria (e.g., Convenience Sampling, Snowball Sampling). It is faster and cheaper but susceptible to selection bias and cannot be used to estimate sampling errors.
B.1.3.2 Exercise 6: Importance of Institutional Review Boards (IRBs)
Institutional Review Boards (IRBs) review and monitor research involving human subjects to protect their rights and welfare. They ensure compliance with ethical principles (autonomy, beneficence, justice) from the Belmont Report, verifying that risks are minimized, informed consent is obtained, and confidentiality is maintained.
B.2 Chapter 2 Solutions
B.2.1 Subsection 2.1: Random Number & Variable Generation
B.2.1.1 Exercise 1: Inverse Transform for Cauchy Distribution
The Cauchy distribution has CDF \(F(x) = \frac{1}{\pi} \arctan\left(\frac{x - x_0}{\gamma}\right) + \frac{1}{2}\). Setting \(U = F(X)\) and solving for \(X\): \[U - \frac{1}{2} = \frac{1}{\pi} \arctan\left(\frac{X - x_0}{\gamma}\right)\] \[\arctan\left(\frac{X - x_0}{\gamma}\right) = \pi \left(U - \frac{1}{2}\right)\] \[\frac{X - x_0}{\gamma} = \tan\left(\pi \left(U - \frac{1}{2}\right)\right)\] \[X = x_0 + \gamma \tan\left(\pi \left(U - \frac{1}{2}\right)\right) \quad \blacksquare\]
B.2.1.2 Exercise 2: Accept-Reject for Beta(2, 2)
The target Beta(2, 2) density is \(f(x) = 6x(1-x)\) for \(x \in [0, 1]\). We use Uniform(0, 1) proposal \(g(x) = 1\). The bounding constant \(c\) is: \[c = \max_x \frac{f(x)}{g(x)} = \max_x 6x(1-x)\] Differentiating \(6x - 6x^2\) gives \(6 - 12x = 0 \implies x = 0.5\). So the maximum is \(c = 6(0.5)(0.5) = 1.5\). The Accept-Reject algorithm is:
- Draw \(Y \sim \text{Unif}(0, 1)\) and \(U \sim \text{Unif}(0, 1)\).
- Accept \(Y\) if \(U \le \frac{6Y(1-Y)}{1.5 \times 1} = 4Y(1-Y)\), else repeat.
B.2.2 Subsection 2.2: Random Vectors & Permutations
B.2.2.1 Exercise 3: Why Cholesky Factor is Used to Introduce Covariance
If \(\mathbf{Z} \sim N(\mathbf{0}, \mathbf{I})\) is a vector of independent standard normals, its covariance matrix is \(\mathbf{I}\). Let \(\mathbf{X} = \mathbf{L}\mathbf{Z}\) where \(\boldsymbol{\Sigma} = \mathbf{L}\mathbf{L}^T\). The covariance of \(\mathbf{X}\) is: \[\text{Cov}(\mathbf{X}) = \text{Cov}(\mathbf{L}\mathbf{Z}) = \mathbf{L} \text{Cov}(\mathbf{Z}) \mathbf{L}^T = \mathbf{L} \mathbf{I} \mathbf{L}^T = \mathbf{L}\mathbf{L}^T = \boldsymbol{\Sigma}\] Thus, multiplying by the Cholesky factor \(\mathbf{L}\) maps independent components to have the exact target covariance structure.
B.2.2.2 Exercise 4: Fisher-Yates Permutations in R with Uniformity Check
fisher_yates_permute <- function(x) {
n <- length(x)
for (i in n:2) {
j <- sample(1:i, 1)
# Swap
temp <- x[i]
x[i] <- x[j]
x[j] <- temp
}
return(x)
}
# Uniformity check: simulate permutations of c(1, 2, 3)
perms <- replicate(10000, paste(fisher_yates_permute(c(1, 2, 3)), collapse=""))
table(perms) / 10000 # Each of the 6 combinations should be approx 1/6 = 0.1667## perms
## 123 132 213 231 312 321
## 0.1643 0.1721 0.1627 0.1706 0.1662 0.1641B.2.3 Subsection 2.3: Processes & Influence Functions
B.2.3.1 Exercise 5: Non-Homogeneous Poisson Process R Script
We simulate via thinning: we generate arrivals from a homogeneous process with rate \(\lambda_{\text{max}} = \max_{t \in [0, T]} \lambda(t)\) and accept them with probability \(\lambda(t) / \lambda_{\text{max}}\):
simulate_nhpp <- function(T_end = 5) {
# max rate for lambda(t) = 3 + 2t on [0, T_end] is at t = T_end
lambda_max <- 3 + 2 * T_end
arrivals <- numeric(0)
t <- 0
while (t < T_end) {
u1 <- runif(1)
t <- t - log(u1) / lambda_max
if (t >= T_end) break
u2 <- runif(1)
if (u2 <= (3 + 2 * t) / lambda_max) {
arrivals <- c(arrivals, t)
}
}
return(arrivals)
}
simulate_nhpp(5)## [1] 1.031922 1.293753 1.525183 1.572654 1.697949 2.074189 2.132589 2.153012
## [9] 2.266981 2.592818 2.599207 2.757676 2.804549 2.900881 3.113941 3.215905
## [17] 3.265490 3.521151 3.572516 3.678333 3.735882 3.784005 3.832224 3.838813
## [25] 3.967317 3.970884 3.983267 3.996382 4.231266 4.268463 4.314231 4.625139
## [33] 4.642964 4.727450 4.770743 4.843408B.3 Chapter 3 Solutions
B.3.1 Subsection 3.1: Variance Reduction Techniques
B.3.1.1 Exercise 1: Optimal Control Coefficient for Multiple Control Variates
Let \(\mathbf{Y} = [Y_1, Y_2]^T\) be a vector of control variates with known expectations. The control estimator is: \[X_c = X + \mathbf{c}^T(\mathbf{Y} - E[\mathbf{Y}])\] The variance-minimizing vector of coefficients is: \[\mathbf{c}^* = -\boldsymbol{\Sigma}_{YY}^{-1} \boldsymbol{\Sigma}_{YX}\] where \(\boldsymbol{\Sigma}_{YY}\) is the covariance matrix of the control variates, and \(\boldsymbol{\Sigma}_{YX}\) is the vector of covariances between the target estimator \(X\) and the control variates \(\mathbf{Y}\).
B.3.1.2 Exercise 2: Stratified sampling program in R (Neyman vs Proportional)
neyman_allocation <- function(N, W, sigmas) {
return(round(N * (W * sigmas) / sum(W * sigmas)))
}
proportional_allocation <- function(N, W) {
return(round(N * W))
}
# Example: N=1000, W=c(0.3, 0.7), sigmas=c(1, 4)
neyman_allocation(1000, c(0.3, 0.7), c(1, 4))## [1] 97 903## [1] 300 700B.3.2 Subsection 3.2: Importance Sampling & HMMs
B.3.2.1 Exercise 3: Optimal Proposal Density in Importance Sampling
The optimal proposal minimizes the variance. The proof is detailed in Appendix A: the optimal density \(g^*(x)\) is proportional to \(|h(x)|f(x)\).
B.3.2.2 Exercise 4: Bootstrap Filter in R for Non-linear Transitions
Let state transition be \(X_t = 0.5 X_{t-1} + \sin(X_{t-1}) + V_t\) and measurement be \(Y_t = X_t^2 + W_t\) where \(V_t, W_t \sim N(0, 1)\).
run_nonlinear_particle_filter <- function(y, N = 500) {
T_steps <- length(y)
particles <- matrix(0, nrow = N, ncol = T_steps)
weights <- matrix(0, nrow = N, ncol = T_steps)
particles[, 1] <- rnorm(N, mean = 0, sd = 1)
weights[, 1] <- dnorm(y[1], mean = particles[, 1]^2, sd = 1)
weights[, 1] <- weights[, 1] / sum(weights[, 1])
for (t in 2:T_steps) {
resample_idx <- sample(1:N, size = N, replace = TRUE, prob = weights[, t - 1])
particles_resampled <- particles[resample_idx, t - 1]
# Non-linear state transition
particles[, t] <- 0.5 * particles_resampled + sin(particles_resampled) + rnorm(N, mean = 0, sd = 1)
# Non-linear measurement update
weights[, t] <- dnorm(y[t], mean = particles[, t]^2, sd = 1)
weights[, t] <- weights[, t] / sum(weights[, t])
}
return(colSums(particles * weights))
}B.3.3 Subsection 3.3: Resampling & Likelihoods
B.3.3.1 Exercise 5: Double Bootstrap Standard Error Bias Estimation
Double bootstrap runs a second bootstrap nested inside each primary bootstrap sample to estimate the bias or variability of the bootstrap statistic:
double_bootstrap_se_bias <- function(data, R1 = 200, R2 = 100) {
n <- length(data)
boot1_se <- numeric(R1)
for (i in 1:R1) {
sample1 <- sample(data, size = n, replace = TRUE)
# Nested bootstrap
boot2_stats <- numeric(R2)
for (j in 1:R2) {
sample2 <- sample(sample1, size = n, replace = TRUE)
boot2_stats[j] <- mean(sample2)
}
boot1_se[i] <- sd(boot2_stats)
}
return(mean(boot1_se) - sd(data)/sqrt(n))
}B.3.3.2 Exercise 6: Empirical Likelihood CI vs Student’s t-interval for Median
Empirical Likelihood (EL) makes no distributional assumptions, which makes it highly robust to skewness. In contrast, the Student’s t-interval assumes normality of the sample mean. For skewed distributions (e.g., Log-normal), the EL confidence interval will be asymmetric and reflect the true shape of the distribution, whereas the Student’s t-interval will remain symmetric, leading to lower coverage probability on the tail side.
B.3.3.3 Exercise 7: Cross-Validation
- R Code Implementation: We write a function that performs \(k\)-fold cross-validation on linear models in R:
run_kfold_cv <- function(data, k = 5) {
n <- nrow(data)
folds <- cut(seq(1, n), breaks = k, labels = FALSE)
folds <- sample(folds) # Randomize fold assignments
mse <- numeric(k)
for (i in 1:k) {
test_idx <- which(folds == i)
train <- data[-test_idx, ]
test <- data[test_idx, ]
fit <- lm(y ~ x, data = train)
preds <- predict(fit, newdata = test)
mse[i] <- mean((test$y - preds)^2)
}
return(mean(mse))
}B.4 Chapter 4 Solutions
B.4.1 Detailed Balance & Stationarity
B.4.1.1 Exercise 1: Detailed Balance is Sufficient but Not Necessary
Let’s construct a 3-state transition matrix that is stationary but not in detailed balance. Consider the cyclic transition matrix: \[\mathbf{P} = \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}\] The stationary distribution is \(\boldsymbol{\pi} = [1/3, 1/3, 1/3]^T\) since: \[\boldsymbol{\pi}^T \mathbf{P} = [1/3, 1/3, 1/3] \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix} = [1/3, 1/3, 1/3] = \boldsymbol{\pi}^T\] However, detailed balance is violated since: \[\pi(1) P(1, 2) = \frac{1}{3} \times 1 = \frac{1}{3}\] \[\pi(2) P(2, 1) = \frac{1}{3} \times 0 = 0\] Since \(\frac{1}{3} \ne 0\), detailed balance does not hold. \(\blacksquare\)
B.4.1.2 Exercise 2: Microscopic Reversibility and Detailed Balance
Detailed balance represents the concept of microscopic reversibility at equilibrium. In physical chemistry, it asserts that any molecular process and its reverse occur at equal average rates. In Markov chain theory, it guarantees that the probability flux from \(x\) to \(y\) matches the reverse flux from \(y\) to \(x\) under the stationary distribution \(\pi\).
B.4.2 Metropolis-Hastings
B.4.2.1 Exercise 1: Metropolis-Hastings Detailed Balance Proof
The continuous case algebraic proof is detailed in Appendix A.
B.4.2.2 Exercise 2: Show Symmetric Proposal Reduces to Metropolis
If proposal is symmetric: \(q(x, y) = q(y, x)\). Hastings’ acceptance ratio is: \[\alpha(x, y) = \min\left(1, \frac{\pi(y) q(y, x)}{\pi(x) q(x, y)}\right)\] Since \(q(x, y) = q(y, x)\), they cancel out: \[\alpha(x, y) = \min\left(1, \frac{\pi(y)}{\pi(x)}\right)\] which is exactly the original Metropolis acceptance probability. \(\blacksquare\)
B.4.3 Gibbs Sampling
B.4.3.1 Exercise 1: Gibbs Sampler for Simple Linear Regression
For regression model \(Y_i \sim N(\beta_0 + \beta_1 X_i, \sigma^2)\), with prior \(p(\beta) \sim N(\mu_0, \boldsymbol{\Sigma}_0)\), the full conditionals are conjugate normals:
Update \(\boldsymbol{\beta} \mid \sigma^2, \mathbf{y} \sim N(\mathbf{m}_n, \mathbf{V}_n)\) where: \[\mathbf{V}_n = \left(\boldsymbol{\Sigma}_0^{-1} + \frac{1}{\sigma^2}\mathbf{X}^T\mathbf{X}\right)^{-1}, \quad \mathbf{m}_n = \mathbf{V}_n \left(\boldsymbol{\Sigma}_0^{-1}\mu_0 + \frac{1}{\sigma^2}\mathbf{X}^T\mathbf{y}\right)\]
Update \(\sigma^2 \mid \boldsymbol{\beta}, \mathbf{y} \sim \text{Inv-Gamma}(a_n, b_n)\) where: \[a_n = a_0 + n/2, \quad b_n = b_0 + \frac{1}{2}\sum (y_i - \mathbf{x}_i \boldsymbol{\beta})^2\]
B.4.3.2 Exercise 2: Bivariate Normal Full Conditionals
For \(X = (X_1, X_2) \sim N(\boldsymbol{\mu}, \boldsymbol{\Sigma})\): \[X_1 \mid X_2 = x_2 \sim N\left( \mu_1 + \rho \frac{\sigma_1}{\sigma_2}(x_2 - \mu_2), \sigma_1^2(1 - \rho^2) \right)\] \[X_2 \mid X_1 = x_1 \sim N\left( \mu_2 + \rho \frac{\sigma_2}{\sigma_1}(x_1 - \mu_1), \sigma_2^2(1 - \rho^2) \right)\]
B.4.4 Ising and Potts Models
B.4.4.1 Exercise 1: Potts \(q=2\) is Equivalent to Ising Model
For \(q=2\), Potts Hamiltonian states \(\sigma_i \in \{1, 2\}\), and we can map these to Ising spins \(s_i \in \{+1, -1\}\) using transformation \(s_i = 2\sigma_i - 3\). The Potts interaction term \(\delta(\sigma_i, \sigma_j)\) is 1 if spins are equal, 0 if different. The Ising term \(s_i s_j\) is 1 if spins are equal, -1 if different. Thus: \[\delta(\sigma_i, \sigma_j) = \frac{s_i s_j + 1}{2}\] Substituting this into the Potts Hamiltonian yields the Ising Hamiltonian scaled by a constant factor \(1/2\) and shifted by a constant energy offset.
B.4.4.2 Exercise 2: Boundary Conditions and Spatial Coalescence
- Fixed Boundary Conditions: Force the boundaries to stay at a specific spin (e.g., all +1). This provides external alignment that speeds up coalescence to that specific phase but suppresses natural fluctuations.
- Periodic Boundary Conditions: Wrap the grid edges (connecting left-to-right, top-to-bottom). This preserves translation invariance, allows realistic cluster growth, but can slow down coalescence because clusters have no fixed boundaries to seed alignment.
B.4.5 Perfect Sampling (CFTP)
B.4.5.1 Exercise 1: Proof that Monotonicity Guarantees Coalescence
The induction proof showing that the minimum and maximum chains sandwich all intermediate paths is detailed in Appendix A.
B.4.5.2 Exercise 2: Forward-Coupling Bias in CFTP
If we start chains at time 0 and run forward to time \(T\), checking for coalescence at time \(T\), the coalesced state will be biased toward states that have shorter paths to coalescence. This is because the coalescence time \(T_{\text{coalesce}}\) is dependent on the random increments. In contrast, starting at \(-T\) in the past and running to 0 guarantees that the chain has reached the stationary distribution at time 0 regardless of the specific starting time.
B.4.6 Simulated Annealing
B.4.6.1 Exercise 1: Balancing Exploration and Exploitation
- High Temperature \(T\) (Exploration): When \(T\) is high, the acceptance probability \(\alpha = \min(1, e^{-\Delta E / T})\) for uphill moves (\(\Delta E > 0\)) is close to 1. This allows the chain to make random walks across energy barriers and explore the entire parameter space, preventing early trapping in local minima.
- Low Temperature \(T\) (Exploitation): As \(T \to 0\), the probability of accepting uphill moves decays to 0. The algorithm behaves like a greedy gradient descent (local exploitation), only accepting moves that decrease the objective function, locking the state into the global minimum.
B.4.6.2 Exercise 2: Simulated Annealing in R for Rastrigin Function
- R Code Implementation:
rastrigin_1d <- function(x) x^2 + 10 * (1 - cos(2 * pi * x))
optimize_rastrigin <- function(init = 2.0, t0 = 10, gamma = 0.9, max_iter = 200) {
curr_x <- init
curr_e <- rastrigin_1d(curr_x)
t <- t0
for (k in 1:max_iter) {
prop_x <- rnorm(1, mean = curr_x, sd = 0.2)
prop_e <- rastrigin_1d(prop_x)
delta_e <- prop_e - curr_e
alpha <- exp(-delta_e / t)
if (delta_e < 0 || runif(1) < alpha) {
curr_x <- prop_x
curr_e <- prop_e
}
t <- t * gamma
}
return(curr_x)
}B.5 Chapter 5 Solutions
B.5.1 Sensitivity Analysis & Score Function
B.5.1.1 Exercise 1: Score Function for Gamma Distribution Scale Parameter
Let \(f(x; \theta) = \frac{1}{\Gamma(k)\theta^k} x^{k-1} e^{-x/\theta}\). Taking the natural logarithm: \[\log f(x; \theta) = -\log \Gamma(k) - k \log \theta + (k - 1) \log x - \frac{x}{\theta}\] Differentiating with respect to \(\theta\): \[\frac{\partial}{\partial \theta} \log f(x; \theta) = -\frac{k}{\theta} + \frac{x}{\theta^2} = \frac{x - k\theta}{\theta^2}\] Thus the score function is \(S(X; \theta) = \frac{X - k\theta}{\theta^2}\).
B.5.1.2 Exercise 2: Expected Value of the Score Function
The proof that \(E_\theta[S(X; \theta)] = 0\) is detailed in Appendix A.
B.5.1.3 Exercise 3: Manual Walkthrough Sensitivity of \(E[X^2]\) for \(X \sim \text{Exp}(\lambda)\)
We want to estimate the derivative of \(E_\lambda[X^2]\) at \(\lambda = 2.0\) using uniform draw \(u = 0.5\). The inverse transform is \(X = -\frac{1}{\lambda} \ln(U)\).
Calculate Sample: \[x = -0.5 \ln(0.5) \approx 0.3466\]
Find Score Function: The density is \(f(x) = \lambda e^{-\lambda x}\). \[\log f(x) = \log \lambda - \lambda x \implies \frac{\partial}{\partial \lambda} \log f(x) = \frac{1}{\lambda} - x\]
Evaluate Single Sample Estimate: \[\text{Gradient estimate} = x^2 \left( \frac{1}{\lambda} - x \right) = (0.3466)^2 (0.5 - 0.3466) = 0.1201 \times 0.1534 \approx 0.0184\]
B.5.2 Robbins-Monro Stochastic Approximation
B.5.2.1 Exercise 1: Robbins-Monro R Script for \(\theta^3 - 8 = 0\)
We find the root of \(M(\theta) = \theta^3 - 8 = 0\) (analytical root: 2) with normal noise. We use a scaled step size \(a_n = 0.1 / n\) to prevent the cubic gradient from causing numerical overflow:
robbins_monro_solve_cube <- function(n_iter = 1000, init = 0.0) {
theta <- init
for (n in 1:n_iter) {
y_obs <- theta^3 - 8 + rnorm(1)
# Using c = 0.1 to maintain numerical stability for cubic gradient
a_n <- 0.1 / n
theta <- theta - a_n * y_obs
}
return(theta)
}
# Test root finding
robbins_monro_solve_cube()## [1] 1.996676B.5.2.2 Exercise 2: Why Step Size \(a_n = 1/n^2\) Fails Convergence
For \(a_n = 1/n^2\), the second condition \(\sum a_n^2 < \infty\) is satisfied. However, the first condition requires \(\sum a_n = \infty\). Since \(\sum_{n=1}^\infty \frac{1}{n^2} = \frac{\pi^2}{6} < \infty\), the total step size budget is bounded. If the initial parameter \(\theta_1\) is far from the true root \(\theta^*\), the algorithm will run out of step capacity and freeze before reaching the root, causing convergence to fail.
B.5.2.3 Exercise 3: Compare Robbins-Monro step sizes \(a_n = 1/n\) vs \(a_n = 1/\sqrt{n}\)
- \(a_n = 1/n\): Converges asymptotically to the root and satisfies both conditions of the Robbins-Monro theorem. It suppresses the noise variance effectively at later steps.
- \(a_n = 1/\sqrt{n}\): Fails the second condition because \(\sum a_n^2 = \sum \frac{1}{n} = \infty\). This means the steps do not decrease fast enough to suppress noise, and the estimator will continue to oscillate wildly around the root instead of converging.
B.5.3 Rare-Event Simulation
B.5.3.1 Exercise 1: Fixed Factor vs. Fixed Effort Splitting
- Fixed Factor Splitting: Each trajectory that crosses a boundary is split into a pre-defined factor \(r_k\) of independent paths. The total number of simulated trajectories is random, which can lead to computational “explosions” or “extinction” (no paths crossing).
- Fixed Effort Splitting: Dictates simulating a fixed number of paths \(N\) at each level, resampling starting points only from successful trajectories of the previous level. This controls computational budget exactly and prevents explosions or extinctions.
B.5.3.2 Exercise 2: Bounded Relative Error in Importance Sampling for \(P(Z \ge 5)\)
Let \(Z \sim N(0, 1)\) and target probability be \(p = P(Z \ge 5)\). We select a shifted proposal \(g(x) \sim N(5, 1)\). The likelihood ratio is: \[w(x) = \frac{f(x)}{g(x)} = \frac{e^{-x^2/2}}{e^{-(x-5)^2/2}} = e^{-5x + 12.5}\] The second moment of the estimator is: \[E_g[w(X)^2 I(X \ge 5)] = \int_5^\infty e^{-10x + 25} \frac{1}{\sqrt{2\pi}} e^{-(x-5)^2/2} \, dx < \infty\] Dividing by \(p^2\) yields a bounded coefficient of variation as \(p \to 0\), proving bounded relative error. \(\blacksquare\)
B.5.3.3 Exercise 3: 2-Level Splitting Manual Walkthrough for 3-Step Random Walk
We want to estimate \(P(X_3 \ge 3.0)\) with levels \(L_1 = 1.5\), \(L_2 = 3.0\).
- Level 1: Run \(N=2\) paths of length 3:
- Path 1: \((0, 1.0, 1.2, 1.4)\) -> Max is 1.4 -> Fail.
- Path 2: \((0, 0.8, 1.6, 2.2)\) -> Max is 2.2 -> Success! Crossing occurs at \(\tau_1 = 2\).
- Level 1 probability \(p_1 = 1/2 = 0.5\).
- Level 2: Resample \(N=2\) starts from Path 2 up to \(\tau_1 = 2\): Prefix is \((0, 0.8, 1.6)\).
- Re-simulate step 3 for path 1: \(1.6 + 1.5 = 3.1 \ge 3.0\) -> Success!
- Re-simulate step 3 for path 2: \(1.6 - 0.2 = 1.4 < 3.0\) -> Fail.
- Level 2 probability \(p_2 = 1/2 = 0.5\).
- Result: Estimated probability is \(p_1 \times p_2 = 0.25\).
B.5.4 Adaptive MCMC
B.5.4.1 Exercise 1: Constant Adaptation Bias
If the proposal variance adaptation step size does not vanish (\(\gamma_t \not\to 0\)), the transition kernel remains dependent on the entire history of the chain indefinitely. This violates the Markov property, and the chain does not have a unique stationary distribution, biasing the samples away from the target \(\pi(x)\).
B.5.4.2 Exercise 2: Why Optimal Acceptance Rate is 0.234 in High Dimensions
In 1D, the optimal acceptance rate is \(0.44\) because we can easily make large moves without getting rejected. As the dimension \(d \to \infty\), the geometry of high-dimensional space restricts us; to avoid getting trapped or rejected constantly, the optimal step size scales as \(O(d^{-1/2})\), which asymptotically yields the optimal Roberts-Gelman-Gilks acceptance rate of \(0.234\).