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.
 A: 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 }

A: 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)
  }
}

