# Estimating sample mean from a biased sample (whose generative process is known)

I'm working on a problem where I'm trying to estimate some property of a dataset from a small non-uniform sample. (Let's just take the population mean because it's simple.)

Formally, assume I have data points $i_1, i_2 \cdots i_n$ each with an associated real-valued score $s_1, s_2, \cdots s_n$. So $i_k$ could represent different people and $s_k$ could be their height for example.

To generate an estimation sample, we are given that each data point $i_k$ also has a known probability $0 \leq p_k \leq 1$ which denotes their probability of membership in the sample. To generate a sample $S$, we start with an empty set and for $k = 1 \cdots n$, item $i_k$ is added to $S$ with probability $p_k$ and rejected otherwise. We also assume that $p_k$ is independent of $s_1 \cdots s_n$, and that the estimator knows the exact membership probabilities $p_1 \cdots p_n$.

Now, my task is to have an unbiased estimate of the population mean $$\frac{1}{n}\sum_n s_n$$ from the sample $S$. My question is, is it possible to get an unbiased sample of the population mean given $S$ and the above information? If so, how should I do this? Intuitively I think the weights of the items in $S$ have to be adjusted to reflect their differing probabilities of membership but I'm not sure if $n$ is important (it is unknown).

Edit: In this case, I don't know $n$, so are there any ways of getting around that? Horvitz-Thompson seems to require $n$ for an estimate. Also I don't necessarily know all the $p_k$, just the ones in the sample.

• It isn't clear to me how your sample is biased. The fact that your sample isn't uniform doesn't seem that relevant. If the probabilities differ & relate to $s_i$ (you say they don't), that would induce bias. But then you could just use the $p_i$'s (you say they're known) as weights & do a weighted average instead of a straight average. Commented Sep 29, 2015 at 21:48
• Indeed I think I used the word bias wrongly (I am not a statistician). Thanks for the tip. Commented Sep 29, 2015 at 22:00

Formalizing gung's suggestion, you can estimate the sample mean by inverse probability weighting, also known as the Horvitz-Thompson estimator. It is admissible in the class of unbiased estimators.

The H-T estimator can be used to estimate the sum $S = \sum_{i=1}^n y_i$ of sample values in a population using a random subsample, as well as the mean. Let's examine the sum estimator first. To model the subsampling, let $B_i \sim \text{Bernoulli}(p_i)$. Then the sum of the random subsample is $$\sum_{i=1}^n y_i B_i$$ the H-T estimator $\hat{S}$ of the population sum is $$\hat{S} = \sum_{i=1}^n y_i B_i / p_i$$ It is easy to see that $\hat{S}$ is unbiased: $$\mathbb{E}[\hat{S}] = \sum_{i=1}^n y_i \mathbb{E}[B_i] / p_i = \sum_{i=1}^n y_i p_i / p_i = S$$ To estimate the mean $S/n$ we can simply use $\hat{S}/n$ if $n$ is known. Otherwise $n$ can be estimated using inverse probability weighting once again: $$\hat{n} = \sum_{i=1}^n B_i/p_i$$ Both $\hat{S}$ and $\hat{n}$ are unbiased, but $\hat{S}/\hat{n}$ may have some bias. However, it should be small when the variance of numerator and denominator are well-controlled - for example in the large sample limit, provided the $p_i$ are not too small.

Here's some R code that shows how the H-T mean estimator works. We assume $n$ is known and compute $\hat{S}/n$, but it is easy to make it calculate $\hat{S}/\hat{n}$ instead.

n=1000
pop = 66+2*rnorm(n)
incl_prob = runif(n)

nTrial = 500
ht_est=numeric(nTrial)
for (i in 1:nTrial) {
included = as.logical(rbinom(n,1,incl_prob))
ht_est[i] = 1/n * sum(pop[included] / incl_prob[included])
}
print(paste0('population mean: ',round(mean(pop),2)))
print(paste0('average Horvitz-Thompson estimate: ',round(mean(ht_est),2)))
print(paste0('standard error in Horvitz-Thompson estimate: ',round(sd(ht_est),2)))


This code makes a single random population of 1000 subjects, subsamples from that population with a subject-dependent probability, then computes the H-T estimator. It does the subsampling & H-T estimation 500 times on the same population to help illustrate the estimator's accuracy. Here's a sample run:

[1] "population mean: 65.94"
[1] "average Horvitz-Thompson estimate: 65.9"
[1] "standard error in Horvitz-Thompson estimate: 5.09"


The first number is the population mean. The population is random, but is generated once at the beginning of the code and is fixed thereafter. Each of the 500 estimation trials takes a different random subsample from this single population.

The second number is the average of 500 Horvitz-Thompson estimates of the population mean, each from a different random subsample pop[included] of the same fixed population pop. Notice how close it is to the population mean, illustrating the unbiased property of the H-T estimator.

The third number is the standard deviation of those 500 estimates. It is an estimate of the standard error for any given H-T estimate of the population mean.

You might wonder why the average H-T estimate is so much closer to the population mean than the standard error would suggest. This is because we have averaged 500 H-T estimates together, and the error in these estimates is roughly $\sigma / \sqrt{T}$ where $\sigma$ is the standard deviation (in this case 5.09) and $T$ is the number of trials. In our code $T = 500$ so $\sigma / \sqrt{T} = 0.22$, which is on the order of the actual deviation, $0.4$, between the population mean and the 500 averaged H-T estimates.

• Great, I think this is just what I was looking for. I guess I can view it as stratified sampling with each individual item as its own strata, is that right? Commented Sep 29, 2015 at 22:01
• +1, it might help if you expanded on your explanation a little, &/or provided a simple numerical example. Commented Sep 29, 2015 at 22:08
• @gung I've added a simulation study!
– Paul
Commented Sep 29, 2015 at 23:04
• Hi Paul, thanks for the info. However in this case I don't know the value of n, the total population items. Is there some way to get around this? Commented Sep 29, 2015 at 23:33
• The population size can itself be estimated using the H-T estimator: $n \approx \sum_{i=1}^v 1/p_i$ where $v < n$ is the size of the subsample. See these notes for more info.
– Paul
Commented Sep 29, 2015 at 23:46