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.


2 Answers 2


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 enter image description here

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))

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

# Make plot
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)
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ben
    Nov 10, 2021 at 22:26
  • 1
    $\begingroup$ @Ben Good idea -- code has been updated. $\endgroup$
    – knrumsey
    Nov 10, 2021 at 22:41
  • 2
    $\begingroup$ 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. $\endgroup$
    – Ben
    Nov 10, 2021 at 22:44

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.