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$ Commented Oct 13, 2014 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
    Commented Oct 13, 2014 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$ Commented Oct 13, 2014 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
    Commented Oct 14, 2014 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
    Commented Oct 14, 2014 at 2:27

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

  • $\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
    Commented Oct 15, 2014 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
    Commented Oct 15, 2014 at 14:14

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.