We begin by defining two widely used estimators in survey sampling: the Horvitz–Thompson (HT) estimator (Horvitz and Thompson, 1952)1¹

Also called Narain-Horvitz-Thompson estimator after Narain (1951) by Rao et al. (1999); Chauvet (2014).

An equivalent framework arises in the context of missing data, i.e., when subjects are observed with nonuniform probabilities. We follow the notations of Khan and Ugander (2023) to describe it here. Here our goal is to estimate the mean ψ of Y 1 , … , Y n , with the observations missing at random. The Bernoulli indicators R k , k = 1 , … , n indicate whether Y k , k = 1 , … , n were observed or not. We assume R k ∼ ind Bernoulli ( p k ) for k = 1 , … , n, which reflects non-response bias in sample surveys. Then, the Horvitz–Thompson and Hájek estimator of ψ = n − 1 ∑ k = 1 n Y k are: (2.4) S ˆ = ∑ k = 1 n Y k R k p k , n ˆ = ∑ k = 1 n R k p k (2.5) ψ ˆ H T = S ˆ n , ψ ˆ Hájek = S ˆ n ˆ .

A notable generalization of the aforementioned estimators is the Trotter–Tukey estimator (Trotter and Tukey, 1956). Recently, Khan and Ugander (2023) rediscovered the Trotter–Tukey estimator as the adaptive normalization (AN) idea by letting the data choose the tuning parameter λ. (2.6) ψ ˆ T T = S ˆ ( 1 − λ ) n + λ n ˆ , λ ∈ R (Trotter–Tukey) (2.7) ψ ˆ A N = S ˆ ( 1 − λ ˆ ) n + λ ˆ n ˆ . (Adaptive Normalization) Furthermore, Khan and Ugander (2023) show that the ‘adaptive normalization’ idea dating back to Trotter and Tukey (1956) leads to a lower asymptotic variance than both while generalizing these estimators, in particular, V ( ψ ˆ A N ) is lower than both V ( ψ ˆ H T ) and V ( ψ ˆ Hájek ). Khan and Ugander (2023) also provide empirical evidence that their adaptive normalization scheme leads to a lower mean squared error of IPW estimators in different application areas in average treatment effect estimation and policy learning.

Horvitz–Thompson estimators possess many desirable properties: they are unbiased, admissible and consistent. Ramakrishnan (1973) provides a simple proof that the HT estimator is admissible in a class of all unbiased estimators of a finite population total. Isaki and Fuller (1982) proves that they achieve consistency under suitable conditions, such as, specifically requiring boundedness of population and weight sequences, bounded product of inclusion probabilities and second moments, and a variance that is O ( n t − 1 ). Delevoye and Sävje (2020) provides conditions such as boundedness of moments for the outcome and inclusion probabilities and weak-design dependence. Delevoye and Sävje (2020) also point out that these conditions might be violated in ill-behaved setting with heavy-tailed outcomes or skewed sampling designs, common in natural experiments, where the practitioner has little control over the design.

It is worth noting here that IPW estimators, in particular the HT estimator, can be looked at as a weighted estimator, where the weights are not model-based, but design-based, as pointed out by several authors from Neyman (1934) to Little (2008). For example, in the HT estimator (2.2), the kth unit is assigned a weight proportional to the inverse of the selection probability as it ‘represents the 1 / p k units of the population’.

One way to look at this connection between sampling strategies and integral approximation is to represent the former as a missing data problem. Suppose that our goal is to estimate: ψ = ∫ 0 1 y ( x ) d x , which can be thought of as a limiting value of ψ n = n − 1 ∑ i = 1 n y i as n → ∞, and can be approximated up to any degree of precision by ψ N = N − 1 ∑ i = 1 N y i for a sufficiently large N. Given a random sample of size n ≪ N, with sample probabilities π attached to y i , we can view estimating θ as a problem of estimating the population quantity ψ N by ψ n . The usual importance sampling estimator in this case is akin to using E π [ ∑ i = 1 n y i / π i ] = ψ, the usual HT estimator. Just like the HT estimator, the variance of importance sampling could blow up for poor choices of g ( · ), while the bias may shrink to zero. As we show below, these connections have been exploited by several authors to propose alternative importance sampling strategies.

