Variance of the antithetic variates estimator with within-pair correlation

Signature

Var_anti = sigma2 / n * (1 + rho)

Function

Antithetic variates estimator of an expected model output

Computational function

• Computational function: antithetic pairs of exponential survival times with a pair-mean standard error

  Takes a hazard and a vector of uniform random numbers, one per pair, and returns the antithetic estimate of mean survival and its standard error. For each uniform number U it samples a survival time by the inverse distribution function and a partner from 1-U, averages the two, and treats the pair means as an independent sample, so it combines HE-FM-ANTV-001 with the pair-mean standard error of HE-FM-ANTV-002. The inputs therefore differ from the formula's variables: sigma2 and rho are not entered but follow from the exponential distribution, with rho equal to 1-pi^2/6.

  Inputs and outputs: h: Constant hazard of death; required, above zero. Unit: events per person-year.; U: Vector of m uniform random numbers, one per pair, each strictly above 0 and below 1. Unit: none.; T: Survival time sampled from U, -ln(U)/h. Unit: years.; T_anti: Survival time of the antithetic partner, -ln(1-U)/h. Unit: years.; Tbar: Pair mean, (T+T_anti)/2. Unit: years.; mu_anti: Mean of the pair means, the estimate of mean survival 1/h. Unit: years.; SE_anti: Sample standard deviation of the pair means divided by the root of m. Unit: years.

  Assumption: Survival times are exponential with a constant hazard and are sampled by the inverse distribution function, so U and 1-U give mirrored quantiles. Each patient has its own random number, used for the same purpose in both runs. Survival time falls as U rises, so it is monotone in U and the within-pair correlation is negative.

  Worked example (Pair drawn at U = 0.20): A uniform number of 0.20 gives a long survival time and its mirror a short one; the pair mean of 9.163 years sits closer to the true 10 years than either draw. h = 0.1; U = 0.20; T = 16.094; T_anti = 2.231; Tbar = 9.163

  Worked example (Pair drawn at U = 0.90): A uniform number of 0.90 gives a short survival time and its mirror a long one, averaging 12.040 years. h = 0.1; U = 0.90; T = 1.054; T_anti = 23.026; Tbar = 12.040

  Worked example (Two pairs combined): The two pairs above give an estimate of about 10.601 years with a standard error of about 1.438 years from the two pair means. h = 0.1; U = (0.20, 0.90); mu_anti = 10.601; SE_anti = 1.438

  Excel: =-LN(A2)/Hazard in B2, =-LN(1-A2)/Hazard in C2 and =(B2+C2)/2 in D2, filled down to row 501, with =RAND() in A2:A501 and the hazard in a cell named Hazard. =AVERAGE(D2:D501) returns mu_anti and =STDEV.S(D2:D501)/SQRT(COUNT(D2:D501)) returns SE_anti for 500 pairs. Pasting the uniform numbers as values fixes the draws so that the run can be reproduced.

  R: anti_exp <- function(u, h) { t1 <- -log(u)/h; t2 <- -log(1-u)/h; pm <- (t1+t2)/2; c(mu_anti = mean(pm), se_anti = sd(pm)/sqrt(length(pm))) } anti_exp(c(0.20, 0.90), 0.1) returns about 10.601 and 1.438, and anti_exp(runif(500), 0.1) runs 500 pairs.

  Python: def anti_exp(u, h): pm = [(-math.log(x)/h-math.log(1-x)/h)/2 for x in u]; mu = sum(pm)/len(pm); return pm, mu, statistics.stdev(pm)/math.sqrt(len(pm)) Uses the math and statistics modules, returns the pair means, the estimate and its standard error, and needs at least two pairs.

  Test (Midpoint draw gives identical antithetic partners): A uniform number of 0.5 is its own mirror image, so both members of the pair have the same survival time. Expected result: TRUE. Excel check: =ABS(-LN(0.5)/Hazard+LN(1-0.5)/Hazard)<1E-12

  Test (Sample correlation close to its theoretical value): Over 500 pairs the sample correlation between T and T_anti lies within 0.1 of 1-pi^2/6, about -0.6449. Expected result: TRUE. Excel check: =ABS(CORREL(B2:B501,C2:C501)-(1-PI()^2/6))<0.1

  Test (Pair-mean standard error below the naive value): The standard error from the 500 pair means is smaller than the standard error from all 1,000 survival times treated as independent. Expected result: TRUE. Excel check: =STDEV.S(D2:D501)/SQRT(500)<STDEV.S(B2:C501)/SQRT(1000)

  Common error (Reflecting the survival time instead of the uniform number): Mirroring a sampled survival time about its mean of 10 years, instead of applying 1-U before the transformation, gives a partner of -3.026 years for a draw of 23.026 years, which is impossible. The complement has to be taken on the uniform scale.

  Source: Owen AB. Monte Carlo Theory, Methods and Examples. 2013. Chapter 8, section 8.2, which defines the antithetic counterpart of a uniform point as one minus each coordinate and estimates the variance of the antithetic estimate from the pair means.

  T = -log(U) / h; T_anti = -log(1 - U) / h; Tbar = (T + T_anti) / 2; mu_anti = mean(Tbar); SE_anti = sd(Tbar) / sqrt(m)

