How to calculate mean and variance from very small function proportional to the density?

Question

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.

Answer 1

This is exactly one of the reasons why people use Monte Carlo integration rather than deterministic numerical integration: you avoid the normalizing constant, which can be rather large or rather small.

Notice also that $\theta>0$, as it is the parameter of a Poisson distribution. You can use a half-normal distribution as a prior, but you need to take into consideration that this is a positive parameter.

If you have the extreme desire of using deterministic numerical integration, then you need to play with the integration range and other specifications. You are using the automatic choice -Inf, Inf, which works in simple problems (and actually it should be 0, Inf), but you cannot always trust automatic command specifications (computers are not as smart as people assume sometimes). For instance, check the following example:

# Simulated data
set.seed(123)
n <- 100
y = rexp(n)

# Unnormalized function
g <- Vectorize(function(theta) theta^(n*mean(y))*exp(-n*theta)*exp(-0.5*theta^2))
# Finding a good range
curve(g,0,2)

# Normalizing constant
norc <- integrate(g,0,3, subdivisions = 1000L)$value
norc 

# Normalized function
f <- Vectorize(function(theta) g(theta)/norc)
integrate(f,0,4)

# Mean
psi1 <- Vectorize(function(theta) theta*f(theta))

meanf <- integrate(psi1,0,5, subdivisions = 1000L)$value
meanf 

# Variance
psi2 <- Vectorize(function(theta) (theta-meanf)^2*f(theta))

varf <- integrate(psi2,0,5, subdivisions = 1000L)$value
varf

The idea is to plot the unnormalized function in order to identify an appropriate integration range, and hope that the integral is still within the reach of the precision of your computer. So, for instance, if $n$ is very large and the sample values are also large, you may end up with a normalizing constant of the order $10^{-100000}$, which cannot be calculated in this way unless you have access to NASA's supercomputing power.

Otherwise, MCMC will do a good job:

# Using Metropolis algorithm

library(mcmc)

log.g <- Vectorize( function(theta){
  if(theta>0) return( n*mean(y)*log(theta) -n*theta -0.5*theta^2 )
else return(-Inf)
})

out <- metrop(log.g, 1, 100000)

out$accept

# Histogram
hist(out$batch, probability = T, breaks = 100, ylim = c(0,4))
curve(f,0,1.5,add=T, col = "red", n=1000, lwd = 2)
# Mean
mean(out$batch)
# Variance
var(out$batch)