Given a sample ( y 1 , … , y n ), from either the density f corresponding to F itself or a suitable proposal density g, the usual importance sampling estimates ψ by the empirical average: (2.8) ψ ˆ I S = 1 n ∑ i = 1 n l ( y i ) f ( y i ) g ( y i ) . If the underlying probability measure is easy to sample from, i.e., if we can afford f = g, the above will reduce to the empirical average ψ ˆ = n − 1 ∑ i = 1 n l ( y i ), which by the Law of Large Numbers, converge to the true value of ψ at O ( n − 1 ) rate. It is worth pointing out that replacing the empirical average for the naïve importance sampling estimator in (2.8) by a Riemann sum estimator provides a remarkable improvement in stability and convergence as demonstrated in (Philippe, 1997; Philippe and Robert, 2001) or (Yakowitz et al., 1978). Recently, Datta and Polson (2025) showed that combining Riemann sum estimators with nested sampling ideas can yield sharper convergence rates (up to O ( n − 4 )) for marginal likelihood estimation in Bayesian problems. We refer the reader to Datta and Polson (2025) for a more detailed discussion of vertical likelihood, and how it connects importance sampling with nested inference. This framework also offers a unifying perspective on several sampling-based approaches that aim to improve the stability of likelihood-weighted estimators.

A natural application of IPW estimator is average treatment effect (ATE) estimation in potential outcomes causal inference framework as pointed out in (Delevoye and Sävje, 2020; Khan and Ugander, 2023). We refer the readers to the excellent references in (Cunningham, 2021; Rubin, 1974; Imbens, 2004) for in-depth coverage and historical backgrounds. Here we measure the difference in potential outcomes observed over time points t 2 > t 1 , where only one of the two potential outcomes are observed for any unit. Using Khan and Ugander (2023)’s notation: we have the triplets ( Y k ( 0 ) , Y k ( 1 ) , p k ) for k = 1 , … , n, where we observe Y k ( I k ) for I k ∼ Bernoulli ( p k ). and the parameter of interest is A T E = τ = E [ Y k ( 1 ) − Y k ( 0 ) ] from the observations available to us. The IPW estimators are employed here to estimate the two population means ψ 1 = E [ Y k ( 1 ) ] and ψ 0 = E [ Y k ( 0 ) ] separately and estimating τ by τ ˆ = ψ 1 ˆ − ψ 0 ˆ. Estimating the population means ψ j , j = 1 , 2 turns ATE into a survey sampling problem and one can use all the estimators: HT, Hájek or the various improvements such as the adaptive IPW estimator by Khan and Ugander (2023) based on the Trotter–Tukey idea (Trotter and Tukey, 1956) here. Khan and Ugander (2023) further develop an adaptive estimator by minimizing the variance of between-group differences and show that both the separate AIPW and the joint AIPW estimators achieve lower mean squared errors compared to the usual HT and Hájek.

Discussing Basu (1988), Hájek (1971) points out that the HT estimator’s “usefulness is increased in connection with ratio estimation” and proposed an estimator in presence of auxiliary information A k , related to Y k and with known total: (3.1) Y ˆ Hájek = ∑ k = 1 n A k × ∑ k = 1 n Y k p k ∑ k = 1 n A k p k , where Y k ∝ A k , k = 1 , … , n , which would not be affected in this example like the standard HT estimator. The Hájek IPW estimate for the population mean in (2.3) can be derived from the special case A k ≡ 1 , ∀ k.

