On/Off Problem

In a standard type of experiment in physics and astronomy we have the following situation. We observe \(X\sim Pois(\lambda+\mu)\), where \(\lambda\) is a known background rate and \(\mu\) is the signal rate of interest. We want to find an interval estimate for \(\mu\). As a numeric example let’s say we have \(\lambda=10\) and \(x=27\).

Let’s do a Bayesian analysis, and use a flat prior on \(\mu\). Of course \(\mu\ge 0\), so

\[ \begin{aligned} &g(\mu)= I_{[0, \infty)}(\mu)\\ &f(x|\mu) = \frac{(\lambda+\mu)^x}{x!}e^{-\lambda-\mu}\\ &f(x,\mu) = \frac{(\lambda+\mu)^x}{x!}e^{-\lambda-\mu}I_{[0, \infty)}(\mu)\\ &f(\mu|x) \propto \frac{(\lambda+\mu)^x}{x!}e^{-\lambda-\mu}I_{[0, \infty)}(\mu) \end{aligned} \] Now we want to generate data from this distribution. Let’s try accept-reject with a normal density mean and standard deviation the same as the Poisson:

f <- function(mu) dpois(27, mu+10)
g <- function(mu) dnorm(mu, 27-10, 27-10)
curve(f(x)/g(x), 0, 50)

It seems this function has a maximum, so maybe we can use accept-reject! Now

onoff.mu <- function(B=1e4, x=27, lambda=10) {
  f <- function(t) dpois(x, t+lambda)
  g <- function(t) dnorm(t, x-lambda, x-lambda)
  M <- (-optimize(function(x) {-f(x)/g(x)}, 
              c(0, 2*(x-lambda)))$objective)
  mu <- rep(0, B)
  for(i in 1:B) {
    repeat {
      y <- rnorm(1, x-lambda, x-lambda)
      if(y<0) next
      u <- runif(1, 0, M*g(y))
      if(u<f(y)) {mu[i] <- y; break}
    }
  }
  mu
}
mu <- onoff.mu()

hist(mu, 50, freq = FALSE, main="")

quantile(mu, c(0.025, 0.975))

##      2.5%     97.5% 
##  8.632818 29.210377

Now let’s say we don’t know \(\lambda\), but we can measure \(Y\sim Pois(\lambda)\) (the “Off” region measurement). Say we observe \(Y=10\). Again as a prior on \(\lambda\) we use the uniform. Now the density becomes

\[ f(x,y|\mu,\lambda) = \frac{(\lambda+\mu)^x}{x!}e^{-\lambda-\mu}\frac{\lambda^y}{y!}e^{-\lambda} \]

and so the joint posterior density of \((\mu,\lambda)|X=x,Y=y\) is given by

\[ f(\mu,\lambda|x,y) \propto \frac{(\lambda+\mu)^x\lambda^y}{x!y!}e^{-2\lambda-\mu} \]

Let’s try Metropolis Hastings with a simple random walk proposal:

onoff.MH <- function(B=1e5, eps, x=27, y=10) {
  f <- function(l, m) dpois(x, l+m)*dpois(y, l)
  lm <- matrix(0, B, 2)
  lm[1, ] <- c(y, x-y)
  for(i in 2:B) {
    u <- max(0, lm[i-1, 1]+runif(1, -eps, eps))
    v <- max(0, lm[i-1, 2]+runif(1, -eps, eps))
    if(runif(1)<f(u, v)/f(lm[i-1, 1], lm[i-1, 2])) 
      lm[i, ] <- c(u, v)
    else
      lm[i, ] <- lm[i-1, ]
  }
  lm
}
lm <- onoff.MH(eps=0.5)
showit <- function(x) {
  B <- nrow(x)
  par(mfrow=c(2, 2))
  plot(cumsum(x[, 1])/1:B, 
       type="l", xlab="x", ylab="")
  plot(cumsum(x[, 2])/1:B, 
       type="l", xlab="y", ylab="")
  hist(x[(B/2):B, 1], 50, 
       freq = FALSE, main="", xlab="x")
  hist(x[(B/2):B, 2], 50, 
       freq = FALSE, main="", xlab="y")
  round(quantile(x[(B/2):B, 2], c(0.025, 0.75)), 3)
}
showit(lm)

