I am trying to approximate the expectation of the "complete-data likelihood" with respect to the distribution of some missing data, and I am having some trouble. This expectation can be written as $$ \int p(\mathbf{y} \mid \mathbf{x}, \mathbf{\beta})f(\mathbf{x}) d\mathbf{x}_{1:n} = \int \prod_{i=1}^np(y_i \mid \mathbf{x}_i, \mathbf{\beta})f(\mathbf{x}_i) d\mathbf{x}_{1:n} $$ where

  1. $p(y_i \mid \mathbf{x}_i, \mathbf{\beta}) = \sigma(\mathbf{x}_i^T\beta)^{y_i} [1-\sigma(\mathbf{x}_i^T\beta)]^{1-y_i}$
  2. $f(\mathbf{x}_i)$ is chosen, so let's not make any assumptions about what it is yet. The $\sigma$ function is the inverse logit, or sigmoid function: $\sigma(r) = e^r/(1+e^r)$.

You can use importance sampling by picking an importance density $q(\mathbf{x}^s \mid \mathbf{y})$, simulating $N_s$ times from it, and calculating: $$ \sum_{s = 1}^{N_s} \frac{p(\mathbf{y} \mid \mathbf{x}^s, \mathbf{\beta})f(\mathbf{x}^s)}{q(\mathbf{x}^s \mid \mathbf{y})} \bigg/ N_s. $$

I am interested in reducing relative variance, and not getting all of these $-\infty$s. Are there well-known or popular strategies to accomplish this? Maybe I can simulate $\mathbf{x}^s$ so that all of these $\mathbf{x}_i^T\beta$s are constrained to be close to $0$.

Edit: perhaps I should be more clear about this: this is a logistic regression model with missing data, not a state space model. The phrase “complete data likelihood” could’ve been misleading as it means two different things in these two contexts. Also, it could’ve been misleading that I simplified notation for $f$, as it doesn’t have anything denoting missing and observed subvectors.


1 Answer 1


Importance sampling is notorious for falling afoul of dimensionality, the so called curse of dimension(ality). This means that a direct frontwise importance sampling attack in this latent variable setting cannot work.

The standard approach is to use instead sequential importance sampling, also called sequential Monte Carlo or particle filtering, where an importance sample is created for each $\mathbf{x}_i$ in the sequence $\mathbf{x}_{1:n}$, targeting the density proportional to $$p(y_i \mid \mathbf{x}_i, \mathbf{\beta})f(\mathbf{x}_i |\mathbf{x}_{1:{i-1}})$$and estimating its normalising constant from that sample.

  • 2
    $\begingroup$ Yes, thanks, (+1), I started in on this a few days ago as resampling might help some in this non-time-series context, but I'm still interested in the shape of this specific $p$. That's still an issue with these strategies, too. On the other hand, resampling might just add noise, because there’s nothing to learn from, going row to row. $\endgroup$
    – Taylor
    Commented Sep 20, 2018 at 5:36
  • 1
    $\begingroup$ Oh you're referring to the "data-tempering" thing. Sorry this took me a while to grasp, I'm usually in the state space model setting. Thanks again. $\endgroup$
    – Taylor
    Commented Nov 29, 2018 at 19:24

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.