Although the circus example was intended to be a pathological example, it provides at least two useful insights. First, weighted estimators can lead to nonsensical answers despite having nice large sample properties (Little, 2008). Second, problems of similar nature occur in importance sampling or Monte Carlo estimation of the marginal likelihood where empirical averages or unbiased estimator could have high or infinite variance (Li, 2010; Raftery et al., 2006). As pointed out in (Datta and Polson, 2025), strategies like Riemann sum (Philippe, 1997; Philippe and Robert, 2001), or careful choice of weight functions like vertical likelihood (Polson and Scott, 2014) or nested sampling (Skilling, 2006), or using ratio estimators (Firth, 2011), or adaptive normalization (Khan and Ugander, 2023) can resolve these issues. In particular, the same trick of ordering the draws and applying a Riemann-sum type approach is the key to avoid falling into examples like Basu’s circus.

We first re-state the example which is itself a simplification of an example from (Robins and Ritov, 1997), in a similar spirit. We consider IID samples ( Y i , X i , R i ), i = 1 , … , B such that Y i ’s are generated as a mixture of Bernoulli distributions with individual parameters indexed by the component label X i , and the ‘missingness’ indicator R i denoting whether Y i was observed or not. Let the ‘success’ probabilities associated to R i be known constants p X i satisfying: (3.2) 0 < δ ≤ p j ≤ 1 − δ < 1 , j = 1 , … , B . Note that this strong condition on p j ’s in (3.2) ensures that the HT estimator will not lead to absurd answers like Basu’s paradox stated earlier. The hierarchical model for each draw ( Y i , X i , R i ), i = 1 , … , n, is: (3.3) X i ∼ U ( 1 , … , B ) (3.4) [ R i ∣ X i = x i ] ∼ Bernoulli ( p x i ) (3.5) [ Y i ∣ R i , X i = x i ] ∼ unobserved if R i = 0 Bernoulli ( θ x i ) if R i = 1 . The parameter of interest is the average ψ = ( 1 / B ) ∑ b = 1 B θ b . Wasserman (2004) argues that since the likelihood has little information on most θ j ’s and the known constants B and p j ’s drop from the likelihood, Bayes’ estimates for ψ are going to be poor. On the other hand, the HT estimate: (3.6) ψ ˆ H T = 1 n ∑ i = 1 n R i Y i p X i is easily seen to be unbiased, given E ( R i Y i / p X i ) = E { E ( R i ∣ Y i , X i ) · E ( Y i ∣ X i ) / p X i } = E ( E ( Y i ∣ X i ) ) = ψ since E ( R i ∣ Y i , X i ) = p X i by construction. The HT estimate will also satisfy V ( ψ ˆ H T ) ≤ 1 / n δ 2 , using Hoeffding’s inequality. We would like to note here that the assumption (3.2) that the selection probabilities are between δ and 1 − δ is exploited explicitly in this asymptotic result, i.e., the HT estimator might have infinite variance as δ goes to zero, making the assumptions crucial.

This example of apparent weakness in Bayesian paradigm and the concluding remarks by Wasserman (2004)3³

“Bayesians are slaves to the likelihood function. When the likelihood goes awry, so will Bayesian inference.” (Wasserman, 2004, pp. 189)

We review the existing Bayesian resolutions for the Robins-Ritov-Wasserman problem using Bayesian ideas and argue that, in some special cases, the improvement in accuracy of Bayes’ estimate is attained via exploiting the ‘borrowing strength’ phenomenon in Stein’s shrinkage.

While a closed form analytical expression for the variance of the Li’s estimator (3.9) is difficult, we can use the linearization strategy aka the delta theorem to derive an approximation for the mean and variance of the Bayes estimator (3.9). Doing so, we provide a proof and a closer look at an assertion in Ritov et al. (2014) that Li’s estimator is consistent only if E ( Y ∣ R = 1 ) = E ( Y ), and illustrate this phenomenon via simulation studies.

Recall that for a pair of random variables X, Y with mean μ X and μ Y and variances V ( X ), V ( Y ) and covariance Cov ( X , Y ), we have the following approximations: (3.13) E ( X / Y ) ≈ μ X / μ Y (3.14) V ( X / Y ) ≈ μ X μ Y V ( X ) μ X 2 + V ( Y ) μ Y 2 − 2 Cov ( X , Y ) μ X μ Y

