Software Engineering Online - Everything Up to Quiz #1

The expression "Monte Carlo method" is actually very general. Monte Carlo (MC) methods are stochastic techniques--meaning they are based on the use of random numbers and probability statistics to investigate problems. You can find MC methods used in everything from economics to nuclear physics to regulating the flow of traffic. Of course the way they are applied varies widely from field to field, and there are dozens of subsets of MC even within chemistry. But, strictly speaking, to call something a "Monte Carlo" experiment, all you need to do is use random numbers to examine some problem.
The use of MC methods to model physical problems allows us to examine more complex systems than we otherwise can. Solving equations which describe the interactions between two atoms is fairly simple; solving the same equations for hundreds or thousands of atoms is impossible. With MC methods, a large system can be sampled in a number of random configurations, and that data can be used to describe the system as a whole.
"Hit and miss" integration is the simplest type of MC method to understand, and it is the type of experiment used in this lab to determine the HCl/DCl energy level population distribution. Before discussing the lab, however, we will begin with a simple geometric MC experiment which calculates the value of pi based on a "hit and miss" integration.

The first figure is simply a unit circle circumscribed by a square. We could examine this problem in terms of the full circle and square, but it's easier to examine just one quadrant of the circle, as in the figure below.
If you are a very poor dart player, it is easy to imagine throwing darts randomly at Figure 2, and it should be apparent that of the total number of darts that hit within the square, the number of darts that hit the shaded part (circle quadrant) is proportional to the area of that part. In other words, If you remember your geometry, it's easy to show that.
If each dart thrown lands somewhere inside the square, the ratio of "hits" (in the shaded area) to "throws" will be one-fourth the value of pi. If you actually do this experiment, you'll soon realize that it takes a very large number of throws to get a decent value of pi...well over 1,000. To make things easy on ourselves, we can have computers generate random* numbers.
If we say our circle's radius is 1.0, for each throw we can generate two random numbers, an x and a y coordinate, which we can then use to calculate the distance from the origin (0,0) using the Pythagorean theorem. If the distance from the origin is less than or equal to 1.0, it is within the shaded area and counts as a hit. Do this thousands (or millions) of times, and you will wind up with an estimate of the value of pi. How good it is depends on how many iterations (throws) are done, and to a lesser extent on the quality of the random number generator. Simple computer code for a single iteration, or throw, might be: x=(random#) y=(random#) dist=sqrt(x^2 + y^2) if dist.from.origin (less.than.or.equal.to) 1.0 let hits=hits+1.0

The actual Monte Carlo method used in this lab to determine the population distribution among rotational energy levels is simpler than the two-dimensional example of the estimation of pi, as only one random number is generated for each "throw." This will be apparent shortly.
For this lab, the Boltzmann distribution can be solved analytically, and it is in fact used in determining the Monte Carlo distribution. As such, this is not a particularly informative simulation (you could just solve the Boltzmann equation for however many energy levels you wished and look at those numbers). However, this lab allows you to watch how changing the number of throws affects the results, and it automates the examination of the effects of temperature and isotope on population of energy levels.
The process used by the computer program for this lab is quite simple.

• 
• 
  
• The input information is used to solve the Boltzmann equation for some number of energy levels. For a maximum of J=4, the relative populations might look like the figure to the right.
• 
  
• The rest of the simulation is easier to understand if you imagine laying the peaks in the "spectrum" side-by-side, as shown below. If we generate random numbers along that line, whenever a number falls within the range of a particular J value, it counts as a "hit" for that energy level. Obviously, for shorter lengths (J=0) the number of "hits" will be smaller than for longer lengths (J=2).
  
  Since random number generators typically produce numbers on the range of zero to one, the population distribution is normalized so that the total "length" is equal to 1.0.
• 
  
• For each "throw," or random number produced, the computer determines which energy level range it belongs to, and calls it a "hit" for that J value. After the computer completes all throws, the number of hits (or relative number of hits) for each energy level are given. It is up to you to compare this output to the theoretical distribution produced by the Boltzmann equation, as described in the Lab Report Instructions.

