Problem
Say I have the following function $g(x)$, which is proportional to the density function $f_\theta(\theta)$ of random variable $\theta$, i.e. $g(\theta) \propto f(\theta)$, such that
$$ \begin{aligned} f(\theta|\mathbf{y}) \propto g(\theta) &= \prod_{i=1}^{100} \left( \frac{\exp(-\theta) * \theta^{y_i}}{y_i!} \right) \times \frac{1}{\sqrt{2\pi}}\exp \left( - \frac{1}{2} \theta^2\right) \\[8pt] &= \exp \left( \sum_{i=1}^{100} \left[ - \theta + y_i\log(\theta) - \log(y_i!) \right] -\log(\sqrt{2\pi}) - \theta^2/2 \right) \end{aligned} $$
for $y_1 = 23, y_2 = 28, \cdots y_{100} = 13$ (no patterns, just data between $(10, 30)$.
Find $\mathbb{E}(\theta | \mathbb{y}), \mathrm{Var}(\theta | \mathbb{y})$ by numerical integration.
Try
Of course I could try MCMC to sample and estimate $\mathbb{E}(\theta | \mathbb{y}), \mathrm{Var}(\theta | \mathbb{y})$, but I would like to find them by numerical integration.
In R
, we have function integrate()
, where "adaptive quadrature of functions of one variable over a finite or infinite interval" (quoted from ?integrate
).
Consider
$$\mathbb{E}(\theta | \mathbb{y}) = \int_{-\infty}^\infty f(\theta|\mathbf{y}) \theta d \theta = \frac{1}{\int_{-\infty}^\infty g(\theta) d\theta} \int_{-\infty}^\infty g(\theta) \theta d \theta$$
But if I type integrate(g, -Inf, Inf)
, it returns nearly zero value.
Alternatively, if I type integrate(psi, -Inf, Inf)
, where psi = function(theta) {g(theta) * theta}
, it returns again nearly zero value.
How should I proceed? Any help will be appreciated.