Sobol Sequence

Sobol Sequence in Quasi-Monte Carlo Health Economic Models

In a health economic model, a Sobol sequence is a reproducible list of parameter sets that replaces the random draws of a probabilistic sensitivity analysis, each new point falling in the gaps left by earlier ones across all inputs at once. Averaging model outputs over these points is quasi-Monte Carlo integration, which can estimate expected costs, QALYs and net monetary benefit with much less noise than as many pseudo-random numbers. The gain rests on choices that ordinary Monte Carlo simulation never faces, above all a power-of-two sample size and a random scramble that lets the error be estimated. The name is shared with Sobol indices, the variance shares of global sensitivity analysis, which are a separate idea.

Expected model outputs as integrals over the unit cube

A probabilistic analysis estimates an expectation. Fang and colleagues describe how a uniform value between 0 and 1 becomes a parameter value through the inverse cumulative distribution function, so the model output is a function on the unit cube and its mean an integral over the cube.

$$\mu = \int_{[0,1]^d} f\left(F_1^{-1}(u_1), \ldots, F_d^{-1}(u_d)\right) du, \qquad \hat{\mu}_n = \frac{1}{n} \sum_{i=1}^{n} f\left(F_1^{-1}(u_{i1}), \ldots, F_d^{-1}(u_{id})\right)$$

where $\mu$ is the expected output, $f$ is the model, $F_j^{-1}$ is the inverse cumulative distribution function of input $j$, $u = (u_1, \ldots, u_d)$ is a point in the unit cube, $u_i = (u_{i1}, \ldots, u_{id})$ is the $i$th point used, $n$ is the number of model runs and $d$ is the number of uncertain inputs.

Plain Monte Carlo integration picks the points at random, with a standard error falling with the square root of $n$. The NICE manual (PMG36, section 4.7.11) asks that enough simulations be used to minimise the effect of Monte Carlo error. Fang and colleagues note that random points leave clusters and gaps that reduce efficiency, which quasi-random sequences are designed to avoid.

Measuring evenness with discrepancy

Quasi-Monte Carlo needs a measure of how far a point set is from uniform. Owen defines the star discrepancy as the largest gap, over boxes anchored at the origin, between the share of points in a box and its volume.

$$D_n^{*} = \sup_{a \in [0,1)^d} \left| \frac{1}{n} \sum_{i=1}^{n} \mathbf{1}\left\{u_i \in [0, a)\right\} - \prod_{j=1}^{d} a_j \right|$$

where $a$ is the far corner of the box $[0, a)$, $\mathbf{1}\{\cdot\}$ is 1 for a point inside the box and 0 otherwise, and $a_j$ is the box length in dimension $j$.

Owen calls a sequence low discrepancy when $D_n^{*}$ is of order $(\log n)^d/n$, and the Koksma-Hlawka inequality turns this into an error bound.

$$\left| \hat{\mu}_n - \mu \right| \le D_n^{*} \, V_{\mathrm{HK}}(f)$$

where $V_{\mathrm{HK}}(f)$ is the total variation, in the sense of Hardy and Krause, of the model output taken as a function of $u$ on the unit cube.

The bound shows that quasi-Monte Carlo should eventually beat random sampling when the variation is finite. Owen stresses that neither term is known in practice, so the bound gives no usable interval, and that the worst case is so pessimistic at ten inputs that real errors are better explained by low-dimensional projections.

Building the points from binary digits

Owen records that Sobol published the first multidimensional digital sequences in 1967 under the name LP sequences, and that they are digital sequences in base 2. In the recipe set out by Joe and Kuo, each dimension $j$ gets a primitive polynomial over the binary field and odd integers $m_{k,j}$ extended by a bitwise recurrence.

$$x_{i,j} = i_1 v_{1,j} \oplus i_2 v_{2,j} \oplus i_3 v_{3,j} \oplus \cdots, \qquad v_{k,j} = \frac{m_{k,j}}{2^k}$$

where $x_{i,j}$ is coordinate $j$ of point $i$, $i_k$ is the $k$th binary digit of $i$ counted from the right, $v_{k,j}$ is a direction number, $m_{k,j}$ is the odd integer that defines it, and $\oplus$ is the bitwise exclusive-or of the binary expansions.

In the example of Joe and Kuo, the polynomial $x^3 + x + 1$ gives the coordinates 0, 0.5, 0.75, 0.25, 0.875 and 0.375 for the first six points. The faster Gray code version of Antonov and Saleev reorders the points only within blocks whose size is a power of two, so uniformity is unchanged.

Owen notes that adding a dimension leaves the earlier ones unchanged, and earlier dimensions are better spread. When $n$ is a power of two, each coordinate of the first $n$ points has exactly one value in every interval of width one over $n$. Some pairs of later dimensions show poor patterns, which vary with the direction numbers. Those of Bratley and Fox cover 40 dimensions, while Joe and Kuo, whose set SciPy uses, searched for good two-dimensional projections up to 21,201 dimensions.