* Computer-generated numbers aren't really random, since computers are deterministic. But, given a number to start with--generally called a random number seed--a number of mathematical operations can be performed on the seed so as to generate unrelated (pseudo-random) numbers. The output of random number generators is tested with rigorous statistical tests to ensure that the numbers are random in relation to one another. One caveat: If you use a random number seed more than once, you will get identical random numbers every time. Thus, for multiple trials, different random number seeds must be used. Commercial programs, like Mathematica, pull a random number seed from somewhere within the system--perhaps the time on the clock--so the seed is unlikely to be the same for two different experiments.

Monte Carlo methods provide approximate solutions to a variety of mathematical problems by performing statistical sampling experiments. They can be loosely defined as statistical simulation methods, where statistical simulation is defined in quite general terms to be any method that utilized sequences of random numbers to perform the simulation. Thus Monte Carlo methods are a collection of different methods that all basically perform the same process. This process involves performing many simulations using random numbers and probability to get an approximation of the answer to the problem. The defining characteristic of Monte Carlo methods is its use of random numbers in its simulations. In fact, these methods derive their collective name from the fact that Monte Carlo, the capital of Monaco, has many casinos and casino roulette wheels are a good example of a random number generator.

The Monte Carlo simulation technique has formally existed since the early 1940s, where it had applications in research into nuclear fusion. However, only with the increase in computer technology and power has the technique become more widely used. This is because computers are now able to perform millions of simulations much more efficiently and quickly than before. This is an important factor because it means that the technique can provide an approximate answer quickly and to a higher level of accuracy, because the more simulations that you perform then the more accurate the approximation is (This point is illustrated in the next section when we compare approximate error for different numbers of simulations).

Note that these methods only provide a approximation of the answer. Thus the analysis of the approximation error is a major factor to take into account when evaluating answers from these methods. The attempt to minimize this error is the reason there are so many different Monte Carlo methods. The various methods can have different levels of accuracy for their answers, although often this can depend on certain circumstances of the question and so some method's level of accuracy varies depending on the problem. This is illustrated well in the next section where we investigate four different Monte Carlo methods and compare there answers and the accuracy of their approximations.

One of the most important uses of Monte Carlo methods is in evaluating difficult integrals. This is especially true of multi-dimensional integrals which have few methods for computation and thus are suited to getting an approximation due to their complexity. It is in these situations that Monte Carlo approximations become a valuable tool to use, as it may be able to give a reasonable approximation in a much quicker time in comparison to other formal techniques.

In this section, we look at four different Monte Carlo methods to approach the problem of integral calculation. We investigate each method by giving a basic description of its individual procedure and then attempt to give a more formal mathematical description. We also discuss the example program that is included with this tutorial. It has implemented all four methods that we discuss, and also has results from these different methods in integrating an example function.

Before we discuss our first method it must be brought to the reader's attention the reasons why we are discussing four different Monte Carlo methods. The first point in explaining this is to acknowledge that there are many considerations when using Monte Carlo techniques to perform approximations. one of the chief concerns is to be able to get as accurate an approximation as possible. Thus, for each method we discuss their associated error statistics (the variance is discussed).

There is also the consideration of what technique is most suitable to the problem and so will get the best results. For example, a curve that has many platcaus (flat sections) may be very suitable to a stratified sampling method because the flat sections will have very low variance. Thus, we discuss four different Monte Carlo methods, because they all have their individual requirements and benefits. Probably the most important consideration, for the methods are the accuracy of their approximations because the more approximation, the less samples are needed to reach a specific level of accuracy.

