Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 76
Текст из файла (страница 76)
Oncethe initialization of the direction numbers iv has been done, the jth point can be obtaineddirectly by XORing together those direction numbers corresponding to nonzero bits in theGray code of j, as described above.The Latin HypercubeWe might here give passing mention the unrelated technique of Latin square orLatin hypercube sampling, which is useful when you must sample an N -dimensionalspace exceedingly sparsely, at M points. For example, you may want to test thecrashworthiness of cars as a simultaneous function of 4 different design parameters,but with a budget of only three expendable cars. (The issue is not whether this is agood plan — it isn’t — but rather how to make the best of the situation!)The idea is to partition each design parameter (dimension) into M segments, sothat the whole space is partitioned into M N cells.
(You can choose the segments ineach dimension to be equal or unequal, according to taste.) With 4 parameters and 3cars, for example, you end up with 3 × 3 × 3 × 3 = 81 cells.Next, choose M cells to contain the sample points by the following algorithm:Randomly choose one of the M N cells for the first point. Now eliminate all cellsthat agree with this point on any of its parameters (that is, cross out all cells in thesame row, column, etc.), leaving (M − 1)N candidates. Randomly choose one ofthese, eliminate new rows and columns, and continue the process until there is onlyone cell left, which then contains the final sample point.The result of this construction is that each design parameter will have beentested in every one of its subranges.
If the response of the system under test isdominated by one of the design parameters, that parameter will be found withthis sampling technique. On the other hand, if there is an important interactionamong different design parameters, then the Latin hypercube gives no particularadvantage. Use with care.CITED REFERENCES AND FURTHER READING:Halton, J.H. 1960, Numerische Mathematik, vol.
2, pp. 84–90. [1]Bratley P., and Fox, B.L. 1988, ACM Transactions on Mathematical Software, vol. 14, pp. 88–100. [2]Lambert, J.P. 1988, in Numerical Mathematics – Singapore 1988, ISNM vol. 86, R.P. Agarwal,Y.M. Chow, and S.J. Wilson, eds. (Basel: Birkhaüser), pp. 273–284.Niederreiter, H. 1988, in Numerical Integration III, ISNM vol. 85, H. Brass and G. Hämmerlin,eds. (Basel: Birkhaüser), pp. 157–171.Sobol’, I.M.
1967, USSR Computational Mathematics and Mathematical Physics, vol. 7, no. 4,pp. 86–112. [3]Antonov, I.A., and Saleev, V.M 1979, USSR Computational Mathematics and MathematicalPhysics, vol. 19, no. 1, pp. 252–256. [4]Dunn, O.J., and Clark, V.A. 1974, Applied Statistics: Analysis of Variance and Regression (NewYork, Wiley) [discusses Latin Square].316Chapter 7.Random Numbers7.8 Adaptive and Recursive Monte CarloMethodsThis section discusses more advanced techniques of Monte Carlo integration. Asexamples of the use of these techniques, we include two rather different, fairly sophisticated,multidimensional Monte Carlo codes: vegas [1,2] , and miser [3]. The techniques that wediscuss all fall under the general rubric of reduction of variance (§7.6), but are otherwisequite distinct.Importance SamplingThe use of importance sampling was already implicit in equations (7.6.6) and (7.6.7).We now return to it in a slightly more formal way.
Suppose that an integrand f can be writtenas the product of a function h that is almost constant times another, positive, function g. Thenits integral over a multidimensional volume V is&&&f dV = (f /g) gdV = h gdV(7.8.1)In equation (7.6.7) we interpreted equation (7.8.1) as suggesting a change of variable toG, the indefinite integral of g. That made gdV a perfect differential.
We then proceededto use the basic theorem of Monte Carlo integration, equation (7.6.1). A more generalinterpretation of equation (7.8.1) is that we can integrate f by instead sampling h — not,however, with uniform probability density dV , but rather with nonuniform density gdV . Inthis second interpretation, the first interpretation follows as the special case, where the meansof generating the nonuniform sampling of gdV is via the transformation method, using theindefinite integral G (see §7.2).More directly, one can go back and generalize the basic theorem (7.6.1) to the caseof nonuniform sampling: Suppose that points xi are chosen within the volume V with aprobability density p satisfying&p dV = 1(7.8.2)The generalized fundamental theorem is that the integral of any function f is estimated, usingN sample points xi , .
. . , xN , by*3 4&&ff 2 /p2 − f /p2f±pdV ≈(7.8.3)I ≡ f dV =ppNwhere angle brackets denote arithmetic means over the N points, exactly as in equation(7.6.2). As in equation (7.6.1), the “plus-or-minus” term is a one standard deviation errorestimate. Notice that equation (7.6.1) is in fact the special case of equation (7.8.3), withp = constant = 1/V .What is the best choice for the sampling density p? Intuitively, we have alreadyseen that the idea is to make h = f /p as close to constant as possible. We can be morerigorous by focusing on the numerator inside the square root in equation (7.8.3), which isthe variance per sample point.
Both angle brackets are themselves Monte Carlo estimatorsof integrals, so we can write3 2 4 3 42 & 2&2 & 2&2fffffS≡≈pdV−=(7.8.4)pdVdV−fdV−p2pp2ppWe now find the optimal p subject to the constraint equation (7.8.2) by the functional variation&&2&δf2(7.8.5)dV −f dV + λ p dV0=δpp7.8 Adaptive and Recursive Monte Carlo Methods317with λ a Lagrange multiplier.
Note that the middle term does not depend on p. The variation(which comes inside the integrals) gives 0 = −f 2 /p2 + λ or|f ||f |p= √ = '|f| dVλ(7.8.6)where λ has been chosen to enforce the constraint (7.8.2).If f has one sign in the region of integration, then we get the obvious result that theoptimal choice of p — if one can figure out a practical way of effecting the sampling — isthat it be proportional to |f |. Then the variance is reduced to zero. Not so obvious, but seento be true, is the fact that p ∝ |f | is optimal even if f takes on both signs. In that case thevariance per sample point (from equations 7.8.4 and 7.8.6) is&2 &2S = Soptimal =|f | dV−f dV(7.8.7)One curiosity is that one can add a constant to the integrand to make it all of one sign,since this changes the integral by a known amount, constant × V .
Then, the optimal choiceof p always gives zero variance, that is, a perfectly accurate integral! The resolution ofthis seeming paradox (already mentioned at the end' of §7.6) is that perfect knowledge of pin equation (7.8.6) requires perfect knowledge of |f |dV , which is tantamount to alreadyknowing the integral you are trying to compute!If your function f takes on a known constant value in most of the volume V , it iscertainly a good idea to add a constant so as to make that value zero. Having done that, theaccuracy attainable by importance sampling depends in practice not on how small equation(7.8.7) is, but rather on how small is equation (7.8.4) for an implementable p, likely only acrude approximation to the ideal.Stratified SamplingThe idea of stratified sampling is quite different from importance sampling.
Let usexpand our notation slightly and let f denote the true average of the function f overthe volume V (namely the integral divided by V ), while f denotes as before the simplest(uniformly sampled) Monte Carlo estimator of that average:&11 f ≡f dVf ≡f (xi)(7.8.8)VN iThe variance of the estimator, Var (f ), which measures the square of the error of theMonte Carlo integration, is asymptotically related to the variance of the function, Var (f ) ≡f 2 − f 2, by the relationVar (f ) =Var (f )N(7.8.9)(compare equation 7.6.1).Suppose we divide the volume V into two equal, disjoint subvolumes, denoted a and b,and sample N/2 points in each subvolume. Then another estimator for f , different fromequation (7.8.8), which we denote f , is1,(7.8.10)f ≡f a + f b2in other words, the mean of the sample averages in the two half-regions.
The variance ofestimator (7.8.10) is given by,- 1,-,Var f =Var f a + Var f b4Varb (f )1 Vara (f )(7.8.11)+=4N/2N/21=[Vara (f ) + Varb (f )]2N318Chapter 7.Random NumbersHere Vara (f ) denotes the variance of f in subregion a, that is, f 2a − f 2a , andcorrespondingly for b.From the definitions already given, it is not difficult to prove the relation11(7.8.12)[Vara (f ) + Varb (f )] + (f a − f b)224(In physics, this formula for combining second moments is the “parallel axis theorem.”)Comparing equations (7.8.9), (7.8.11), and (7.8.12), one sees that the stratified (into twosubvolumes) sampling gives a variance that is never larger than the simple Monte Carlo case— and smaller whenever the means of the stratified samples, f a and f b, are different.We have not yet exploited the possibility of sampling the two subvolumes with differentnumbers of points, say Na in subregion a and Nb ≡ N − Na in subregion b.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