Try this function

Implementations

• Excel

  Antithetic estimator variance in one cell

  With named cells Sigma2 for the variance of one output, Evaluations for the total number of model runs and PairCorr for the within-pair correlation, the formula returns the variance of the antithetic estimate.

  =Sigma2/Evaluations*(1+PairCorr)

Assumptions

• Mirrored inputs keep their distribution in antithetic pairs

  Each member of a pair has the same distribution as an ordinary draw. This holds when the complement 1-U is applied to the uniform numbers before an inverse distribution function. Generators that use rejection steps do not map one uniform number to one value in a monotone way, so feeding them the complement does not give a mirrored draw.

• Independent pairs and aligned random number streams for antithetic sampling

  The n/2 pairs are independent of each other, and both runs in a pair use each random number for the same purpose. In a discrete event model the number of random numbers consumed depends on which events occur, so a separate stream for each patient, and for each stochastic process within a patient, keeps the draws aligned.

• Sign of the within-pair correlation in antithetic sampling

  A negative rho is guaranteed when the output is monotone in every random input, with the direction free to differ between inputs. Otherwise the sign and size of rho are estimated from a pilot run of the actual model for the quantity that matters to the decision.

Worked examples

• Antithetic variance for exponential survival at 1,000 evaluations

  In the article's illustrative discrete event simulation, survival times with a hazard of 0.1 per year have variance 100 years squared and a within-pair correlation of 1-pi^2/6, about -0.6449. With 500 pairs the variance is 0.1 × 0.3551 = 0.0355, a standard error of about 0.188 years against 0.316 years for 1,000 independent patients.

  sigma2 = 100; n = 1000; rho = -0.6449; Var_anti = 0.0355

• Uncorrelated antithetic pairs match ordinary sampling

  With rho equal to zero the pairing has no effect, and the variance equals that of 1,000 independent evaluations of the same survival model, 100/1000 = 0.1.

  sigma2 = 100; n = 1000; rho = 0; Var_anti = 0.1

• Perfectly correlated antithetic pairs double the variance

  In the article's illustrative extreme for the expected value of perfect information, a V-shaped maximum gives both members of a pair the same value, so rho equals 1 and the variance doubles. The variance and run size of the survival example are reused here for comparison.

  sigma2 = 100; n = 1000; rho = 1; Var_anti = 0.2

Common errors

• Comparing antithetic pairs with independent runs at unequal cost

  The fair comparison holds the number of model evaluations fixed. Comparing 1,000 antithetic pairs, which need 2,000 evaluations, with 1,000 independent runs gives an apparent variance factor of about 0.178 instead of 0.3551 in the survival example, overstating the gain twofold.

• Expecting antithetic gains for a symmetric model output

  Outputs built from maxima, absolute values or squared deviations can be symmetric in their inputs, and pairing then doubles their contribution to the variance. In the article's extreme case for the expected value of perfect information, rho equals 1 and the variance is twice that of ordinary sampling, at a near-tie between strategies where the quantity matters most. A pilot run checks the sign of rho before the method is adopted.

• Expecting large antithetic gains when rare events drive the model

  When an event with probability p is simulated as occurring if U is below p, the two members of a pair cannot both have the event for p below 0.5, and the within-pair correlation is minus p divided by one minus p. For an illustrative annual probability of 0.02 the correlation is about -0.020 and the variance factor about 0.980, a gain of only 2%.

Sources

• Antithetic estimate and its variance in a variance reduction chapter

  Owen AB. Monte Carlo Theory, Methods and Examples. 2013. Chapter 8, Variance reduction, section 8.2 (Antithetics), equations 8.2 and 8.3, which give the antithetic estimate and its variance sigma2 (1 + rho)/n, the best and worst cases, the result for monotone functions and the even and odd decomposition.

  View source →

• Original paper introducing antithetic variates

  Hammersley JM, Morton KW. A new Monte Carlo technique: antithetic variates. Mathematical Proceedings of the Cambridge Philosophical Society. 1956;52(3):449-475. The paper that introduced antithetic variates.

  View source →

Canonical Identity