# How $Var[e^{\frac{-1}{X+a}}]$ varies with $n$ where $X \sim Bin(n,p)$?

I have a binomial random variable $$X \sim Bin(n,p)$$. I am interested in the variance of a function $$f(X)$$ given by :

$$f(X)=e^{\frac{-1}{X+a}}$$. Here $$a>0$$.

Specifically, I would like to know how $$Var[f(X)]$$ varies with $$n$$?

I am unable to come with a way to proceed. Any help is appreciated.

Since $$X$$ is a discrete random variable with finite support, it is simple to compute the moments of any function of $$X$$ using the law of the unconscious statistician. Applying this law gives:

\begin{align} m_k(n,p,a) &\equiv \mathbb{E} \Bigg( \exp \bigg( - \frac{k}{X+a} \bigg) \Bigg) \\[6pt] &= \sum_{x=0}^n \exp \bigg( - \frac{k}{x+a} \bigg) \cdot \text{Bin}(x|n, p) \\[6pt] &= \sum_{x=0}^n \exp \bigg( - \frac{k}{x+a} \bigg) \cdot {n \choose x} x^p (1-x)^{n-p}. \\[6pt] \end{align}

Using this moment function you then have:

$$\mathbb{V} \Bigg( \exp \bigg( - \frac{1}{X+a} \bigg) \Bigg) = m_2(n,p,a) - m_1(n,p,a)^2.$$

You can implement the variance function in R as shown below. This version is vectorised on the parameter n to allow you to produce the variance over a range of values. The computation is done in log-space to make it more stable for cases involving small probability values.

#Generate variance function
var.special <- function(n, p, a) {
logm1 <- rep(0, length(n))
logm2 <- rep(0, length(n))
for (i in 1:length(n)) {
x    <- 0:n[i]
logf <- -1/(x+a)
logm1[i] <- matrixStats::logSumExp(  logf + dbinom(x, n[i], p, log = TRUE))
logm2[i] <- matrixStats::logSumExp(2*logf + dbinom(x, n[i], p, log = TRUE))  }
V <- exp(logm2) - exp(2*logm1)
V }


To supplement the answer given by @Ben, consider the following plot

which is produced with the following R code.

# Moment function
m <- function(k, n, p, a){
res <- rep(NA, length(n))
for(i in seq_along(n)){
xx <- 0:n[i]
ff <- exp(-1/(xx+a))
xf <- k*log(ff) + dbinom(xx, n[i], p, log=TRUE) #For computational stability
res[i] <- exp(matrixStats::logSumExp(xf))
}
return(res)
}

# Variance function
v <- function(n, p, a){
m(2, n, p, a) - m(1, n, p, a)^2
}

# Make plot
par(mfrow=c(2,2))
nn <- 0:30
aa <- c(0.01, 0.1, 1, 10)
pp <- c(0.25, 0.5, 0.75)

for(i in 1:4){
plot(NULL, xlim=c(0, 30), ylim=c(0, max(v(nn, pp[1], aa[i]))),
xlab="sample size", ylab="Var(f(X))",
cex.lab=1.4, main=paste0("a = ", aa[i]), cex.main=1.4)
if(i == 1){
legend('topright', c("p=0.25", "p=0.50", "p=0.75"),
lty=1:3, lwd=2, bty='n', cex=1.4)
}
for(j in 1:3){
lines(nn, v(nn, pp[j], aa[i]), col=i, lty=j, lwd=2)
}
}

• Nice work (+1). However, I'd recommend amending your code for m to compute in log-space and then return to standard space at the end of the computation. That will give you a much more stable function that can compute the result for much higher values of n.
– Ben
Nov 10, 2021 at 22:26
• @Ben Good idea -- code has been updated. Nov 10, 2021 at 22:41
• Unfortunately using a sum of exponentials converts out of log-space in the intermediate calculation (which does not give stable computation). I've taken the liberty to replace sum(exp(...)) with exp(matrixStats::logSumExp(...)) to fix this. Feel free to revert if you prefer another method.
– Ben
Nov 10, 2021 at 22:44