# How to estimate the accuracy of an integral?

An extremely common situation in computer graphics is that the colour of some pixel is equal to the integral of some real-valued function. Often the function is too complicated to solve analytically, so we're left with numerical approximation. But the function is also often very expensive to compute, so we're greatly constrained in how many samples we can compute. (E.g., you can't just decide to take one million samples and leave it at this.)

In general then, what you want to do is evaluate the function at randomly-chosen points until the estimated integral becomes "accurate enough". Which brings me to my actual question: How do you estimate the "accuracy" of the integral?

More specifically, we have $f : \mathbb{R} \rightarrow \mathbb{R}$, which is implemented by some complicated, slow computer algorithm. We want to estimate

$$k = \int_a^b f(x) \ dx$$

We can compute $f(x)$ for any $x$ we desire, but it's expensive. So we want to choose several $x$-values at random, and stop when the estimate for $k$ becomes acceptably accurate. To do this, of course, we need to know how accurate the current estimate actually is.

I'm not even sure what statistical tools would be appropriate for this kind of problem. But it appears to me that if we know absolutely nothing about $f$, then the problem is unsolvable. For example, if you compute $f(x)$ a thousand times and it's always zero, your estimated integral would be zero. But, knowing nothing about $f$, it's still possible that $f(x) = 1,000,000$ everywhere except the points you happened to sample, so your estimate is horribly wrong!

Perhaps, then, my question should have started with "what do we need to know about $f$ to make it possible to estimate the accuracy of our integral?" For example, often we know that it is impossible for $f$ to ever be negative, which would appear to be a highly relevant fact...

Edit: OK, so this seems to have generated lots of responses, which is good. Rather than reply to each of them individually, I'm going to try to fill in some additional background here.

When I say we know "nothing" about $f$, I mean that we can compute $f$, but we don't know anything further about it. I would expect (and the comments seem to agree) that having more knowledge allows us to use better algorithms. It seems that knowing bounds on $f$ and/or the first derivative of $f$ would be useful.

In most of the problems I'm thinking about, $f$ changes depending on scene geometry and the location within the scene under consideration. It's not some nice, tidy piece of algebra that you can analytically solve. Typically $f$ represents light intensity. Obviously light intensity can't ever be negative, but there is no limit to how large its positive values can be. And finally, object edges usually result in sharp discontinuities in $f$, and usually you can't predict where these are.

In short, $f$ is damned fiddly, so my first port of call was to ask what we can do with it given no further information. It appears that without at least some upper and lower bounds, the answer is "not a hell of a lot"... So it looks like I need to start making up some assumptions to make any headway here.

Also, given the number of times "Monte Carlo" has come up, I'm guessing that's the technical term for this sort of integration?

-
When you say "if we know absolutely nothing about $f$", what do you mean exactly? We can calculate $f$, right? –  Macro Jun 26 '12 at 14:34
To the close voters, it appears that this problem is about monte carlo integration based on the statement "So we want to choose several x-values at random, and stop when the estimate for k becomes acceptably accurate." which seems distinctly on topic since it involves random number generation/sampling, and parameter estimation so I'm not quite sure what the rationale for the "off topic" votes is. –  Macro Jun 26 '12 at 14:38
Typically, when you integrate over a known function, you can do much better than Monte Carlo integration. Monte Carlo converges to the true value at a rate of $1/\sqrt{N}$, where $N$ is the number of evaluation points. Other algorithms, e.g., quadrature-based, will converge at a rate $1/N$ or even faster (e.g., for a function that's periodic over the region of integration), assuming some level of smoothness of the function. Still others, based on quasi-random sequences (e.g., Sobol' sequences), will converge at intermediate rates, e.g., $(\ln N)^n/N$ for an $n$-dimensional integration. –  jbowman Jun 26 '12 at 14:48
This has clear but unrewarding answers. The answer to the second question is "nothing": the sole requirement is that $f$ be measurable, which is implicit in asking for its integral. But then the only thing you can do amounts to random sampling. With additional assumptions one can do much better at estimating the integral and assessing the accuracy. So a better question is "what improvements in accuracy estimation can be achieved with which assumptions." But this is overly broad. Therefore, please tell us what kind of functions you are currently dealing with. –  whuber Jun 26 '12 at 15:06
@whuber, with only the assumption that $f$ is finite on $(a,b)$ you can a) generate $X$ from a uniform(a,b) distribution, b) calculate $f(X)$ c) repeat $n$ times d) average the results. Then you have an estimate of the integral since it's proportional to $E(f(X))$ and an estimate of the variance ${\rm var}(f(X))/n$, and you can make statements about the error by invoking the central limit theorem. Why do you think that's not rewarding? –  Macro Jun 26 '12 at 15:16

For simplicity assume f(x)>=0 for all x in[a,b] and we know M such that f(x) < M for all x in [a,b]. The the integral I of f over [a,b] can be enclosed in the rectangle with width b-a and hight M. The integral of f is the proportion of the rectangle that falls under the function f multiplied by M(b-a). Now if you pick points in the rectangle at random and count the point as a success if it falls under the curve and as a failure otherwise you have set up a Bernoulli trial . The sample fraction of points inside is a binomial proportion and hence has mean p and variance p(1-p)/n where n is the number of points taken. Hence you can construct a confidence interval for p and since I =p M(b-a) a confidence interval for I also since for the estimate I^ =p^ M (b-a), Var(I^)= M$^2$(b-a)$^2$ p(1-p)/n. So to use statistics to determine the smallest n for which the integral is accurate enough you could specify an upper limit S on the variance of I^. Note p(1-p)/n <=1/(4n) for every 0<=p<=1. So set S=M$^2$(b-a)$^2$ /(4n) or n = smallest integer> M$^2$(b-a)$^2$ /(4S).

-
This would work under the assumptions you laid out in the first sentence but based on the problem description it seems unlikely that you can, a priori, bound the function values between $0$ and $M$. It appeared that all you're given is the ability to calculate $f$ and nothing else. –  Macro Jun 26 '12 at 14:58
@Macro Without knowing anything about f I do not see how one could say anything about the statistical accuracy of an estimate of the integral based on evaluating it at a fixed finite set of points. My assumptions are rather minimal. If f is bounded on the interval [a,b] there should be some M large enough that it could be used as an upper bound on f. –  Michael Chernick Jun 26 '12 at 15:05
I certainly agree with your first sentence, which begins to address the OP's second question. But, the method you've described requires that you, a priori, know $M$, which is not a particularly minimal assumption. –  Macro Jun 26 '12 at 15:12
It is an assumption. I used the term mimimal to say that i am making as few assumptions as possible to reach a definitive answer. –  Michael Chernick Jun 26 '12 at 15:23
What an ingenious idea... You're right, it doesn't work without bounds on $f$, but it looks like you can't do much of anything without that information. –  MathematicalOrchid Jun 26 '12 at 20:43
This is a non-trivial question that involves issues like total variation of $f$ and its reasonable multivariate extensions. The Stanford statistician Art Owen has worked on this using randomized quasi-Monte Carlo techniques. The regular Monte Carlo allows direct estimation of the accuracy of the integral, but each individual evaluation is not that accurate. Quasi-Monte Carlo produces more accurate estimates, but it is a fully deterministic technique, and as such does not allow to estimate the variance of your result. He showed how to combine the two approaches, and his paper is very lucid, so I won't try to reproduce it here.