1

I read this sentence in a book:

"... therefore this method is particularly useful for Bayesian inference since it doesn't require a normalization constant"

The method is a computational algorithm that takes a density (possibly un-normalized) and does some stuff with it.

I've seen this sentiment elsewhere too, that in Bayesian stats, we need methods that don't require normalization constants ... but I don't understand what the big deal is??

If we have a density $f$ on $R$ but we don't know the normalization constant ... give me two seconds to type this in R:

integrate(f, -Inf, Inf)

and out pops the (inverse of the) normalization constant.

So what is the big deal? Why do we care whether or not a method needs normalization constants when it takes 2 seconds to calculate them in R?

fewf
  • 11
  • The integration you do in R is still an approximation and carries some error. If you have a method that can avoid you this, why not? – polettix Sep 20 '19 at 15:58
  • Every computational method is an approximation dude... saying "hey u can't use that approximated normalization constant in that approximative algorithm" sounds a bit strange. I would agree with your statement if we wanted to use some exact, analytical formula or something, but that's not the case. Everything we do is an approximation. Besides, R's integrate is pretty darn accurate. – fewf Sep 20 '19 at 16:13
  • 3
    The problem is if there are many parameters/random variables in the density. High dimensional numerical integration is a very tough computational problem. Making it faster and more reliable in various settings remains a very active area of research. – CloseToC Sep 20 '19 at 16:26
  • 2
    Did you really just dismiss numerical integration of arbitrary functions in high dimensional spaces as trivial – Chris Haug Sep 20 '19 at 17:02

2 Answers2

4

The normalization constant involves an integration, and integration can be intractable in high dimensional problems.

In one dimension, you can evaluate $f$ at a set of points $X = \{x_1, x_2, \dots, x_n\}$ and use something like trapezoid rule to approximate the integral.

This set $X$ grows exponentially with increase in dimensions. eg. for two variables $x,y$, this becomes $\{(x_1,y_1), (x_1,y_2),\dots, (x_1,y_n), (x_2,y_1),\dots,(x_n,y_n)\}$.

See curse of dimensionality and numerical integration on wikipedia

fool
  • 2,440
2

To illustrate the problem, I decided to simulate the simplest problem that I encounter in statistical estimation. The sample is $\{-6,-1,1\}$. It is drawn from $$\frac{\sigma}{\pi[\sigma^2+(x-\mu)^2]}.$$ What is the posterior, and more importantly for me, what is the predictive distribution?

I used:

rm(list=ls())
options(digits = 12)

library(stats)
library(cubature)

f<-function(x){
  result<-1/(pi**3)*(x[1]**3)/(x[1]**2+(-1-x[2])**2)/(x[1]**2+(-6-x[2])**2)/(x[1]**2+(1-x[2])**2)
  return(result)
}

result<-cubintegrate(f,c(-Inf,-Inf),c(Inf,Inf),fDim = 2,method = "hcubature")
print(result)

There was no meaningful result because the error was too large relative to the answer. That is just three observations. It doesn't take much at all to break down multiple integrations over an arbitrary posterior, even of low dimension.

It is a big deal because it is difficult. It appears to look simple. It isn't. It is even worse when the posterior is located in a very narrow region, so small that you start having significant digit issues. Add just a little dimensionality or a lot of data and your built-in function is going to fall apart.

Finally, and this might be a small issue depending on your goals and your specific problem, but tools like integrate or cubintegrate can be biased estimators of the constant of integration.

Dave Harris
  • 7,630