Agricultural Statistics & Data Science · Statistical computing, data science and practice · also Statistical Computing, Computational Statistics, Advanced Statistical Computing (PhD)
The algorithms inside modern statistics: random-number generation, Monte Carlo integration, the bootstrap and jackknife, permutation tests, the EM algorithm, MCMC and numerical linear algebra, with attention to accuracy and reproducibility.
In the exams: AICE-JRF Computer Application Unit 10 (present); GATE DA (present); PhD course in most programmes; light in national exams
Modern statistics runs on algorithms. Confidence intervals for awkward statistics come from the bootstrap; p-values for small or unusual designs come from permutation tests; maximum likelihood estimates with missing data or mixtures come from the EM algorithm; Bayesian posterior distributions for hierarchical models of crop trials or animal breeding come from Markov chain Monte Carlo; and new methods are evaluated by simulation studies before anyone trusts them on real data. All of these rest on pseudo-random numbers, numerical linear algebra and, increasingly, parallel computing. A statistician who uses these tools without understanding them risks intervals that do not cover, chains that have not converged, and simulations whose conclusions are artefacts of too few replicates.
This manual covers pseudo-random number generation and tests of randomness; generating random variates by inversion, acceptance–rejection and transformation; Monte Carlo integration and variance reduction; the bootstrap, with percentile and BCa intervals; the jackknife; permutation and randomisation tests; the EM algorithm and its applications to missing data and mixtures; MCMC algorithms in practice; numerical linear algebra for statistics, especially QR and Cholesky decompositions; the design and reporting of simulation studies; and parallel and high-performance computing. The aim is to understand what each algorithm does, when it works and how to check that it has worked.
Doctoral examinations in computer applications for agriculture include numerical and computational statistics; graduate entrance examinations in data science test simulation, resampling and computational complexity; and statistical computing is a standard doctoral course in most statistics programmes. Typical questions concern linear congruential generators, the inversion method, acceptance rates, Monte Carlo standard errors, bootstrap and jackknife estimates, permutation counts, EM properties, Metropolis–Hastings acceptance ratios, effective sample size, decompositions used for least squares, and Amdahl's law.
Computers generate pseudo-random numbers: deterministic sequences that behave, for practical purposes, like independent draws from the uniform distribution on (0, 1). Determinism is an advantage, because setting the seed reproduces a simulation exactly. The classical linear congruential generator (LCG) is xₙ₊₁ = (axₙ + c) mod m, with uₙ = xₙ/m. With a = 5, c = 3, m = 16 and seed 7, the next values are (35 + 3) mod 16 = 6 and then (30 + 3) mod 16 = 1. The Hull–Dobell theorem gives conditions for an LCG with c ≠ 0 to reach the full period m: c and m are coprime, a − 1 is divisible by every prime factor of m, and a − 1 is divisible by 4 if m is; the example satisfies all three, so its period is 16. LCGs have known weaknesses: successive tuples lie on a small number of hyperplanes (Marsaglia's observation), and low-order bits have short periods, so they are unsuitable for serious simulation. Modern generators, such as the Mersenne Twister, with period 2¹⁹⁹³⁷ − 1, permuted congruential generators and counter-based generators, have long periods and good high-dimensional uniformity; for cryptographic uses, different, cryptographically secure generators are required.
Generators are tested by checking whether their output departs from properties of truly random sequences. Simple tests include the chi-square frequency test (do values fall equally into k bins?), the Kolmogorov–Smirnov test of uniformity, the serial test (are pairs or triples uniformly distributed in squares or cubes?), the runs test (are runs up and down of the expected lengths?), the gap test, the poker test and autocorrelation tests. Large batteries of tests, such as the TestU01 suite and the earlier Diehard tests, apply hundreds of such tests to very long sequences; good modern generators pass them, while many older ones fail. A single test cannot prove randomness; it can only fail to detect non-randomness. In practice, statisticians should use the well-tested default generator of their software, record the generator and seed, and, in parallel computing, use streams designed to be independent rather than consecutive seeds.
If U is uniform on (0, 1) and F is a continuous distribution function, then X = F⁻¹(U) has distribution F (the probability integral transform). For the exponential distribution with rate λ, F(x) = 1 − e^(−λx), so X = −ln(1 − U)/λ, or equivalently −ln(U)/λ since 1 − U is also uniform; with λ = 2 and U = 0.5, X ≈ 0.347. Inversion works for any distribution whose quantile function is available: the Weibull, Pareto, logistic and Cauchy (X = tan[π(U − ½)]) have closed forms. For discrete distributions, the generalised inverse is used: generate U and return the smallest x with F(x) ≥ U, searching the cumulative probabilities, which for a Bernoulli(p) variable means returning 1 if U < p. Inversion preserves the order of uniforms, which is useful for variance reduction and for generating correlated variables, but numerical inversion can be slow for distributions such as the normal and gamma, for which other methods are often used.
When F⁻¹ is not available, acceptance–rejection samples from a target density f using a proposal density g that is easy to sample and satisfies f(x) ≤ c g(x) for all x and some constant c ≥ 1. Draw Y from g and U uniform on (0, 1); accept X = Y if U ≤ f(Y)/[c g(Y)], and otherwise repeat. Accepted values have density f exactly. The probability of acceptance on each try is 1/c, so the number of tries per accepted value is geometric with mean c: with c = 1.5, two-thirds of proposals are accepted and on average 1.5 tries are needed. Efficiency therefore requires a proposal close in shape to the target, with tails at least as heavy. Squeeze methods add easily computed bounds to avoid evaluating f often, and adaptive rejection sampling builds envelopes for log-concave densities, which made it useful inside early Gibbs samplers.
Many variates are produced by transforming others. The Box–Muller method turns two independent uniforms U₁ and U₂ into two independent standard normals: Z₁ = √(−2 ln U₁) cos(2πU₂) and Z₂ = √(−2 ln U₁) sin(2πU₂); the polar (Marsaglia) method avoids the trigonometric functions. Other relationships: a chi-square variable with k degrees of freedom is a sum of k squared standard normals; a gamma variable with integer shape k is a sum of k exponentials; if X and Y are independent gammas with the same rate, X/(X + Y) is beta; a t variable is Z/√(V/k) with V chi-square; Poisson counts can be generated by counting exponential inter-arrival times in a unit interval; and a binomial is a sum of Bernoulli variables. Multivariate normal vectors with mean μ and covariance Σ are generated as μ + LZ, where L is the Cholesky factor of Σ (Σ = LL′) and Z a vector of independent standard normals. Mixture distributions are sampled by first choosing a component with its mixing probability and then drawing from that component.
In practice statisticians rarely write their own generators for standard distributions: the built-in functions of statistical software use carefully tested algorithms, often specialised acceptance–rejection or ratio-of-uniforms methods for the normal, gamma, beta and binomial, chosen for speed and accuracy across all parameter values. Writing a custom generator is justified for non-standard distributions, such as a fitted mixture, a truncated distribution or a posterior known up to a constant. Truncated distributions are conveniently sampled by inversion restricted to the truncation interval, X = F⁻¹(F(a) + U[F(b) − F(a)]), which never wastes draws, whereas naive rejection of values outside [a, b] becomes very inefficient when the interval has small probability. Every custom generator should be checked by comparing a large sample with the target distribution, using histograms, quantile plots and goodness-of-fit tests, before it is used in a simulation study.
Monte Carlo integration estimates an expectation θ = E[h(X)] = ∫h(x)f(x)dx by the sample mean θ̂ = (1/n)Σh(Xᵢ) of independent draws from f. By the law of large numbers θ̂ converges to θ, and by the central limit theorem its standard error is σ/√n, where σ² = Var[h(X)] is estimated from the draws. The error therefore falls as 1/√n regardless of the dimension of the integral, which is why Monte Carlo methods dominate in high dimensions, but halving the error needs four times as many draws. A classic example estimates π: of n points uniform in the unit square, the proportion falling inside the quarter circle estimates π/4; with 7,850 of 10,000 inside, π̂ = 3.14, and the standard error is 4√(0.785 × 0.215/10,000) ≈ 0.016. Monte Carlo is used to compute probabilities and expectations under complex models, power of tests, distributions of statistics and posterior quantities, and every Monte Carlo estimate should be reported with its Monte Carlo standard error.
Variance reduction techniques obtain the same precision with fewer draws. Antithetic variates pair each uniform U with 1 − U (or each normal Z with −Z); if h is monotone, the paired values are negatively correlated, and their average has lower variance than two independent values. Control variates use a variable C correlated with h(X) whose mean is known, and estimate θ by the mean of h(X) − β(C − E[C]); with the optimal β = Cov(h, C)/Var(C), the variance is reduced by the factor 1 − ρ², where ρ is their correlation, so a control with ρ = 0.9 reduces the variance by 81%. Stratified sampling divides the range of the input into strata and samples each in proportion, removing the between-strata component of variance. Importance sampling draws from a different density g and reweights by f/g: θ̂ = (1/n)Σh(Xᵢ)f(Xᵢ)/g(Xᵢ); it is essential for estimating small tail probabilities, such as the probability of a rare extreme event, where plain Monte Carlo would see few or no successes, but a poorly chosen g with lighter tails than f can make the variance infinite. Common random numbers, using the same random stream to compare two systems or methods, reduce the variance of their difference. Quasi-Monte Carlo replaces random points with low-discrepancy sequences, such as Sobol' sequences, which fill space more evenly and can converge faster for smooth integrands in moderate dimensions.
The bootstrap, introduced by Efron in 1979, estimates the sampling distribution of a statistic by resampling the observed data. In the non-parametric bootstrap, B samples of size n are drawn with replacement from the n observations, the statistic θ̂* is computed for each, and the distribution of the θ̂* values approximates the sampling distribution of θ̂. The bootstrap standard error is the standard deviation of the θ̂*, and the bootstrap bias estimate is the mean of the θ̂* minus θ̂. Each bootstrap sample omits some observations: the probability that a given observation is left out of a sample of size n is (1 − 1/n)ⁿ, which tends to e⁻¹ ≈ 0.368, so about 63.2% of distinct observations appear in each resample. A few hundred resamples usually suffice for a standard error, but confidence intervals, which depend on the tails, need one or two thousand or more. The parametric bootstrap instead simulates from a fitted model, which is more efficient when the model is right. Resampling must follow the structure of the data: for a designed experiment, residuals within the model or whole blocks may be resampled; for clustered data such as plots within farms, clusters are resampled; and for time series, blocks of consecutive observations are resampled to preserve dependence.
The percentile interval takes the α/2 and 1 − α/2 quantiles of the bootstrap distribution, for example the 2.5th and 97.5th percentiles for a 95% interval. It is simple and transformation-respecting (the interval for g(θ) is g applied to the interval for θ) but can have poor coverage when the statistic is biased or its variance depends on the parameter. The normal interval θ̂ ± z × SE_boot assumes approximate normality. The basic (reverse percentile) interval reflects the percentiles around θ̂. The bootstrap-t interval studentises each resample using its own standard error and is second-order accurate when a good standard error is available. The bias-corrected and accelerated (BCa) interval adjusts the percentiles for median bias (through the proportion of θ̂* below θ̂) and for skewness or changing variance (through an acceleration constant, usually estimated by the jackknife); it is second-order accurate and transformation-respecting, and is generally recommended for routine use. The bootstrap fails or needs modification for some statistics, such as the maximum of a sample, and for very small samples, where the empirical distribution poorly represents the population.
For regression and designed experiments there are two main bootstrap schemes. Case (pairs) resampling draws whole observations, the response together with its covariates, and is robust to heteroscedasticity and model misspecification, but it can produce resamples in which some treatment levels are missing. Residual resampling fits the model, resamples the residuals and adds them to the fitted values to create new responses, keeping the design fixed; it is more efficient when the model and constant-variance assumption hold, and the wild bootstrap, which multiplies each residual by a random sign or weight, retains heteroscedasticity. For mixed models with blocks or sites, resampling must happen at the level of the random effects, for example resampling whole sites with their plots. Bootstrap hypothesis tests require resampling from a distribution that satisfies the null hypothesis, for instance by centring the data at the hypothesised mean before resampling; resampling the raw data and checking whether the interval excludes the null value is the interval-inversion alternative. Bootstrap estimates of prediction error, such as the .632 estimator, adjust for the optimism of evaluating a model on the data that fitted it.
The jackknife, introduced by Quenouille for bias reduction and extended by Tukey to variance estimation, computes the statistic n times, each time leaving out one observation, giving θ̂₍ᵢ₎ for i = 1, …, n, with mean θ̂₍.₎. The jackknife estimate of bias is (n − 1)(θ̂₍.₎ − θ̂), and the bias-corrected estimate is θ̂ − bias = nθ̂ − (n − 1)θ̂₍.₎. The jackknife standard error is √[((n − 1)/n) Σ(θ̂₍ᵢ₎ − θ̂₍.₎)²]; the factor (n − 1)/n, rather than 1/(n − 1), compensates for the leave-one-out estimates being much closer to each other than independent estimates would be. For the sample mean, the jackknife standard error equals the usual s/√n exactly. Pseudo-values, θ̃ᵢ = nθ̂ − (n − 1)θ̂₍ᵢ₎, behave approximately like independent observations, and their mean and variance give the bias-corrected estimate and its variance.
The jackknife needs only n recomputations, fewer than the bootstrap, and is deterministic. It works well for smooth statistics, such as means, variances, correlations, ratios and regression coefficients, and is used in survey sampling for variance estimation of complex estimators, where the delete-a-group jackknife leaves out whole primary sampling units or groups of them, and in species richness estimation (jackknife estimators of the number of species). It fails for non-smooth statistics such as the median and other quantiles, where the leave-one-out values change in jumps and the variance estimate is inconsistent; the delete-d jackknife, leaving out d observations at a time with d growing with n, or the bootstrap, is used instead. The jackknife also supplies the acceleration constant in BCa bootstrap intervals and influence values that show which observations most affect an estimate. Leave-one-out cross-validation, used to estimate prediction error, follows the same idea.
A permutation test obtains the null distribution of a test statistic by rearranging the data in all the ways that are equally likely under the null hypothesis. To compare two treatment groups, if treatments have no effect, the labels attached to the observed values are arbitrary, so every reassignment of labels to values is equally likely. The test statistic, such as the difference in means, is computed for each reassignment, and the p-value is the proportion of reassignments giving a statistic at least as extreme as the observed one. With two groups of five, there are C(10, 5) = 252 distinct reassignments, so the smallest attainable one-sided p-value is 1/252 ≈ 0.004 and the smallest two-sided p-value is 2/252. When complete enumeration is impractical, a random sample of B permutations gives a Monte Carlo approximation, and the p-value is computed as (number of permuted statistics at least as extreme + 1)/(B + 1), counting the observed arrangement, which keeps the test valid.
In a randomised experiment, Fisher's randomisation argument gives the permutation test a design basis: the reference set consists of the treatment allocations that the randomisation could have produced, so the test needs no assumption of normality or random sampling from a population. The permutations must respect the design: in a randomised complete block design, treatments are permuted only within blocks; in a split-plot, main-plot treatments are permuted among main plots and subplot treatments within main plots; and in a completely randomised design, all allocations are permuted. Permuting freely across blocks would ignore the restriction and give a wrong reference distribution. Permutation tests are exact under the null hypothesis of no treatment effect on any unit, robust to non-normal errors and outliers when a robust statistic is used, and applicable to any statistic. Their limitations are that they test a specific null of exchangeability; when groups have unequal variances, the null of equal distributions is not the same as equal means, and the test may be invalid for the means question. Permutation methods are also used for testing correlation (permuting one variable), for multivariate tests such as PERMANOVA and the Mantel test, and for controlling family-wise error across many tests with maximum-statistic methods.
The expectation–maximisation (EM) algorithm, formalised by Dempster, Laird and Rubin in 1977, finds maximum likelihood estimates when data are incomplete, either genuinely missing or augmented with unobserved latent variables that make the likelihood simple. Each iteration has two steps. In the E-step, the expected value of the complete-data log-likelihood is computed, given the observed data and the current parameter estimate θ⁽ᵗ⁾; this function is called Q(θ | θ⁽ᵗ⁾). In the M-step, Q is maximised over θ to give θ⁽ᵗ⁺¹⁾. For exponential-family complete-data models, the E-step reduces to computing expected sufficient statistics and the M-step to the usual complete-data estimates using them. The key property is that the observed-data likelihood never decreases from one iteration to the next, so EM is stable and needs no step-size choices. Its drawbacks are that convergence is only linear, and can be very slow when the proportion of missing information is large; that it may converge to a local maximum or saddle point, so several starting values should be tried; and that it does not directly give standard errors, which are obtained by methods such as Louis's formula, the supplemented EM algorithm or the bootstrap. Variants include generalised EM (the M-step only increases Q), ECM (conditional maximisation in steps) and Monte Carlo EM (the E-step by simulation).
With missing data, EM is valid when data are missing at random (MAR), meaning that the probability of being missing depends only on observed values. For a multivariate normal sample with some values missing, the E-step replaces missing values' sufficient statistics (sums and cross-products) by their conditional expectations given the observed values, including the conditional variances, and the M-step recomputes means and covariances; simply filling in predicted values and treating them as observed would understate variability. In finite mixture models, such as a mixture of two normal distributions describing a population of plants from two genotypes or fish from two cohorts, the latent variable is each observation's component membership. The E-step computes each observation's posterior probability of belonging to each component (its 'responsibility') given current parameters, and the M-step updates the mixing proportions as the average responsibilities and the component means and variances as responsibility-weighted means and variances. The same structure underlies estimation in hidden Markov models (the Baum–Welch algorithm), latent class models, random-effects models and the estimation of allele frequencies when genotypes are not fully observed. In designed experiments, EM gives a principled version of the classical method of estimating missing plot values so that the analysis of variance can proceed.
EM gives maximum likelihood estimates under a specified model, but many analyses need a completed data set that can be analysed by standard methods. Multiple imputation, developed by Rubin, creates m completed data sets by drawing missing values from their predictive distribution given the observed data, with parameter uncertainty propagated, for example by Bayesian draws or by bootstrapping before imputing. Each data set is analysed in the usual way, and the m estimates are combined by Rubin's rules: the overall estimate is their mean, and its variance is the average within-imputation variance plus (1 + 1/m) times the between-imputation variance. Imputing a single predicted value, by contrast, treats guesses as data and understates standard errors. Like EM, multiple imputation is valid under missingness at random; neither method fixes data missing for reasons related to the unobserved values themselves, which require sensitivity analysis with explicit assumptions about the missingness mechanism.
Markov chain Monte Carlo (MCMC) draws dependent samples from a distribution π, typically a Bayesian posterior known only up to a constant, by constructing a Markov chain whose stationary distribution is π. In the Metropolis–Hastings algorithm, from the current state x a candidate y is drawn from a proposal q(y | x) and accepted with probability α = min{1, [π(y)q(x | y)]/[π(x)q(y | x)]}; otherwise the chain stays at x. The unknown normalising constant cancels in the ratio. For a symmetric proposal, such as a random walk y = x + ε with ε normal, the ratio reduces to π(y)/π(x) (the original Metropolis algorithm): moves to higher density are always accepted, and moves to lower density with probability equal to the density ratio, so a move to a point with half the density is accepted half the time. The proposal scale must be tuned: very small steps are almost always accepted but explore slowly, and very large steps are mostly rejected. Theory for random-walk proposals suggests acceptance rates of about 0.44 in one dimension and about 0.23 in high dimensions. The Gibbs sampler updates each parameter, or block of parameters, in turn by drawing from its full conditional distribution given all the others; it is a special case of Metropolis–Hastings in which every draw is accepted, and it is convenient for hierarchical models with conjugate structure, such as the variance-component models of animal and plant breeding. Hamiltonian Monte Carlo uses gradients of the log-density to make long, efficient moves and underlies modern probabilistic programming software.
MCMC output must be checked before use. The initial part of the chain, the burn-in or warm-up, which depends on the starting values, is discarded. Trace plots of parameters against iteration should look like stationary noise without trends. Running several chains from dispersed starting points and comparing within-chain and between-chain variances gives the Gelman–Rubin statistic R̂, which should be close to 1; values above 1.01 to 1.1, depending on the convention, indicate lack of convergence. Successive draws are autocorrelated, so the chain contains less information than an independent sample of the same length; the effective sample size, ESS = N/(1 + 2Σρₖ), where ρₖ are the autocorrelations at lag k, measures this: 10,000 draws with Σρₖ = 4.5 are worth about 1,000 independent draws. Monte Carlo standard errors of posterior means should be based on the ESS. Thinning, keeping every k-th draw, saves storage but does not improve precision for a given run length. Poor mixing is often due to strong posterior correlations between parameters and can be improved by reparameterisation, such as centring covariates or using non-centred parameterisations for hierarchical models, and by block updates. Convergence diagnostics can reveal problems but never prove convergence.
Statistical formulas are written with matrix inverses, such as β̂ = (X′X)⁻¹X′y and the generalised least-squares estimator (X′V⁻¹X)⁻¹X′V⁻¹y, but explicitly inverting matrices is slower and less accurate than solving systems through factorisations. Computers store numbers in floating point with about 16 significant digits in double precision, and every operation rounds, so algorithms are judged by numerical stability, how much rounding error they let accumulate, as well as speed. The condition number κ of a matrix bounds the amplification of relative errors; forming X′X squares the condition number of X, so a regression whose X has κ = 10⁶ leads to normal equations with κ = 10¹², leaving only about four reliable digits. This is the main reason for using orthogonal decompositions.
The QR decomposition writes an n × p matrix X of full column rank as X = QR, where Q (n × p) has orthonormal columns and R (p × p) is upper triangular. Then the least-squares problem reduces to Rβ = Q′y, solved by back substitution, without forming X′X; the residual sum of squares is ‖y − QQ′y‖², since the fitted values are QQ′y. QR is computed by Householder reflections, the standard and most stable method, by Givens rotations, which suit sparse matrices and updating, or by modified Gram–Schmidt. Because (X′X)⁻¹ = R⁻¹(R⁻¹)′, the variance–covariance matrix of the estimates is obtained from the small triangular R. Column pivoting in QR detects rank deficiency, identifying aliased columns in, for example, over-parameterised factorial models, which is how standard regression software reports coefficients 'not defined because of singularities'.
For a symmetric positive definite matrix A, the Cholesky decomposition A = LL′, with L lower triangular and positive diagonal, exists and is unique, costs about n³/3 operations (half that of LU) and is stable without pivoting. It solves systems Ax = b by two triangular solves, gives log-determinants as 2Σ ln Lᵢᵢ, needed in multivariate normal and mixed-model likelihoods, generates multivariate normal variates as μ + LZ, and solves normal equations when X is well conditioned. Failure of the Cholesky algorithm signals that a matrix is not positive definite, a useful check on estimated covariance matrices. Other tools include the eigendecomposition and singular value decomposition for principal components, rank determination and pseudo-inverses; sparse factorisations for the very large, sparse mixed-model equations of breeding evaluations; and iterative methods such as conjugate gradients. Numerical care also includes working on the log scale for likelihoods, using the log-sum-exp trick, log Σ exp(aᵢ) = m + log Σ exp(aᵢ − m) with m = max aᵢ, to avoid overflow and underflow, and computing variances by stable algorithms.
Simulation studies evaluate statistical methods by applying them to many data sets generated from known models, so that estimates can be compared with the truth. Morris, White and Crowther (2019) set out a structured plan summarised as ADEMP: Aims (what question the study answers, such as whether a mixed-model test controls type I error with few blocks); Data-generating mechanisms (the models, parameter values, sample sizes and departures from assumptions to be simulated, usually in a factorial arrangement); Estimands (the targets of inference, such as a treatment effect); Methods (the procedures to be compared); and Performance measures (how they will be judged). Data-generating mechanisms should include realistic scenarios drawn from the application, such as the sizes, variance components and missingness patterns of actual crop trials, and also challenging ones that probe the method's limits.
Common performance measures are bias (mean of estimates minus the true value), empirical standard error (standard deviation of the estimates), mean squared error (bias² plus variance), the average model-based standard error compared with the empirical one, coverage of confidence intervals (proportion containing the true value), type I error and power of tests, and convergence failure rates. Every performance measure is itself estimated from a finite number of replicates n_sim and has a Monte Carlo standard error, which must be reported. For coverage with true probability p, the Monte Carlo SE is √[p(1 − p)/n_sim]: with p = 0.95 and 1,000 replicates, about 0.0069, so an observed coverage of 0.94 is not clear evidence of under-coverage. The number of replicates is chosen to achieve a required precision: for a coverage SE of 0.005, n_sim = 0.95 × 0.05/0.005² = 1,900. For bias, the Monte Carlo SE is the empirical SE divided by √n_sim.
Simulations should set and record random seeds, store the generated data or the states of the generator so that results can be reproduced and failures investigated, use independent random streams for parallel runs, handle non-convergence explicitly (reporting how often it occurred and how it was treated rather than silently dropping failed replicates), and use the same simulated data sets for all methods compared, which reduces the Monte Carlo error of differences. Results are best presented graphically, with performance measures and their Monte Carlo SEs across scenarios, and the code should be shared. Simulation studies by method developers tend to favour their own methods, so neutral comparison studies with pre-specified designs are valuable.
Many statistical computations consist of independent tasks: bootstrap resamples, permutation replicates, simulation replicates, cross-validation folds, separate MCMC chains, and analyses of many traits, sites or genes. These are 'embarrassingly parallel' and can run simultaneously on several processor cores or machines with almost no communication, giving speed-ups close to the number of workers. Other computations, such as a single large matrix factorisation or a single MCMC chain, are harder to parallelise. Amdahl's law limits the speed-up from N workers when only a fraction p of the work can be parallelised: speed-up = 1/[(1 − p) + p/N]. With p = 0.9 and eight cores, the speed-up is 1/(0.1 + 0.1125) ≈ 4.7, and no number of cores can exceed 1/(1 − p) = 10. So the serial part, often data input, setup and combining results, ultimately limits performance.
Parallel computing uses multicore processors on a single machine (shared memory), clusters of machines connected by a network (distributed memory, often programmed with message-passing interfaces or higher-level frameworks), and graphics processing units (GPUs), which have thousands of simple cores suited to the same arithmetic on many data, such as matrix operations and the training of neural networks. High-performance computing clusters are shared through job schedulers, to which users submit batch jobs specifying cores, memory and time; job arrays run many similar tasks, such as simulation scenarios. Statistical languages provide parallel tools that distribute independent function calls over cores or nodes, and optimised numerical libraries already use multiple cores for linear algebra. Cloud computing provides such resources on demand.
Parallel random number generation requires streams that are independent; giving each worker a consecutive seed of a simple generator can produce correlated streams, whereas generators with designed parallel streams or jump-ahead avoid this, and results remain reproducible only if the assignment of streams to tasks is fixed. Overheads of starting workers and moving data can exceed the gains for small tasks, so work should be divided into chunks large enough to amortise them. Memory can be the binding constraint: copying a large data set to every worker may exhaust memory. Before parallelising, code should be profiled and vectorised, since a better algorithm or vectorised code often gives larger gains than more cores. Finally, results from parallel runs should be identical to serial runs with the same streams, which is a useful test.
Computational results should be reproducible by others and by the analyst at a later date. That requires recording the software versions, the random number generator and its seeds, the hardware-dependent settings that matter, and the exact code and data used, preferably in a version-controlled repository with a project environment that fixes package versions. Floating-point results can differ slightly between machines, compilers and numbers of threads, because the order of summation changes rounding; results of simulations and MCMC runs should therefore be judged by their Monte Carlo error rather than by exact agreement, while deterministic analyses should agree to the precision reported. Containers and workflow tools that rerun the full pipeline from raw data to tables and figures make reproduction routine, and they protect long-running studies, such as multi-year trial series, from silent changes when software is updated.
A postgraduate student uses this subject to write simulation code with seeds, generate variates from non-standard distributions, compute bootstrap and BCa intervals for ratios and heritabilities, run permutation tests that respect the randomisation of a field trial, implement EM for a mixture of two normal distributions, fit a simple Bayesian model by Metropolis–Hastings and check its convergence, and design a small simulation study with Monte Carlo standard errors. Examination questions test the formulas and properties: LCG periods, inversion, acceptance rates, (1 − 1/n)ⁿ, jackknife SEs, permutation counts, EM monotonicity, acceptance ratios, ESS and Amdahl's law.
A statistician five years into practice uses resampling to give uncertainty for quantities without easy formulas, such as ratios of yields or stability indices; uses permutation tests for small experiments with skewed responses; fits mixed and Bayesian models with MCMC and checks diagnostics before reporting; handles missing data with likelihood or EM-based methods; and parallelises long simulations on a cluster. Typical decisions are how many resamples or iterations are enough, how to resample clustered data, and whether a computational result is trustworthy.
A senior statistician with twenty-five years of experience sets computational standards for a group: approved generators and seeding practices, required diagnostics for MCMC, minimum resample numbers, documentation of simulation studies, and code review. They decide when computational methods are acceptable in regulatory or official reports, invest in computing infrastructure, and train colleagues to recognise failed convergence and unreliable resampling. They decide when a computational shortcut is acceptable and when the full method must be run, and they ensure that critical results are reproduced independently before publication.
A researcher with a PhD and long experience develops computational methods: efficient MCMC and variational algorithms for large hierarchical models in genomics and spatial statistics, bootstrap methods for complex survey and dependent data, scalable EM variants, and well-designed simulation studies that compare methods fairly. Such researchers write and maintain software implementing these methods, review computational claims in papers, and teach doctoral courses in statistical computing. They also mentor students in writing efficient, tested code for new methods.
This manual explains the main algorithms and their properties; it does not cover sequential Monte Carlo, variational inference, approximate Bayesian computation, stochastic optimisation or the details of Hamiltonian Monte Carlo tuning. Thresholds such as acceptance-rate targets, R̂ cut-offs and resample numbers are conventions and guidelines, not exact rules. Software-specific tools for parallel computing are described generically.
Original practice questions written to the syllabus and pattern of the exams named in each question's examTags. They are not past papers and do not reproduce any board's questions.
These notes are free to read here, with no account.
What the Laisha app adds:
Free for ten days, then a subscription. We would rather write that here than have you find it later.