Scrambling and replicate error estimates

An unscrambled sequence is deterministic, so it has no standard error. Randomised quasi-Monte Carlo applies a scramble that keeps the net structure while making each point uniform on the cube, and Owen shows that the estimate is then unbiased. Running $R$ independent scrambles gives an error estimate.

$$\hat{\mu}_{\mathrm{pool}} = \frac{1}{R} \sum_{r=1}^{R} \hat{\mu}_r, \qquad \mathrm{SE}\left(\hat{\mu}_{\mathrm{pool}}\right) = \sqrt{\frac{\sum_{r=1}^{R} \left(\hat{\mu}_r - \hat{\mu}_{\mathrm{pool}}\right)^2}{R(R-1)}}$$

where $\hat{\mu}_r$ is the estimate from scramble $r$ of $n$ points, $R$ is the number of scrambles and $\hat{\mu}_{\mathrm{pool}}$ is the pooled estimate, with an interval from the t distribution on $R-1$ degrees of freedom.

Owen suggests 5 or 10 replicates when a rough error estimate will do, and Fang and colleagues cite 8 to 32 as empirically sufficient. Owen reports that for smooth outputs scrambled nets reach a root mean squared error of order $n^{-3/2}$ apart from logarithmic factors.

Worked example with two and then three uncertain inputs

All numbers in this example are illustrative. Incremental net monetary benefit is $\mathrm{INMB} = \lambda p g - C$, with $\lambda$ = GBP 25,000 per QALY, one of the two values PMG36 (section 4.2.16) uses for net health benefit. First, the response probability $p$ is uniform between 0.20 and 0.60, the QALY gain per responder $g$ uniform between 0.50 and 1.00, and the cost $C$ fixed at GBP 3,000, so the true mean is 25,000 × 0.40 × 0.75 − 3,000 = GBP 4,500.

1. Take the first eight unscrambled points. These are the points shown in the SciPy documentation, with $p = 0.20 + 0.40 u_1$ and $g = 0.50 + 0.50 u_2$.

2. Check the balance. Each coordinate takes each multiple of 0.125 below 1 once, so every eighth of each range holds one point, and each quarter of the square holds two.

3. Run the model. For point 3, 25,000 × 0.50 × 0.625 = 7,812.50, less GBP 3,000, gives GBP 4,812.50. The eight values sum to GBP 30,687.50, so the estimate is GBP 3,835.94, an error of −GBP 664.06. The bias arises because the unscrambled points sit at the left ends of their intervals, so each coordinate averages 0.4375, not 0.5. A random sample of eight would have a standard error of GBP 931.7, since the variance of $p$ times $g$ is 1/90, giving an output variance of 625,000,000 × 1/90 = 6,944,444.

4. Move to a fuller model. The same formula now takes $p$ from a beta distribution Beta(8, 12), $g$ from a gamma distribution with shape 9 and scale 1/12, and $C$ from a gamma distribution with shape 16 and scale 312.5. The means are 0.40, 0.75 and GBP 5,000, so the true mean INMB is 25,000 × 0.40 × 0.75 − 5,000 = GBP 2,500. With independent inputs, the variance of $p$ times $g$ is 3/175, so the output variance is 625,000,000 × 3/175 = 10,714,286 plus 1,562,500 for $C$, a total of 12,276,786. One random sample of 1,024 runs has a standard error of $\sqrt{12{,}276{,}786/1{,}024}$ = GBP 109.5.

5. Run ten scrambles of 1,024 points. A script using SciPy 1.18 gave replicate means of 2,499.4, 2,501.1, 2,500.9, 2,500.3, 2,497.3, 2,500.3, 2,499.9, 2,497.9, 2,499.7 and 2,498.0. They sum to 24,994.8, so the pooled mean is GBP 2,499.48. The squared deviations sum to 15.656, so $\mathrm{SE} = \sqrt{15.656/90}$ = GBP 0.42, and with $t = 2.262$ the 95% interval is 2,499.48 ± 0.94. One random sample of all 10,240 runs would have a standard error of GBP 34.6.

6. Repeat for the probability of cost-effectiveness. The output is now 1 when INMB is positive and 0 otherwise, with a true mean of about 0.754. Over 500 repetitions at 1,024 runs, scrambled points cut the variance of the mean INMB about 3,600-fold against random sampling, but that of the probability only about 10-fold.

7. Interpret. This small, smooth model is close to ideal for quasi-Monte Carlo, and a large Markov or patient-level model will usually gain far less. The jump at zero INMB explains the weaker gain on the probability, since Owen notes that an indicator usually has infinite variation unless its set has edges parallel to the axes. Both designs estimate the same expected INMB.

