Next: The Veto Algorithm Up: Monte Carlo Techniques Previous: Monte Carlo Techniques   Contents

## Selection From a Distribution

The situation that is probably most common is that we know a function which is non-negative in the allowed range . We want to select an at random' so that the probability in a small interval around a given is proportional to . Here might be a fragmentation function, a differential cross section, or any of a number of distributions.

One does not have to assume that the integral of is explicitly normalized to unity: by the Monte Carlo procedure of picking exactly one accepted value, normalization is implicit in the final result. Sometimes the integral of does carry a physics content of its own, as part of an overall weight factor we want to keep track of. Consider, for instance, the case when represents one or several phase-space variables and a differential cross section; here the integral has a meaning of total cross section for the process studied. The task of a Monte Carlo is then, on the one hand, to generate events one at a time, and, on the other hand, to estimate the total cross section. The discussion of this important example is deferred to section .

If it is possible to find a primitive function which has a known inverse , an can be found as follows (method 1):

 (2)

The statement of the first line is that a fraction of the total area under should be to the left of . However, seldom are functions of interest so nice that the method above works. It is therefore necessary to use more complicated schemes.

Special tricks can sometimes be found. Consider e.g. the generation of a Gaussian . This function is not integrable, but if we combine it with the same Gaussian distribution of a second variable , it is possible to transform to polar coordinates

 (3)

and now the and distributions may be easily generated and recombined to yield . At the same time we get a second number , which can also be used. For the generation of transverse momenta in fragmentation, this is very convenient, since in fact we want to assign two transverse degrees of freedom.

If the maximum of is known, in the range considered, a hit-or-miss method will always yield the correct answer (method 2):

14.
select an with even probability in the allowed range, i.e. ;
15.
compare a (new) with the ratio ; if , then reject the value and return to point 1 for a new try;
16.
otherwise the most recent value is retained as final answer.
The probability that is proportional to ; hence the correct distribution of retained values. The efficiency of this method, i.e. the average probability that an will be retained, is . The method is acceptable if this number is not too low, i.e. if does not fluctuate too wildly.

Very often does have narrow spikes, and it may not even be possible to define an . An example of the former phenomenon is a function with a singularity just outside the allowed region, an example of the latter an integrable singularity just at the and/or borders. Variable transformations may then be used to make a function smoother. Thus a function which blows up as for , with an close to 0, would instead be roughly constant if transformed to the variable .

The variable transformation strategy may be seen as a combination of methods 1 and 2, as follows. Assume the existence of a function , with over the range of interest. Here is picked to be a simple' function, such that the primitive function and its inverse are known. Then (method 3):

17.
select an according to the distribution , using method 1;
18.
compare a (new) with the ratio ; if , then reject the value and return to point 1 for a new try;
19.
otherwise the most recent value is retained as final answer.
This works, since the first step will select with a probability and the second retain this choice with probability . The total probability to pick a value is then just the product of the two, i.e. .

If has several spikes, method 3 may work for each spike separately, but it may not be possible to find a that covers all of them at the same time, and which still has an invertible primitive function. However, assume that we can find a function , such that over the range considered, and such that the functions each are non-negative and simple, in the sense that we can find primitive functions and their inverses. In that case (method 4):

20.
select an at random, with relative probability given by the integrals

21.
for the selected, use method 1 to find an , i.e.

22.
compare a (new) with the ratio ; if , then reject the value and return to point 1 for a new try;
23.
otherwise the most recent value is retained as final answer.
This is just a trivial extension of method 3, where steps 1 and 2 ensure that, on the average, each value picked there is distributed according to : the first step picks with relative probability , the second with absolute probability (this is one place where one must remember to do normalization correctly); the product of the two is therefore and the sum over all gives back .

We have now arrived at an approach that is sufficiently powerful for a large selection of problems. In general, for a function which is known to have sharp peaks in a few different places, the generic behaviour at each peak separately may be covered by one or a few simple functions , to which one adds a few more to cover the basic behaviour away from the peaks. By a suitable selection of the relative strengths of the different 's, it is possible to find a function that matches well the general behaviour of , and thus achieve a reasonable Monte Carlo efficiency.

The major additional complication is when is a multidimensional variable. Usually the problem is not so much itself, but rather that the phase-space boundaries may be very complicated. If the boundaries factorize it is possible to pick phase-space points restricted to the desired region. Otherwise the region may have to be inscribed in a hyper-rectangle, with points picked within the whole hyper-rectangle but only retained if they are inside the allowed region. This may lead to a significant loss in efficiency. Variable transformations may often make the allowed region easier to handle.

There are two main methods to handle several dimensions, each with its set of variations. The first method is based on a factorized ansatz, i.e. one attempts to find a function which is everywhere larger than , and which can be factorized into , where . Here each may in its turn be a sum of functions , as in method 4 above. First, each is selected independently, and afterwards the ratio is used to determine whether to retain the point.

The second method is useful if the boundaries of the allowed region can be written in a form where the maximum range of is known, the allowed range of only depends on , that of only on and , and so on until , whose range may depend on all the preceding variables. In that case it may be possible to find a function that can be integrated over through to yield a simple function of , according to which is selected. Having done that, is selected according to a distribution which now depends on , but with through integrated over. In particular, the allowed range for is known. The procedure is continued until is reached, where now the function depends on all the preceding values. In the end, the ratio is again used to determine whether to retain the point.

Next: The Veto Algorithm Up: Monte Carlo Techniques Previous: Monte Carlo Techniques   Contents
Stephen Mrenna 2005-07-11