Calculation of the expectation of a posterior distribution using numerical integration methods

Question

I want to calculate the expectation of the following posterior distribution:

$$E( \theta \mid {\bf u} ) = \int\limits_{ - \infty }^\infty \theta \cdot g(\theta \mid {\bf u} )\,d\theta $$

and if possible the variance of it too:

$$\mathop{\rm var} ( \theta \mid {\bf u} ) = \int\limits_{ - \infty }^\infty \theta ^2 \cdot g( \theta \mid {\bf u} ) \, d\theta - \left[ E( \theta \mid {\bf u} ) \right]^2$$

${\bf u}$ is a vector of 1's and 0's, which are known. And $\theta \in (-\infty, \infty).$

$g( \theta \mid {\bf u})$ is a usual posterior distribution defined as:

$$g(\theta \mid {\bf u}) = \frac{L( \theta \mid {\bf u})g(\theta)}{\int\limits_{-\infty}^\infty L( \theta \mid {\bf u})g(\theta) \, d\theta }$$

where $L( \theta \mid {\bf u} )$ is the likelihood defined as:

$$L( \theta \mid {\bf u} ) = \prod\limits_{i = 1}^n P_i^{u_i}{( 1 - P_i )^{1 - u_i}}$$

and $P_i$ is defined as:

$$P_i = P_i( u_i = 1\mid \theta ) = {\beta_{3i}} + \frac{\beta_{4i} - \beta_{3i} }{{1 + {e^{ -\beta_{1i}( \theta - \beta_{2i} )}}}}$$

In this equation all $\beta$ values can be treated as known constants. $g(\theta)$ is the prior distribution of $\theta$. It can be any continuous distribution (such as normal, beta, uniform, etc.) but generally a normal distribution is used:

$$g(\theta) = \frac{1}{\sigma \sqrt {2\pi } }{e^{ - \frac{( \theta - \mu)^2}{2\sigma ^2}}}$$

where $\mu$ and $\sigma$ are also known.

What I want is to calculate this expectation and variance fast. These computations will be used in a long simulation study. So I eliminate methods like Monte-Carlo integration. I also want them to be precise to a given level (such as 0.001, or even smaller). In couple articles that I've read they calculated these integrals using Gauss-Hermite quadrature. But they did not tell the specifics of the calculation, how they integrate different prior distributions, etc. My limited understanding of Gauss-Hermite quadrature (HERE) tells me that I have to reparameterize these integrals to obey its form. But I cannot able to do that.

Any help will be much appreciated.

Note: I will be using R in simulations, but I don't want to use any package.

Answer 1

Gauss-Hermite is a good method for solving this problem. In your problem, the posterior mean can be written as: $$ E(\theta | \mathbf{u}) = \frac{\int \theta L(\theta | \mathbf{u}) g(\theta) d\theta}{\int L(\theta | \mathbf{u}) g(\theta) d\theta} $$ This is a ratio of two integrals. You need to apply Gauss-Hermite to each integral separately. Let's start with the denominator since it is the simplest. The denominator integral is equivalent to $E[L(\theta | \mathbf{u})]$. The Wikipedia page that you linked to explains how to approximate any such expectation. Let's look at the last formula on that page: $$ E[h(y)] \approx \frac{1}{\sqrt{\pi}} \sum_{i=1}^n w_i h(\sqrt{2} \sigma x_i + \mu) $$ In your problem, $y$ is $\theta$, $h(y)$ is $L(\theta | \mathbf{u})$, $\mu$ is the prior mean of $\theta$, and $\sigma$ is the prior standard deviation of $\theta$. Plug in the values of $x_i$ and $w_i$ from a table and you have your approximation.

The numerator integral is $E[\theta L(\theta | \mathbf{u})]$, so just change $h(y)$ to be $\theta L(\theta | \mathbf{u})$ and apply the formula again. For the variance, you also need the second moment of $\theta$ whose numerator is $E[\theta^2 L(\theta | \mathbf{u})]$ and the denominator is the same as above. So in total you apply Gauss-Hermite approximation three times to get the posterior mean and variance of $\theta$.