The first method for discussion is crude Monte Carlo. We are going to use this technique to solve the integral 1. (see card #1 - Side A).

A basic description of this method is that you take a number, N, of random samples where a i=s (sample value) i=b. For each random sample s, we find the function value f(s) for the function f(x). We sum all of these values, and divide by N to get the mean value from our samples. We then multiply this value by the interval (b-a) to get the integral. This can be represented as: (see card #1 - Side B)

Now, the next part to describe is the accuracy of this approximation technique, because otherwise the answer will not be meaningful without a description of its uncertainty. In my implementation, I did a very simple, one-dimensional example. I found the variance of my sample by finding the mean of all of my N-group of samples. I then substituted the information into the equation to determine the variance. The sample variance equation is: (see card #2 - Side A).

Using this information you can determine the confidence intervals of your method and determine how accurate your answer is. This information about our implementation example, is discusses further on in the tutorial.

The next method for discussion is Acceptance - Rejection Monte Carlo. This is the easiest technique to explain and understand. However, it is also the technique with the least accurate approximation out of the four discussed.

A basic description of the Acceptance-Rejection method is that you have your integral as before. In the (a, b) interval for any given value of 'x' the function you find the upper limit. You then enclose this interval with a rectangle that is high enough to be above the upper limit, so we can be sure that the entire function for this interval is within the rectangle. We now begin taking random points within the rectangle. We now begin taking random points within the rectangle and evaluate this point to see if it is below the curve or not. If the random point is below the curve then it is treated as a successful sample. Thus, you take 'N' random points and perform this check remembering to keep count of the number of successful samples there have been. Now, once you have finished sampling, you can approximate the integral for the interval (a, b) by finding the area of the surrounding rectangle. You then multiply this area by the number of successful samples over the total number of samples, and this will give you an approximation of the integral for the interval (a, b). This is illustrated more clearly in the diagram of the integral for the interval (a, b). This is illustrated more clearly in the diagram.

Mathematically, we can say that the ratio of the area below the function f(x) and the whole area of the rectangle (Max(f(x)) * (b-a)) is approximately the ratio of the successful samples (k) and the whole number (N) of samples taken. Therefore: (see Card #2 - Side B).

To find the accuracy of the approximation, I have used the same variance technique as for the Crude Monte Carlo method. This analysis shows that the Acceptance-Rejection method gives a less accurate approximation than crude Monte Carlo.

The last method we look at is importance sampling. This is the most difficult technique to understand out of the four techniques. In fact, before any attempt to explain the basic principles behind this method, we will discuss a sort of link from the last technique (stratified sampling) to the concepts behind importance sampling.

Now in the figure we see that the four sub-intervals are quite different. I1 sees the value of f(x) staying very constant. however, as we progress across to B, we see the value of f(x) becoming larger than the fairly steady curve of I1. Now these larger values of f(x) are going to have more impact on the value of the integral. So should we not do more samples in the area where there is the highest values. By doing this, we will get a better approximation of a sub-interval that contributes more to the integral than the other sub-intervals. We won't skew the contributes more to the integral than the other sub-intervals. We won't skew the results either, because you still have to get the total integral by adding every sub-intervals together, all you have is a more accurate approximation of an important sub-interval of the curve.

This leads in to the method of importance of sampling. This method gets its name because it attempts to do more samples at the areas of the function that are more important. The way it does this is by bringing in a probability distribution function (pdf). All this is, is a function that attempts to say which areas of the function in the interval should get more samples. It does this by having a higher probability in that area. This diagram can be used as reference during my explanation. (see Card #3 - Side A & B).

Now first of all, note that we can define the integral equation as the following because they are equal - (see Card #4 - Side A).
p(x) is the probability distribution function. Note that the integral of p(x) over (a, b) is always equal to 1 and that for no value of 'x' within the interval (a, b) will p(x) evaluate to 0. The question is what do we do with p(x) now that we have put it into the equation. Well, we can use it so that when perform our samples of the curve f(x) within the interval (a,b), we can make the choices taking into account the probability of that particular shample getting selected. For example, accouding to the example graph about half of the pobability curve area is in the last quarter of the interval (a, b). Therefore, when we choose samples we should do it in a way so that half of the samples get taken in this area of the interval.

Now, you may be wondering why we should do this? and how we can do it without jeopardising the accuracy of our approximation? Well, the reason why we should do this, is that if we have chosen a good probability distribution function, then it should have a higher probability for samples to be selected at the important parts of the interval (the parts where the values are highest). Thus, we will spend more effort getting an accurate approximation of the important parts of the curve. But, this will affect the accuracy of our approximation because we will hopefully have a set of samples that focuses on certain parts of the curve. However, we counteract this by giving the value of f(x) / p(x) for every individual sample. This acts as a counterbalance to our unbalanced sampling technique. Thus the end result is a Monte Carlo method that effectively samples the important parts of the curve (as long as it a good probability distribution function) and then scales this sampling to give an approximation of the integral of f(x). Note again that the success of this method in getting a more accurate approximation is entirely dependent on selecting a good p(x). It has to be one that makes it more likely that a sample will be in an area of the interval (a, b) where the curve f(x) has a higher than the average value (and is thus more important to the approximation).

This method is effective in reducing error when it does have a good pdf because it samples the important parts of the curve more (due to the increased probability of a sample being selected in an important area). Thus it can get a good approximation of these important parts which lowers variance because these important parts are defined so because they have a larger effect on the overall approximation value.

These methods that are discussed previously are all important methods that do have some key differences. The last two methods are able to improve the accuracy of the approximations greatly, however they do need to have suitable conditions. In the case of importance sampling, it needs a good probability distribution function to come up with an effective approximation. In the case of stratified sampling it can come up with a much more accurate approximation if the shape of the curve is suitable and can be broken up into sections, with some being relatively flat thus allowing a very accurate approximation for that sub-interval. The first two methods, which are really the two basic Monte Carlo methods, are important to know as they both are used as the basis in more complex techniques.

The program implemented all four algorithms, and used them on the function illustrated in the following figure.

Our results did agree with our predictions. The acceptance-rejection method was the most inaccurate. Crude Monte Carlo was the next least effective approximation model. Then the importance sampling model was next with the Stratified Model being the most efficient model in my program. On reflection, this is a sensible result, because the Stratified Model splits the interval into four sub-intervals and two of these sub-intervals have a constant value of 1 (thus a variance of 0). Another contributing factor is in the fact that my probability distribution function models the function effectively enough but a better pdf would have resulted in more accurate results.

Note that a full print out of the results and the pogram code is in the Appendix. Here is the summary table of the variance results from the program. Note that in my implementation, I also created a hybrid method which was based on the one that I mentioned as a lead in to the discussion on Imporance Sampling. however, the results for this method were disappointing, although I think it was a fault in my implementation, I also created a hybrid method which was based on the one that I mentioned as a lead in to the discussion on Importance sampling. However, the results for this method were disappointing, although I think it was a fault in implementation because the method sounds feasible. I have included in the results table, although I am sure it should be able to approximate better than the standard stratified sampling model.

The previous section went into detail about the use of various Monte Carlo methods to evaluate integrals. However Monte Carlo techniques can be applied to many different forms of problems. In fact, Monte Carolo techniques are widely used in physics and chemistry to simulate complex reactions and interactions. This section is to illustrate the use of the basic Monte Carlo algorithm in another form of problem. It is a fairly simplistic example, however it illustrates Monte Carlo being used from a different perspective. This form of problem can also be seen in the focus questions in the appendix.

In this example, we have to imagine a coconut shy. We want to determine the probability that if we take 10 shots at the coconut shy we will have an even number of hits. The only fact that we know is that there is a 0.2 probability of having a hit with a single shot. We can work out the answer to this question using Monte Carlo by performing a large number of simulations of taking 10 shots at the coconut shy. We can then count all of the simulations that have an even number of hits and put that number over the toal number of simulations. This gives us an approximation of the probability of getting an even number of hits when we take 10 shys at the coconuts.

This example illustrates the use of the Monte Carlo algorithm for a different sort of problem. You may wonder where the use of random numbers is involved in this process. Well, as we know the probability of getting a hit with a single shot, we can use some random process to determine if each shot in each simulation is successful. For example, we could use a random number between 0 and 1 to simulate a shot at the coconuts. If the random number is between 0 and 0.2 then we can call it a hit, otherwise we can count it as a loss. (Note that the probability still stays at 0.2 to score a hit with a single shot). Thus now we can do all of the simulations just by generating random numbers and using these as our single shots at the coconuts.

Thus this is a simple illustration of using the principle behind Monte Carlo methods and applying it to a different form of problem to come out with an effective approximation of the answer. Note that this example is so simple that Monte Carlo techniques would not be a sensible chouice in this situation because the actual answer can be worked out with much less effort than performing a few hundred thousand simulations. However Monte Carlo techniques are more valuable when the problem is highly complex and where the effort required to get the actual answer is larger than the amount of effort required to get a reasonable approximation through simulation.

Two of the main reasons why we use Monte Carlo methods are because of their anti-aliasing properties and their ability to approximate quickly an answer that would be very time-consuming to find out the answer too if we were using methods to determine the exact answer.

This last point refers to the fact that Monte Carlo methods are used to simulate problems that are too difficult and time-consuming to use other methods for. An example is in the use of Monte Carlo tecniques in intergrating very complex multi-dimensional integrals. This is a task that other processes can not handle well, but which Monte Carlo can.

The first point refers to the fact that since Monte Carlo methods involve a random component in the algorithm, then this goes some way to avoiding the poblems of anti-aliasing (only for certain applications). An example that ws brought to my attention was that of finding the area of the black squares on a chess board. Now, if I was using an acceptance-rejection method to attack this problem, I should come out with a fiar approximation, due to the fact that I would be going to radom points on the chessboard. However what would happen if I was trying to do the same process but I used an algorithm that movged to a certain next point a set distance away and then continuted to do this, thus not having any random point selection. Well, the potential problem is that this person may have a bad step size and may overevaluate or underevaluate the number of successful trials he has, thus inevitably giving a poor approximation.

These are two solid reasons why people use Monte Carlo techniques. Other possible reasons could include its ease simulating complex physical systems in the fields of physics, engineering and chemistry.

Now from the above descriptions, we can see the value of Monte Carlo methods are their ability to give reasonable approximations for problems that can be very complex and time and resource consuming to solve. But how can this ability be used in the field of Computer Vision?

The area of computer vision that I first found Monte Carlo methods being mentioned was object tracking. The article that I found used a Monte Carlo technique discussed an algorithm called CONDENSATION - Conditional Density Propagation for Visual Tracking.

This technique was invented to handle the problem of tracking curves in dense visual clutter. Thus it tracks outlines and features of foreground objects, modeled as curves, as they move in substantial clutter, and in fact does this at fairly quick, efficient rate.

Basically, the algorithm has a cycle in which it is trying to predict the state that the object is in. It models the state an object is in and then as the picture moves ahead a frame, it tries to guess what likely states the object is now in. The big advantage with this technique is that it keeps multiple hypothesis of the object state open, thus allowing for a bit more flexibility with fixing incorrect guesses as to the object state.

Now, the computer holds these guesses as to what the state is going to be as probability distribution functions. In this case the areas of the curve with higher probabilty are the states the computer thinks are going to be next. The last step is where Monte Carlo Methods come into the procedure. This is where the computer readomly samples from this probability distribution function (which is guessing at the state of the object). It then checks these samples with image data to see if they support any of the states of the object at all. The more accurate the sample, then the more successfully weighted the sample response is.

Thus once you are finsihed taking the samples you now have the new probabilty distributions for the new frame of movement which will relfect the probable parameter state that the object is in.

Monte Carlo methods are very broad area of mathematics. They allow us to get reasonable approximations of very difficult problems through simulation and the use of random numbers. I discussed four of these methods that can be used in the evaluation of integrals. I then discussed my implementations of these methods and discussed my program results. Finally, I gave two examples of Monte Carlo methods being used within the field of Computer Science.

#include #include #include #include const int SAMPLES = 100000; /* /zeng/joy/mclab/mcintro.html */ int main(int argc, char *argv[]) { double h, x, y, z; int i; int countIn, countOut; countIn = 0; countOut = 0; /* initialize random seed: */ srand ( time(NULL) ); for (i =1; i

The Monte Carlo method can also be used to calculate values for definite integrals which either cannot be evaluated analytically, or are too complex to be calculated easily. Shown in the top figure to the right is a plot of the Gaussian function F(u)=2pi^1/2e^-u2. The integral, or area under the curve, of this function from zero to some number x you should recognize as the error function erf(x) of x. This is illustrated in the figure to the right for a value of x=1. As discussed in previous lessons values for the error function cannot be obtained analytically, and so one must resort to other methods to find the error function. One of these is the Monte Carlo method. To find the definite integral over the Gaussian in the figures above we first draw a box containing the function from the lower to the upper boundaries of the integral; for this example we will use x=0 to x=1. The is illustrated to the left of the figure below. If we pick random points within this box, the ratio of the number of points that fall below the curve to the total number of points chosen is equal to the ratio of the integral to the area of the box. Since we know the size of the sides of the box and hence its area, we can estimate the integral.
As in the calculation of pi, the accuracy of the value of the integral depends on the number of points that are used. With the speed of modern computers and the simplicity of this calculation, one can in practice use a very large number of points and obtain a relatively accurate integral. The answer ultimately depends, however, on how random the numbers actually are that the computer picks, as discussed in more detail below.