3
$\begingroup$

I am experiencing what seems to be a bias in importance sampling, which, given that it's an unbiased procedure, should not be there.

Consider linear regression $$ y = X\beta+\epsilon $$ where there are $n$ observations and $k$ regressors, and $\epsilon \sim N(0, 1)$. Let the prior for $\beta$ be independent $N(0,1)$, so that the posterior is $$ \beta|Y\sim N(V^{-1}X'y, V^{-1}), $$ where $V=I+X'X$. The marginal distribution of $y$ is $N(0, I+XX')$.

What I am doing is to estimate the marginal likelihood by importance sampling. Let $\bar{\beta}$ and $V_\beta$ be the mean and covariance based on $R$ draws from the posterior distribution, and $q$ be the $N(\bar{\beta}, V_\beta)$ distribution. I generate $L$ samples from $q$, and then evaluate $$ \frac{1}{L}\sum_{l=1}^L\frac{p(y|\beta^{(l)})p(\beta^{(l)})}{q(\beta^{(l)})}, $$ which by the usual properties of importance sampling is an unbiased estimator of $p(y)$. It seems to be so when $k$ is low (5, say). For higher dimensions, a bias appears to be introduced. But I fail to see why.

Here are three boxplots displaying the difference between the importance sampling estimates and the true marginal likelihood. They are generated using 20 replications, 1000 posterior and importance sampling draws and $(n, k)\in \{(20, 5), (200, 50), (2000, 500)\}$.

enter image description hereenter image description hereenter image description here

I am aware that off-the-shelf importance sampling may fail miserably for high-dimensional problems, but my understanding is that this is due to low efficiency. I am a little puzzled by this bias, and I'm not sure whether this is fully attributable to high variance. (I understand that the number of importance draws is low here, but that should be a problem of efficiency rather than bias.) If anyone can offer some insights or references for this, I would be very thankful. The R code is below.


library(mvtnorm)
k <- 5
n <- 20
reps <- 20
res <- matrix(NA, reps, 2)
set.seed(100)
for (i in 1:reps) {
  # Generate data
  beta_true <- rnorm(k)
  X <- matrix(rnorm(k*n), n, k)
  y <- X %*% beta_true + rnorm(n)

  # Posterior sampling
  beta_cov <- solve(diag(k) + crossprod(X))
  beta_mean <- beta_cov %*% crossprod(X, y)
  beta_post <- rmvnorm(1000, beta_mean, beta_cov)

  # Moments for importance sampling
  beta_m <- colMeans(beta_post)
  beta_c <- diag(diag(cov(beta_post)))
  beta_is <- rmvnorm(1000, beta_m, beta_c)

  is_likelihood <- sapply(1:nrow(beta_is), function(x) sum(dnorm(y, X %*% beta_is[x, ], 1, log = TRUE)))
  is_prior <- sapply(1:nrow(beta_is), function(x) sum(dnorm(beta_is[x, ], log = TRUE)))
  is_q <- dmvnorm(beta_is, beta_m, beta_c, log = TRUE)
  log_ratio <- is_likelihood + is_prior - is_q

  # Marginal likelihoods
  ml_is <- max(log_ratio) + log(mean(exp(log_ratio - max(log_ratio))))
  ml_true <- dmvnorm(matrix(y, nrow = 1), matrix(0, n, 1), diag(n) + tcrossprod(X), log = TRUE)
  res[i, ] <- c(ml_true, ml_is)
}
$\endgroup$
  • 1
    $\begingroup$ Have you checked the variance is finite? Otherwise it should indeed be unbiased. Or it may be that, since you compare the logs and not the marginals, the estimators get sufficiently biased to show the difference. $\endgroup$ – Xi'an Dec 14 '17 at 14:17
  • $\begingroup$ @Xi'an As far as I can tell, in the $k=500$, $n=2000$ it seems to behave nicely in the sense that if I do 100,000 draws the log ratio terms range between -4840 and -4721, so there are no outliers distorting the result. However, I think your point about the log scale might be the key here: by Jensen's inequality, I believe $E(\ln \hat{p}(y))\leq \ln[E\hat{p}(y)]=\ln p(y)$. Given that log marginal likelihoods are incredibly common, is this ever acknowledged (or corrected for)? Or am I missing something? $\endgroup$ – hejseb Dec 14 '17 at 14:50
  • $\begingroup$ The fact that the log transform does not preserve unbiasedness is well-known if this is what you ask. $\endgroup$ – Xi'an Dec 14 '17 at 15:05
  • $\begingroup$ @Xi'an Yes, but what I mean is that I rarely see this mentioned in the literature and instead the log-scale is used almost by default in many cases. But I'm guessing that there are no easy (and/or useful) ways to correct for it. My real problem that led to the question was that I needed a way to simulate the true value for benchmarking other methods against, for which I used importance sampling with 500,000 draws. But that was before I realized the apparent problems with the log transform... $\endgroup$ – hejseb Dec 14 '17 at 15:15
  • 1
    $\begingroup$ There are ways to turn a sequence of biased estimators into an unbiased one, check Russian Roulette as a key word. And thermodynamic integration. $\endgroup$ – Xi'an Dec 14 '17 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.