Journal of Undergraduate Research
Volume 7, Issue 6 - July/August 2006

Monte Carlo Calculations for Rapidly Changing Distributions

Arif H. Ahmed and John R. Klauder

ABSTRACT

Many dimensional integrals are frequently approximated numerically by using importance sampling in Monte Carlo calculations. However, certain distributions of interest involve rapidly changing behavior and conventional Monte Carlo methods are unsuitable in such cases. For this work, we show that a modified Monte Carlo procedure produces satisfactory results.

INTRODUCTION

The Monte Carlo method can be best described in terms of a relationship between geometric and algebraic aspects of theoretical science. Over the years, powerful methods expressing solutions of probabilistic problems in analytical terms have been developed. With the advancement of computational power, now the statistical problems can be modeled directly, allowing solving of any problem with an analytical expression resembling a probability problem by constructing a realization of the random variable problem.

Contemporary problems in physics often require working with systems with a large number of degrees of freedom. Among these are many atoms in a chunk of condensed matter, the many electrons in an atom, or the infinitely many values of a quantum field at all points in a region of space-time. These systems are often described by the evaluation of integrals of very high dimension. The straightforward evaluation of such integrals is almost impossible in conventional ways. The Monte Carlo method [1-3] provides ways of efficiently evaluating integrals of high dimension.

The Monte Carlo Method

Since it is not very clear how one can actually calculate an integral with random numbers, we can start out with a simple two dimensional problem. Consider a square with corners at (±1/2, ±1/2) and we want to calculate the area of a circle inscribed within it. If we select a number of uniformly distributed points in the square at random, we know that the probability that they will fall within the circle is equal to the ratio of the areas, which is π/4. The fraction of points that fall within the circle will represent the area of the circle.

Consider an integral of the form

(1)

By normalization, we can modify the problem such that

(2)

Now g(x) can be regarded as a probability distribution function, over which we can compute the average of the function f(x) (denoted as <f>).

(3) 

If we draw a large number, N, of random values of x from the distribution described by g(x), we can directly estimate the average value of f(x) by

(4)

This equation forms the basis of Monte Carlo evaluations of integrals, and the set of N xi are the realization of the probability distribution function, g(x).   

Even though computationally intensive, in the Monte Carlo method the computation time increase is only linear (sometimes quadratic) with additional dimensions, while for a classical approach, it is exponential.

The Metropolis Algorithm

Even though the basic Monte Carlo method itself improves the multidimensional integral problem we are trying to solve, it is not the most efficient technique to sample. The Metropolis algorithm represents the function in a new way, leading to a more efficient and accurate way of solving physical problems.

The Metropolis algorithm is based on the concept of a random walk. In the two dimensional classical random walk, a sample point starts at the origin of an x-y coordinate system, and with successive iteration, is allowed to move one unit in any coordinate direction in the plane with equal probability. In principle, with much iteration, this will cover the entire 2-dimensional space. But for an efficient approach, we want the walker to spend more time in regions where the function to be sampled is large, rather than covering the entire space. This will provide a better realization of the probability density function since the probability of finding a sample walker in a given region should be the same as g(x₁, x₂ … … x_α).  

Difficulties in evaluating rapidly changing probability density functions and possible solutions

One of the main principles of a Monte Carlo calculation for integrals is that we have a distribution function g(x) that represents the weight function. To describe how g(x) acts, one obvious question that arises is how we choose the step size of the walker. If the step size is large, then most of the trial steps will be rejected. If the step size is very small, then most trial steps will be accepted, but the random walker will never move too far, thus generating a poor sample of the function of interest.

Also, not all distribution functions are expected to have a smooth shape that can be easily captured by numerical simulation. Most real life problems involve discontinuities or sharp changes.

In our specific application, the function of interest [4] is:

(5)

which has the following probability distribution function:

(6)

where m=1, r=0, and our main goal was to observe the effect of ε in evaluating our function of interest.

To help visualize the rapid change in g, we can consider the plot of the function g(x) with ε = 0 and ε = 0.1.

We notice that for smaller ε values, the change is much sharper than for larger ε values, confirming our expectations.

In the modified Metropolis algorithm, random G(x) values were picked and corresponding x values were used to find the proper walk for the simulation. From the plots in figure 2 we can observe that much of the rapid changes occurring near the vertical axis were captured well, thus providing an accurate representation of the function of interest.

Calculation of average quantities for conditions of interest

For our particular problem, we considered three average quantities of interest. Measurements for 5 different ε values for 3 different dimensions were evaluated. For each walk associated with a particular dimension, 10,000 steps were taken to reach the best g-distributed value, and the following tables listing the average quantities were calculated.

CONCLUSION

From tables 1a, b and c, we can conclude that even though the calculation involved random numbers, the answers to the average quantities were very precise for the same parameters. This shows that we can depend on the calculation method for consistent results. We also note that for particular ε values, the average quantities for different dimensions have a consistent ratio representing a linear trend.

One of our main purposes of the study was to observe whether the Monte Carlo method can detect the rapidly changing characteristics of the distribution function. From the results obtained above, we notice that the different ε values resulted in different values for the average quantities. We can conclude that the algorithm does indeed take the rapid change of the probability density function into account.

Even though the modified Metropolis algorithm of the Monte Carlo method was successful in our case, we found it to be computationally intensive. For larger models, having thousands of dimensions, it might not be a suitable procedure with current resources available. Further investigation could be carried out to improve the computation time, for example, by taking fewer walks, but still obtaining a statistically significant answer.

---

REFERENCES

1. Metropolis, Nicholas, Rosenbluth, Marshall. “Equation of State Calculations by Fast Computing Machines.” The Journal of Chemical Physics. 21 (1953): 1087-1092.
2. Gibbs, William R. “Introduction to Monte Carlo.” Computation in Modern Physics. World Scientific (1994). 28-52.
3. Koonin, Steven E. “Monte Carlo Methods.” Computational Physics. Benjamin/Cummings (1986).185-215.
4. Klauder, John R. “Ultralocal Models: Functional integral formulation.” Beyond Conventional Quantization. Cambridge University Press (2000). 260-263.

---

--top--

Back to the Journal of Undergraduate Research