Choices that change the result

Most errors come from habits that are harmless with random draws. Owen compares quasi-Monte Carlo with antibiotics, since in both it is best to finish the course.

• Sample size. Owen warns that 100,000 points could be much worse than 131,072, which is two to the power 17, and even worse than 65,536. SciPy warns that balance is lost with other sizes, a skipped first point or thinning, and Owen finds that skipping the first scrambled point can raise the error by a factor proportional to $\sqrt{n}$.
• Input order. Because the first coordinates are the best spread, Owen advises giving them to the most important inputs.
• The first point. The unscrambled sequence starts at the origin, where an inverse normal cumulative distribution function is infinite. Owen notes that making each point uniform, as scrambling does, avoids such singularities without knowing where they lie.
• Correlated inputs. PMG36 (section 4.7.19) asks that correlation between parameters be reflected, and Fang and colleagues note that correlated inputs can be handled with a Cholesky decomposition once the uniform points are generated.
• Adding runs. SciPy can extend a design to the next power of two, while Kucherenko and colleagues note that a Latin hypercube cannot take new points one at a time.

Evidence from health economic and test applications

Fang and colleagues tested scrambled Sobol sequences on EVPPI, where the nested calculation is costly. They report a cost of order $\varepsilon^{-2}$ for a mean squared error of $\varepsilon^2$, against $\varepsilon^{-3}$ for nested Monte Carlo. For EVPI in a simplified depression decision tree, the root mean squared error was 2.37 against 3.12 with standard Monte Carlo. They found that quasi-Monte Carlo worked less well in highly complex models such as Markov models.

For metamodelling, Degeling and colleagues note that Sobol sequences may be more efficient than Latin hypercube designs in problems with few inputs. On test functions, Kucherenko and colleagues found errors falling close to one over $n$ when a few inputs or low-order interactions dominated, smaller gains under high-order interactions, and some cases where Latin hypercube sampling did better at small sample sizes. They judged Sobol points the safest choice for functions of unknown type.

Where Sobol points sit among sampling designs

Pseudo-random draws give a simple standard error at the cost of chance clusters and gaps. Latin hypercube sampling, another form of variance reduction, stratifies each input separately but leaves the joint spread to random pairing. Fang and colleagues describe rank-1 lattice rules as the other common route to low discrepancy.

Limitations and common misreadings

The method reduces simulation noise only. It leaves parameter uncertainty and structural choices untouched.

1. Treating Sobol points as random draws. Unscrambled points are a quadrature rule, so their spread is not a sample standard error.
2. Using round sample sizes. Powers of ten or a dropped first point break the balance the method relies on.
3. Expecting the toy gain in every model. Gains shrink when many inputs interact or the output jumps, as decision indicators do.
4. Reading even coverage as less uncertainty. A cost-effectiveness acceptability curve shows the same decision uncertainty under either design, estimated more precisely.

Sources

• Degeling K, IJzerman MJ, Lavieri MS, Strong M, Koffijberg H. Introduction to metamodeling for reducing computational burden of advanced analyses with health economic models: a structured overview of metamodeling methods in a 6-step application process. Medical Decision Making. 2020;40(3):348-363. https://doi.org/10.1177/0272989X20912233
• Fang W, Wang Z, Giles MB, Jackson CH, Welton NJ, Andrieu C, Thom H. Multilevel and quasi Monte Carlo methods for the calculation of the expected value of partial perfect information. Medical Decision Making. 2022;42(2):168-181. https://doi.org/10.1177/0272989X211026305
• Joe S, Kuo FY. Notes on generating Sobol sequences. University of New South Wales; August 2008. https://web.maths.unsw.edu.au/~fkuo/sobol/joe-kuo-notes.pdf
• Kucherenko S, Albrecht D, Saltelli A. Exploring multi-dimensional spaces: a comparison of Latin hypercube and quasi Monte Carlo sampling techniques. arXiv:1505.02350; 2015. https://arxiv.org/abs/1505.02350
• National Institute for Health and Care Excellence. NICE technology appraisal and highly specialised technologies guidance: the manual (PMG36). Published 31 January 2022, last updated 13 October 2026. Sections 4.2.16, 4.7.11 and 4.7.19. https://www.nice.org.uk/process/pmg36
• Owen AB. Monte Carlo Theory, Methods and Examples. Draft chapters 15 to 17 (quasi-Monte Carlo); 2019. https://artowen.su.domains/mc/qmcstuff.pdf
• Owen AB. On dropping the first Sobol' point. arXiv:2008.08051; 2020. https://arxiv.org/abs/2008.08051
• SciPy. scipy.stats.qmc.Sobol, SciPy v1.18.0 reference guide. Accessed 28 October 2026. https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.qmc.Sobol.html