11

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?

MathematicalOrchid
  • 2,430
  • 3
  • 13
  • 15
  • 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
  • 2
    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
  • 1
    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
  • 1
    @Macro That procedure is unrewarding because that's the *worst* you can do. As jbowman points out, very mild assumptions about $f$ can lead to far better estimates. BTW, it is meaningless to stipulate that $f$ is "finite." If it's a well-defined function, all its values are real numbers and *a fortiori* finite. If you meant "bounded," that does you no good unless you know the bounds beforehand. – whuber Jun 26 '12 at 15:20
  • @whuber, yes I meant bounded. I agree it is a ridiculously inefficient algorithm but I don't think it's required to know the bounds beforehand - if you had forever to plug in random uniforms, I don't see how the CLT wouldn't apply there. I guess one problem could be having single points of discontinuity in $f$ where it jumps to an extreme value, which would "break" any ability to quantify the error. – Macro Jun 26 '12 at 15:27
  • 1
    @Macro "Most" functions aren't continuous anywhere! In fact, I do not see how the CLT could possibly apply in general. $f$ could be the inverse CDF of literally any distribution, for instance, in which case your Monte-Carlo draws are sampling from that distribution--for which the CLT need not apply even if the integral itself (i.e., the mean) exists. I think it would be much more fruitful for the OP to narrow the question and respondents to follow jbowman's suggestions. – whuber Jun 26 '12 at 18:08
  • My guess is that your function is probably continuous or differentiable in a lot of places, except at your scene boundaries. Can you identify piecewise differentiable chunks, use quadrature there, and sum the chunks? – Patrick Caldon Jun 27 '12 at 06:53
  • And just to add - this problem could be pretty nasty. "Shaders" are effectively little computational routines. If your function is $f(x) = 0$ if an abstract computer (ie. shader) using a coding of $x$ as its program halts, and $1$ otherwise ... well I would not suggest Monte-Carlo. – Patrick Caldon Jun 27 '12 at 06:55

2 Answers2

8

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.

A supplementary reading to this would of course be Niederreiter (1992) monograph.

StasK
  • 29,235
  • 2
  • 80
  • 165
  • 3
    (+1) [Josef Dick](http://quasirandomideas.wordpress.com/) also has some interesting fairly recent related results. – cardinal Jun 26 '12 at 23:37
2

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

Michael R. Chernick
  • 39,640
  • 28
  • 74
  • 143
  • 3
    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
  • 1
    @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 R. 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
  • 2
    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 R. 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
  • Just wondering... When you evaluate $f(x)$, you then know _exactly_ how much of that vertical slice is under the curve. Can you not use that knowledge somehow to improve the estimate a bit? – MathematicalOrchid Jun 27 '12 at 19:59
  • I don't see how knowing the height of f at a single point x tells you anthing about the area. It does tell you that M must be > f(x). If you have f evaluated at several points you can approximate the integral by rectangles but who knows how good that approximation will be? – Michael R. Chernick Jun 27 '12 at 20:08