##   2.5%    75% 
##  3.865 19.135

Often in these cases we also have the problem of efficiency, that is not all events that actually happen are indeed observed. So if the efficiency is (say) 75% and we have 20 observations, the actual number likely was

round(20/0.75, 1)

## [1] 26.7

Unfortunately the exact efificency is also not known, and often it is modeled as a normal distribution where the standard deviation is a known fraction of the mean. So now our model becomes

\[ f(\mu,\lambda, \epsilon|x,y,z) \propto \frac{(\lambda+\epsilon \mu)^x\lambda^y}{x!y!}e^{-2\lambda-\epsilon \mu}\exp\{-(\epsilon-z)^2/(pz)^2\} \] say we have \(z=0.75,p=1/10\), then

onoff1.MH <- function(B=1e5, eps, x=27, y=10, z=0.75) {
  f <- function(l, m, e) 
    dpois(x, l+e*m)*dpois(y, l)*dnorm(e, z, z/10)
  lm <- matrix(0, B, 3)
  lm[1, ] <- c(y, (x-y)/z, z)
  for(i in 2:B) {
    u <- max(0, lm[i-1, 1]+runif(1, -eps, eps))
    v <- max(0, lm[i-1, 2]+runif(1, -eps, eps))
    w <- min(1, 
        max(0, lm[i-1, 3]+runif(1, -eps/10, eps/10)))
    if(runif(1)<f(u, v, w)/f(lm[i-1, 1], 
                             lm[i-1, 2],
                             lm[i-1, 3])) 
      lm[i, ] <- c(u, v, w)
    else
      lm[i, ] <- lm[i-1, ]
  }
  lm
}
lm <- onoff1.MH(eps=0.5)
showit <- function(x) {
  B <- nrow(x)
  par(mfrow=c(2, 3))
  for(i in 1:3)
    plot(cumsum(x[, i])/1:B, 
       type="l", xlab="", ylab="")
  for(i in 1:3)
  hist(x[(B/2):B, i], 50, 
       freq = FALSE, main="", xlab="")

  round(quantile(x[(B/2):B, 2], c(0.025, 0.75)), 3)
}
showit(lm)

##   2.5%    75% 
##  5.340 26.162

onoff2.MH <- function(B=2*1e5, eps, 
            x=c(27, 30), y=c(10, 13), 
            z=c(0.75, 0.6), p=c(1/10, 1/10) ) {
  n <- length(x)
  f <- function(m, l, e) 
    prod(dpois(x, l+e*m)*dpois(y, l)*dnorm(e, z, p*z))
  lm <- matrix(0, B, 2*n+1)
  lm[1, ] <- c(mean((x-y)/z), y, z)
  u <- rep(0, n)
  v <- u
  for(i in 2:B) {
    w <- lm[i-1, 1]+runif(1, -eps, eps)
    w <- ifelse(w<0, 0, w)
    u <- lm[i-1, 1+1:n]+runif(n, -eps, eps)
    u <- ifelse(u<0, 0, u)
    v <- lm[i-1, n+1+1:n]+runif(n, -eps/10, eps/10)
    v <- ifelse(v<0, 0, v)
    v <- ifelse(v>1, 1, v)
    if(runif(1)<f(w, u, v)/
                f(lm[i-1, 1], lm[i-1, 1+1:n], 
                  lm[i-1, n+1+1:n])) 
      lm[i, ] <- c(w, u, v)
    else
      lm[i, ] <- lm[i-1, ]
  }
  lm
}
lm <- onoff2.MH(eps=0.5)
lm[1, ]

## [1] 25.50 10.00 13.00  0.75  0.60

lm[1e4, ]

## [1] 20.0430581 15.7947019 11.2341037  0.7253795  0.5928673

lm[1e5, ]

## [1] 23.1893881 11.4761880 15.5305271  0.7082571  0.6046551

B <- nrow(lm)
par(mfrow=c(2, 3))
for(i in 1:5)
  plot(cumsum(lm[, i])/1:B, 
       type="l", xlab="", ylab="")
par(mfrow=c(1, 1))

hist(lm[(B/2):B, 1], 50, 
       freq = FALSE, main="", xlab="")

round(quantile(lm[(B/2):B, 1], c(0.025, 0.975)), 1)

##  2.5% 97.5% 
##  13.8  35.7