Calculation of the expectation of a posterior distribution using numerical integration methods 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.
 A: 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$.
