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.

  • 3
    $\begingroup$ Are you asking for help w/ code? If so, this Q belongs on Stack Overflow. Add a reproducible example & we can migrate it for you. $\endgroup$ Oct 13 '14 at 13:17
  • $\begingroup$ Eventually I want a code, but I want to learn the logic of calculation of such integrals so that I can write the code myself. $\endgroup$
    – HBat
    Oct 13 '14 at 13:30
  • $\begingroup$ That's fine, it's just that 'how to get R to do this quickly' will end up being very R-code specific. $\endgroup$ Oct 13 '14 at 13:33
  • $\begingroup$ I'm not clear on your reasoning that because "these computations will be used in a long simulation study," it eliminates Monte Carlo integration. $\endgroup$
    – Sycorax
    Oct 14 '14 at 2:18
  • $\begingroup$ I need to calculate the first expectation and variance say 1 million times with different beta parameters, within a reasonable time. If I use a Monte-carlo integration which itself needs possibly 10k calculations for each single integration, the program will took very long. $\endgroup$
    – HBat
    Oct 14 '14 at 2:27

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$.

  • $\begingroup$ Dear @TomMinka, I can able to solve the problem using the explanation you offered. When I first looked at wikipedia page, I couldn't figured out that h(y) can be equal to the likelihood. I thought that I should write the full equation and try to get a form including e^(-x^2) from there, which I couldn't. Thank you very much for your clear explanation. I also want to change the prior distributions. As far as I understand, if I want to use a uniform distribution I need to use Gauss–Legendre similar to the way you showed, for beta I need to use Gauss–Jacobi, and etc. Right? $\endgroup$
    – HBat
    Oct 15 '14 at 10:54
  • $\begingroup$ Yes, you will need to change the quadrature rule for different priors. Not because the functional form of the prior is different, but mainly because the range of integration is different. $\endgroup$
    – Tom Minka
    Oct 15 '14 at 14:14

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.