To calculate the approximate mean and variance for the Li’s posterior mean estimator (3.9), we first calculate the mean and variances for the numerator and denominator separately as follows. Recall once again that the posterior mean estimator is given by ψ ˆ L i = ( ∑ i R i Y i + 1 ) / ( ∑ i R i + 2 ) assuming α F = 1, i.e., a Beta ( 1 , 1 ) or Uniform prior on ψ but the exact choice of shape parameter α F does not influence the asymptotic nature of the variance. For notational convenience, we define the following quantities: θ ¯ = 1 B ∑ b = 1 B θ b ≐ ψ , p ¯ = 1 B ∑ b = 1 B p b , θ . p ¯ = 1 B ∑ b = 1 B θ b p b , σ θ , p = θ . p ¯ − p ¯ θ ¯ = θ . p ¯ − p ¯ ψ . First, the expectation for the numerator and denominator follows from applying iterated expectations as: (3.15) E ( ∑ i R i Y i + 1 ) = E X { ∑ i = 1 n E ( R i ∣ X i ) E ( Y i ∣ X i ) } + 1 = E X { ∑ i = 1 n θ X i p X i } + 1 = n θ . p ¯ + 1 = n ψ p ¯ + n σ θ , p + 1 (3.16) E ( ∑ i R i + 2 ) = n p ¯ + 2 . Hence, the expectation for the Bayes’ estimator can be approximated by taking the ratio of right hand sides from (3.15) and (3.16) as: E ( ψ ˆ L i ) ≈ n ψ p ¯ + n σ θ , p + 1 n p ¯ + 2 → ψ + σ θ , p / p ¯ , as n → ∞ , showing the ψ ˆ L i is asymptotically unbiased if σ θ , p = 0, i.e., if θ X ⊥ p X ∣ X, or equivalently if E ( Y ∣ R = 1 ) = E ( Y ), as claimed by (Ritov et al., 2014). Now, the variances for the numerator and denominator can be calculated using the conditional variance identity. The expressions are given below. (3.17) V ( ∑ i R i + 2 ) = ∑ i V ( R i ) = ∑ i E [ V ( R i ∣ X i ) ) ] + V [ E ( R i ∣ X i ) ] = ∑ i [ E ( p x i ( 1 − p X i ) ) + V ( p X i ) ] = n ( p ¯ − p ¯ 2 ) . Similarly, (3.18) V ( ∑ i R i Y i + 1 ) = ∑ i V ( R i Y i ) = n ( θ . p ¯ − θ . p ¯ 2 ) , where θ . p ¯ = 1 B ∑ b = 1 B θ b p b . Finally, the covariance term is: (3.19) Cov ∑ i R i Y i + 1 , ∑ i R i + 2 = Cov ∑ i R i Y i , ∑ i R i = ∑ i Cov ( R i Y i , R i ) = n θ . p ¯ ( 1 − p ¯ ) .

Putting everything together from (3.15), (3.16), (3.18), (3.17) and (3.19), we get the following formula for approximate variance for the Bayes estimator as: (3.20) V ( ψ ˆ L i ) ≈ n θ . p ¯ + 1 n p ¯ + 2 2 × n ( θ . p ¯ − θ . p ¯ 2 ) ( n θ . p ¯ + 1 ) 2 + n ( p ¯ − p ¯ 2 ) ( n p ¯ + 2 ) 2 − 2 n θ . p ¯ ( 1 − p ¯ ) ( n p ¯ + 2 ) ( n θ . p ¯ + 1 ) .

A couple of immediate implications of the formula in (3.20) are as follows. (i)

Li’s Bayesian estimator is consistent if θ X i and p X i are uncorrelated, as E ( ψ ˆ L i ) → ψ and V ( ψ ˆ L i ) = O ( 1 / n ) → 0 as the sample size n → ∞.