Monte Carlo Integration

1 AMonte Carlo IntegrationTHE techniques developed in this dissertation are all Monte Carlo methods. Monte Carlomethods are numerical techniques which rely on random sampling toapproximatetheirresults. Monte Carlo Integration applies this process to the numerical estimation of integrals. Inthis appendix we review the fundamental concepts of Monte Carlo Integration upon which ourmethods are based. From this discussion we will see why Monte Carlo methods are a particularlyattractive choice for the multidimensional Integration problems common in computer references for Monte Carlo Integration in the context of computer graphics include Pharrand Humphreys [2004], Dutr et al.

2 [2006], and Veach [1997].The term Monte Carlo methods originated at the Los Alamos National Laboratory inthe late 1940s during the development of the atomic bomb [Metropolis and Ulam, 1949]. Notsurprisingly, the development of these methods also corresponds with the invention of the firstelectronic computers, which greatly accelerated the computation of repetitive numerical [1987] provides a detailed account of the origins of the Monte Carlo Vegas algorithms are another class of method which rely on randomization to computetheir results. However, in contrast to Las Vegas algorithms, which always produce the exact orcorrect solution, the accuracy of Monte Carlo methods can only be analyzed from a statisticalviewpoint.

3 Because of this, we first review some basic principles from probability theory beforeformally describing Monte Carlo BackgroundIn order to define Monte Carlo Integration , we start by reviewing some basic ideas VariablesArandom variableXis a function that maps outcomes of a random process to numbers. Arandom variable can be either discrete ( , the roll of a six-sided die where a fixed set of outcomesis possible,X={1,2,3,4,5,6}), or continuous ( , a person s height, which can take on real valuesR). In computer graphics we more commonly deal with continuous random variables which takeon values over ranges of continuous domains ( , the real numbersRor the sphere of directions ).

4 Distributions and Density FunctionsThecumulative distribution function, or CDF, of a random variableXis the probabilitythat a value chosen from the variable s distribution is less than or equal to some thresoldx:c d f(x)=P r{X x}.( )The correspondingprobability density function, or PDF, is the derivative of the CDF:p d f(x)=dd xc d f(x).( )CDFs are always monotonically increasing, which means that the PDF is always important relationship arises from the above two equations, which allows us to compute theprobability that a random variable lies within an interval:P r{a X b}= bap d f(x)dx.( )151 From this expression it is clear that the PDF must always integrate to one over the full extent of Values and VarianceTheexpected valueorexpectationof a random variableY=f(X) over a domain (x) isdefined asE[Y]= (x)f(x)p d f(x)d (x),( )while itsvarianceis 2[Y]=E[(Y E[Y])2],( )where , thestandard deviation, is the square root of the variance.

5 From these definitions it iseasy to show that for any constanta,E[a Y]=a E[Y],( ) 2[a Y]=a2 2[Y].( )Furthermore, the expected value of a sum of random variablesYiis the sum of their expectedvalues:E[ iYi]= iE[Yi].( )From these properties it is possible to derive a simpler expression for the variance: 2[Y]=E[Y2] E[Y]2.( )152 Additionally, if the random variables areuncorrelated, a summation property also holds for thevariance1: 2[ iYi]= i 2[Yi].( ) Monte Carlo EstimatorThe Basic Carlo Integration uses random sampling of a function to numeri-cally compute an estimate of its integral. Suppose that we want to integrate the one-dimensionalfunctionf(x) fromatob:F= baf(x)dx.

6 ( )We can approximate this integral by averaging samples of the functionfat uniform random pointswithin the interval. Given a set ofNuniform random variablesXi [a,b) with a correspondingPDF of 1/(b a), the Monte Carlo estimator for computingFis FN =(b a)1N 1N i=0f(Xi),( )The random variableXi [a,b) can be constructed by warping a canonical random numberuniformly distributed between zero and one, i [0, 1):Xi=a+ i(b a).( )Using this construction, we can expand the estimator to: FN =(b a)1NN 1 i=0f(a+ i(b a)).( )1 This property is often made with the stronger condition that the variables areindependent; however, it suffices forthem to be FN is a function ofXi, it is itself a random variable.]]]

7 We use this notation to clarify that FN is an approximation , the Monte Carlo estimator in Equation computes the mean value of thefunctionf(x) over the intervalatob, and then multiplies this mean by the length of the interval(b a). By moving (b a) into the summation, the estimator can be thought of as choosing aheight at a random evaluation of the function and averaging a set of rectangular areas computedby multiplying this height by the interval length (b a). These two interpretations are illustratedin Figure Value and ConvergenceIt is easy to show that the expected value of FN is in factF:E[ FN ]=E[(b a)1NN 1 i=0f(Xi)],=(b a)1NN 1 i=0E[f(Xi)],from Equations and (b a)1NN 1 i=0 baf(x)p d f(x)dx,from Equation 1 i=0 baf(x)dx,sincep d f(x)=1/(b a)= baf(x)dx,=F.

8 ( )Furthermore, as we increase the number of samplesN, the estimator FN becomes a closer andcloser approximation ofF. Due to theStrong Law of Large Numbers, in the limit we can guaranteethat we have the exact solution:P r{limN FN =F}=1.( )In practice we are interested in knowing just how quickly this estimate converges to a154 Figure : An illustration of the two interpretations of the basic Monte Carlo estimator in Equation four samples: computing the mean value, or height, of the function and multiplying by the intervallength (top), or computing the average of several rectangular areas (bottom).sufficiently accurate solution.

9 This can be analyzed by determining theconvergence rateof theestimators variance. In the next section we will show that the standard deviation is proportionalto: [ FN ] 1pN.( )Unfortunately, this means that we must quadruple the number of samples in order to reduce theerror by half!Standard Integration techniques exist which converge much faster in one dimension;however, these techniques suffer from thecurse of dimensionality, where the convergence ratebecomesexponentiallyworse with increased dimensions. The basic Monte Carlo estimator abovecan easily be extended to multiple dimensions, and, in contrast to deterministic quadraturetechniques, the convergence rate for Monte Carlo is independent of the number of dimensionsin the integral.

10 This makes Monte Carlo Integration the only practical technique for many highdimensional Integration problems, such as those encountered when computing global IntegrationMonte Carlo Integration can be generalized to use random variables drawn from arbitraryPDFs and to compute multidimensional integrals, such asF= (x)f(x)d (x),( )with the following modification to Equation : FN =1NN 1 i=0f(Xi)p d f(Xi).( )It is similarly easy to show that this generalized estimator also has the correct expected value:E[ FN ]=E[1NN 1 i=0f(Xi)p d f(Xi)],=1NN 1 i=0E[f(Xi)p d f(Xi)],=1NN 1 i=0 f(x)p d f(x)p d f(x)dx,=1NN 1 i=0 f(x)dx,= f(x)dx,=F.