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

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

• what's wrong with nearly zero value? how near to zero it is? May 14, 2019 at 15:24
• @Aksakal I have not mentioned, but it is actually smaller than the reported error. (e.g. 1.5*1e-17 with error 1e-16) May 15, 2019 at 0:45

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)