B.1 - Distribution of the difference of Poisson distributed counts

dPoisDiff <- function(d, lambda1, lambda2) {
   xmax = round(max(lambda1,lambda2)) + 20*sqrt(max(lambda1,lambda2))
   sum( dpois((0+d):xmax, lambda1) * dpois(0:(xmax-d), lambda2) )
}

l1 = 1
l2 = 1
d = -8:8
fd = rep(0,length(d))
for(i in 1:length(d)) fd[i] = dPoisDiff(d[i], l1, l2)
E.d  <- sum(d*fd)
E.d2 <- sum(d^2*fd)
sigma.d <- sqrt(E.d2 - E.d^2)
cat(sprintf(" d: %.3f +- %.3f ", E.d, sigma.d))
cat(sprintf(" (exact:  %.3f +- %.3f)\n", l1-l2, sqrt(l1+l2)))

barplot(fd, names=d, col='cyan', xlab='d', ylab='f(d)')

(The function dPoisDiff() is simple implementation of the reasoning shown in the text. For a more professional function see footnote .)