Given a vector $X=(x_1,\dots x_N) \in R_+^N$ and $N$ i.i.d Bernoulli random variable $A_i$ with $$P(A_i= 1)=p = 1 -P(A_i = 0)$$

First we define the random variable $T$ as $$ T = \frac{\sum_{i}x_i * A_i}{\sum_i A_i} $$

If we let $0/0 = 0$, then $T$ should be an unbiased estimator of the sample mean of $X$, that is $E[T] = \sum_{i} x_i/n$. I am wondering if there is any closed form formula to express the variance of $T$? Any upper bound on the variance in terms of $N$ is also helpful here.



1 Answer 1


In Short

A variation of this problem has been described by Grab and Savage (1954) 'Tables of the Expected Value of 1/X for Positive Bernoulli and Poisson Variables' JASA Vol. 49, No. 265

The following situation often arises in sampling problems. An observation $y_i \, (i=1, \cdots x)$ is made of $x$ individuals and the average, $$Y = (y_1 + y_2 + \dots + y_x)/x,$$ is computed.

If the $y_i$'s are independent observations on a random variable with mean value $\mu$ and variance $\sigma^2$, we find that the mean value of $Y$ is $\mu$. However, if $x$, the sample size, is a random variable, then the variance of $Y$ is not $\sigma^2/x$ but it is $\sigma^2E(1/x)$.


With the definition $0/0 = 0$ it becomes a bit difficult because this will give some probability that $T=0$. Below I will consider the problem as truncated such that you never have $\sum A_i =0$ or $T=0$.

Your problem is equivalent to sampling from $X$ without repetition, and taking the mean of the sample, where the sample size $n$ is determined by a binomial distribution.

If $X$ is very large in comparison to the sample size then you can estimate this by sampling with repetition.

In that approximation you can compute the variance by considering the moments conditional on the sample size $n = \sum A_n$. The variance scales like $\frac{Var(X)}{n}$ and the total variance will be a sum of those variances for different $n$ * $$Var(T) \approx Var(X) \sum_{n=1}^N n^{-1}{Pr}\left[\sum A_i = n\right] = Var(T)E(n^{-1})$$

For this expectation of $n^{-1}$, I have not an easy solution. Of course, you can simply compute it as a sum, but when $N$ becomes large you might want to have some expression to estimate it. I have created a separate question for this: Expectation of negative moment $E[k^{-1}]$ for zero-truncated Poisson distribution


Below is a comparison using the integral expression for $E(n^{-1})$ from the other question. I used $X = \lbrace{1,2,\dots,1000\rbrace}$. It illustrates that it works for smaller $p$. the discrepancy between the computation and simulation at higher $p$ is due to the approximation of sampling without repetition by sampling with repetition.

comparison of simulation with computation for different values of p

x = 1:1000
n = length(x)
p = 0.01
lambda = n*p

#### compute estimate
compVar <- function() {
  var_x = mean(x^2)-mean(x)^2
  E_inv_n = exp(-lambda)/(1-exp(-lambda))*(expint::expint_Ei(lambda)-log(lambda)-0.57721)
  ### result = 94.17509

### Simulate 

sim <- function() {
  a <- rbinom(n,1,p) ## Bernoulli variable
  result <- sum(a*x)/sum(a)

simVar <- function() {
  n_sim <- 10^4
  T_sim <- replicate(n_sim,sim())
  return(mean(T_sim^2, na.rm = TRUE)-mean(T_sim, na.rm = TRUE)^2)

### compare for different values of p

sim_v <- c()
comp_v <- c()

prng <- 10^seq(-2,-0.2,0.2)
for (ps in prng) {
  p = ps
  lambda = n*p
  sim_v <- c(sim_v,simVar())
  comp_v <- c(comp_v,compVar())

plot(prng,comp_v, log = 'xy',
     xlab = expression(P*(A[i] == 0)),
     ylab = "variance", type = "l", ylim = c(3*10^1,10^4), xlim = c(0.01,1),
     xaxt = "n", yaxt = "n")
points(prng,sim_v, pch = 21 , col = 1, bg = 0)

logaxis <- function(pos,order,tck = -0.005, las = 1) {
  axis(pos, at = 10^order, las = las)
  for (subord in order[-1]) {
    axis(pos, at = 0.1*10^subord*c(1:9),labels=rep("",9), tck = tck)

logaxis(pos = 1, order = c(-3:0))
logaxis(pos = 2, order = c(1:4), las = 2)

legend(0.01,100,c("simulation","computed"), lty = c(0,1), pch = c(21,NA))

* (note that this simple computation of the variance of the compound distribution works because the mean is the same I every condition)


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.