Monte Carlo Simulation In Java 3x26p

Monte carlo Simulation in java This section under major construction. In 1953 Enrico Fermi, John Pasta, and Stanslaw Ulam created the first "computer experiment" to study a vibrarting atomic lattice. Nonlinear system couldn't be analyzed by classical mathematics. Simulation = analytic method that imitates a physical system. Monte Carlo simulation = use randomly generated values for uncertain variables. Named after famous casino in Monaco. At essentially each step in the evolution of the calculation, Repeat several times to generate range of possible scenarios, and average results. Widely applicable brute force solution. Computationally intensive, so use when other techniques fail. Typically, accuracy is proportional to square root of number of repetitions. Such techniques are widely applied in various domains including: deg nuclear reactors, predicting the evolution of stars, forecasting the stock market, etc. Generating random numbers. The math library function Math.random generate a pseudo-random number greater than or equal to 0.0 and less than 1.0. If you want to generate random integers or booleans, the best way is to use the library Random. Program RandomDemo.java illustrates how to use it. Random random = new Random(); boolean a = random.nextBoolean(); int b = random.nextInt(); int c = random.nextInt(100); double d = random.nextDouble(); double e = random.nextGaussian(); 1

// // // // //

true or false between -2^31 and 2^31 - 1 between 0 and 99 between 0.0 and 1.0 Gaussian with mean 0 and stddev =

Note that you should only create a new Random object once per program. You will not get more "random" results by creating more than one. For debugging, you may wish to produce the same sequence of pseudo-random number each time your program executes. To do this, invoke the constructor with a long argument. Random random = new Random(1234567L);

The pseudo-random number generator will use 1234567 as the seed. Use SecureRandom for cryptographically secure pseudo-random numbers, e.g., for cryptography or slot machines. Linear congruential random number generator. With integer types we must be cognizant of overflow. Consider a * b (mod m) as an example (either in context of a^b (mod m) or linear congruential random number generator: Given constants a, c, m, and a seed x[0], iterate: x = (a * x + c) mod m. Park and Miller suggest a = 16807, m

= 2147483647, c = 0 for 32-bit signed integers. To avoid overflow, use Schrage's method. Precompute: Iterate:

q = m / a, r = m % a x = a * (x - x/ q) * q) - r * (x / q)

Exercise: compute cycle length. Library of probability functions. OR-Objects contains many classic probability distributions and random number generators, including Normal, F, Chi Square, Gamma, Binomial, Poisson. You can the jar file here. Program ProbDemo.java illustrates how to use it. It generate one random value from the gamma distribution and 5 from the binomial distribution. Note that the method is called getRandomScaler and notgetRandomScalar. GammaDistribution x = new GammaDistribution(2, 3); System.out.println(x.getRandomScaler()); BinomialDistribution y = new BinomialDistribution(0.1, 100); System.out.println(y.getRandomVector(5));

Queuing models. M/M/1, etc. A manufacturing facility has M identical machines. Each machine fails after a time that is exponentially distributed with mean 1 / μ. A single repair person is responsible for maintaining all the machines, and the time to fix a machine is exponentially distributed with mean 1 / λ. Simulate the fraction of time in which no machines are operational. Diffusion-limited aggregation. Diffuse = undergo random walk. The physical process diffusion-limited aggregation (DLA) models the formation of an aggregate on a surface, including lichen growth, the generation of polymers out of solutions, carbon deposits on the walls of a cylinder of a Diesel engine, path of electric discharge, and urban settlement. The modeled aggregate forms when particles are released one at a time into a volume of space and, influenced by random thermal motion, they diffuse throughout the volume. There is a finite probability that the short-range attraction between particles will influence the motion. Two particles which come into with each other will stick together and form a larger unit. The probability of sticking increases as clusters of occupied sites form in the aggregate, stimulating further growth. Simulate this process in 2D using Monte Carlo methods: Create a 2D grid and introduce particles to the lattice through a launching zone one at a time. After a particle is launched, it wanders throughout with a random walk until it either sticks to the aggregate or wanders off the lattice into the kill zone. If a wandering particle

enters an empty site next to an occupied site, then the particle's current location automatically becomes part of the aggregate. Otherwise, the random walk continues. Repeat this process until the aggregate contains some pre-determined number of particles. Reference: Wong, Samuel, Computational Methods in Physics and Engineering, 1992. Program DLA.java simulates the growth of a DLA with the following properties. It uses the helper data type Picture.java. Set the initial aggregate to be the bottom row of the N-by-N lattice. Launch the particles from a random cell in top row. Assume that the particle goes up with probability 0.15, down with probability 0.35, and left or right with probability 1/4 each. Continue until the particles stick to a neighboring cell (above, below, left, right, or one of the four diagonals) or leaves the N-by-N lattice. The preferred downward direction is analogous to the effect of a temperature gradient on Brownian motion, or like how when a crystal is formed, the bottom of the aggregate is cooled more than the top; or like the influence of a gravitational force. For effect, we color the particles in the order they are released according to the rainbow from red to violet. Below are three simulations with N = 176; here is an image with N = 600.

Brownian motion. Brownian motion is a random process used to model a wide variety of physical phenomenon including the dispersion of ink flowing in water, and the behavior of atomic particles predicted by quantum physics. (more applications). Fundamental random process in the universe. It is the limit of a discrete random walk and the stochastic analog of the Gaussian distribution. It is now widely used in computational finance, economics, queuing theory, engineering, robotics, medical imaging, biology, and flexible manufacturing systems. First studied by a Scottish botanist Robert Brown in 1828 and analyzed mathematically by Albert Einstein in 1905. Jean-Baptiste Perrin performed experiments to confirm Einstein's predictions

and won a Nobel Prize for his work. An applet to illustrate physical process that may govern cause of Brownian motion. Simulating a Brownian motion. Since Brownian motion is a continuous and stochastic process, we can only hope to plot one path on a finite interval, sampled at a finite number of points. We can interpolate linearly between these points (i.e., connect the dots). For simplicitly, we'll assume the interval is from 0 to 1 and the sample points t0, t1, ..., tN are equally spaced in this interval. To simulate a standard Brownian motion, repeatedly generate independent Gaussian random variables with mean 0 and standard deviation sqrt(1/N). The value of the Brownian motion at time i is the sum of the first i increments.

Geometric Brownian motion. A variant of Brownian motion is widely used to model stock prices, and the Nobel-prize winning Black-Scholes model is centered on this stochastic process. A geometric Brownian motion with drift μ and volatility σ is a stochastic process that can model the price of a stock. The parameter μ models the percentage drift. If μ = 0.10, then we expect the stock to increase by 10% each year. The parameter σ models the percentage volatility. If σ = 0.20, then the standard deviation of the stock price over one year is roughly 20% of the current stock price. To simulate a geometric Brownian motion from time t = 0 to t = T, we follow the same procedure for standard Brownian motion, but multiply the increments, instead of adding them, and incorporate the drift and volatility parameters. Specifically, we multiply the current price by by (1 + μΔt + σsqrt(Δt)Z), where Z is a standard Gaussian and Δt = T/N Start with X(0) = 100, σ = 0.04. construction of BM. Black-Scholes formula. Move to here? Ising model. The motions of electrons around a nucleus produce a magnetic field associated with the atom. These atomic magnets act much like conventional magnets. Typically, the magnets point in random directions, and all of the forces cancel out leaving no overall magnetic field in a macroscopic clump of matter. However, in some materials (e.g., iron), the magnets can line up producing a

measurable magnetic field. A major achievement of 19th century physics was to describe and understand the equations governing atomic magnets. The probability that state S occurs is given by the Boltzmann probability density function P(S) = e E(S)/kT / Z, where Z is the normalizing constant (partition function) sum e-E(A)/kT over all states A, k is Boltzmann's constant, T is the absolute temperature (in degrees Kelvin), and E(S) is the energy of the system in state S. Ising model proposed to describe magnetism in crystalline materials. Also models other naturally occurring phenomena including: freezing and evaporation of liquids, protein folding, and behavior of glassy substances. Ising model. The Boltzmann probability function is an elegant model of magnetism. However, it is not practical to apply it for calculating the magnetic properties of a real iron magnet because any macroscopic chunk of iron contains an enormous number atoms and they interact in complicated ways. The Ising model is a simplified model for magnets that captures many of their important properties, including phase transitions at a critical temperature. (Above this temperature, no macroscopic magnetism, below it, systems exhibits magnetism. For example, iron loses its magnetization around 770 degrees Celsius. Remarkable thing is that transition is sudden.) reference First introduced by Lenz and Ising in the 1920s. In the Ising model, the iron magnet is divided into an N-by-N grid of cells. (Vertex = atom in crystal, edge = bond between adjacent atoms.) Each cell contains an abstract entity known as spin. The spin si of cell i is in one of two states: pointing up (+1) or pointing down (-1). The interactions between cells is limited to nearest neighbors. The total magnetism of the system M = sum of si. The total energy of the system E = sum of - J si sj, where the sum is taken over all nearest neighbors i and j. The constant J measures the strength of the spin-spin interactions (in units of energy, say ergs). [The model can be extended to allow interaction with an external magnetic field, in which case we add the term -B sum of sk over all sites k.] If J > 0, the energy is minimized when the spins are aligned (both +1 or both -1) - this modelsferromagnetism. if J < 0, the energy is minimized when the spins are oppositely aligned - this models antiferromagnetism. Given this model, a classic problem in statistical mechanics is to compute the expected magenetism. A state is the specification of the spin for each of the N^2 lattice cells. The expected magnetism of the system E[M] = sum of M(S) P(S) over all states S, where M(S) is the magnetism of state S, and P(S) is the probability of state S occurring according to the Boltzmann probability function. Unfortunately, this

equation is not amenable to a direct computational because the number of states S is 2N*N for an N-by-N lattice. Straightforward Monte Carlo integration won't work because random points will not contribute much to sum. Need selective sampling, ideally sample points proportional to e-E/kT. (In 1925, Ising solved the problem in one dimension - no phase transition. In a 1944 tour de force, Onsager solved the 2D Ising problem exactly. His solution showed that it has a phase transition. Not likely to be solved in 3D - see intractability section.) Metropolis algorithm. Widespread usage of Monte Carlo methods began with Metropolis algorithm for calculation of rigid-sphere system. Published in 1953 after dinner conversation between Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller. Widely used to study equilibrium properties of a system of atoms. Sample using Markov chain using Metropolis' rule: transition from A to B with probability 1 if Δ E <= 0, and with probability e-ΔE/kTif Δ E > 0. When applied to the Ising model, this Markov chain is ergodic (similar to Google PageRank requirement) so the theory underlying the Metropolis algorithm applies. Converges to stationary distribution. Program Cell.java, State.java, and Metropolis.java implements the Metropolis algorithm for a 2D lattice. Ising.java is a procedural programming version. "Doing physics by tossing dice." Simulate complicated physical system by a sequence of simple random steps. Measuring physical quantities. Measure magnetism, energy, specific heat when system has thermalized (the system has reached a thermal equilibrium with its surrounding environment at a common temperature T). Compute the average energy and the average magenetization over time. Also interesting to compute the variance of the energy or specific heat = <E2> - <E>2, and the variance of the magnetization or susceptibility<Χ> = <M2> - <M>2. Determining when system has thermalized is a challenging problem - in practice, many scientists use ad hoc methods. Phase transition. Phase transition occurs when temperature Tc is 2 / ln(1 + sqrt(2)) = 2.26918). Tc is known as the Curie temperature. Plot magnetization M (average of all spins) vs. temperature (kT = 1 to 4). Discontinuity of slope is signature of second order phase transition. Slope approaches infinity. Plot energy (average of all spinspin interactions) vs. temperature (kT = 1 to 4). Smooth curve through phase transition. Compare against exact solution. Critical temperature for which algorithm dramatically slows down. Below are the 5000th sample trajectory for J/kT = 0.4 (hot / disorder) and 0.47 (cold / order). The system becomes magnetic as temperature

decreases; moreover, as temperature decreases the probability that neighboring sites have the same spin increasing (more clumping). Experiments.  Start will above critical temperature. State converges to nearly uniform regardless of initial state (all up, all down, random) and fluctuates rapidly. Zero magnetization.  Start well below critical temperature. Start all spins with equal value (all up or all down). A few small clusters of opposite spin form.  Start well below critical temperature. Start with random spins. Large clusters of each spin form; eventually simulation makes up its mind. Equally likely to have large clusters in up or down spin.  Start close to critical temperature. Large clusters form, but fluctuate very slowly.

Exact solution for Ising model known for 1D and 2D; NP-hard for 3d and nonplanar graphs. Models phase changes in binary alloys and spin glasses. Also models neural networks, flocking birds, and beating heart cells. Over 10,000+ papers published using the Ising model.

Java Simulation --- Calculate PI using Random Numbers import java.util.*; public class CalculatePI2 { public static boolean isInside (double xPos, double yPos) { boolean result; double distance = Math.sqrt((xPos * xPos) + (yPos * yPos)); if (distance < 1) result = false; return(distance < 1); } public static double computePI (int numThrows) { Random randomGen = new Random (System.currentTimeMillis()); int hits = 0; double PI = 0; for (int i = 0; i <= numThrows; i++) { double xPos = (randomGen.nextDouble()) * 2 - 1.0; double yPos = (randomGen.nextDouble()) * 2 - 1.0; if (isInside(xPos, yPos)) { hits++; } PI = (4 * (hits/numThrows)); } return PI; } public static void main (String[] args) { Scanner reader = new Scanner (System.in); System.out.println("This program approximates PI using the Monte Carlo method."); System.out.println("It simulates throwing darts at a dartboard."); System.out.print("Please enter number of throws: "); int numThrows = reader.nextInt(); double PI = computePI(numThrows); double Difference = PI - Math.PI;

System.out.println ("Number of throws = " + numThrows + ", Computed PI = " + PI + ", Difference = " + Difference ); } }

Introduction A class of computational algorithms that rely on repeated random sampling to compute their results are called the Monte Carlo methods. Monte Carlo methods are often used when simulating physical and mathematical systems. Because of their reliance on repeated computation and random or pseudo-random numbers, Monte Carlo methods are most suited to calculation by a computer. Monte Carlo methods tend to be used when it is infeasible or impossible to compute an exact result with a deterministic algorithm. Monte Carlo algorithms work based on the Law of Large Numbers. It says that if you generate a large number of samples, eventually, you will get the approximate desired distribution. Monte Carlo methods have three characteristics: 1. 2. 3.

Random sample generation The input distribution is known Numerical experiments The direct output of the Monte Carlo simulation method is the generation of random sampling. Other performance or statistical outputs are indirect methods which depend on the applications. There are many different numerical experiments that can be done, probability distribution is one of them. Probability distribution identifies either the probability of each value of an unidentified random variable (when the variable is discrete), or the probability of the value falling within a particular interval (when the variable is continuous). That is equivalent to saying that for random variables X with the distribution in question, Pr[X = a] = 0 for all real numbers a. That is, the probability that X attains the value a is zero, for any number a. If the distribution of X is continuous, then X is called a continuous random variable. Normal distribution, continuous uniform distribution, beta distribution, and Gamma distribution are well known absolutely continuous distributions. In the sample program, I will be using a Triangular distribution, which is is a continuous probability distribution with a lower limit a, mode c, and upper limit b.

Triangular distribution is therefore often used in business decision making, particularly in simulations. Generally, when not much is known about the distribution of an outcome (say, only its smallest and largest values), it is possible to use uniform distribution. But, if the most likely outcome is also known, then the outcome can be simulated by a Triangular distribution. Whilst Monte Carlo methods can be a valuable aid to decision making, they should not be used without an understanding of the process being modeled. Typically, in any activity, there are a few critical inputs, and it is these which should be the focus of the analysis, not trying to attribute complex distributions to minor variables. So, a Monte Carlo analysis should be preceded by a sensitivity analysis to determine what the important parameters are.

Background There are many ways of approaching a Monte Carlo analysis, but a good starting point is a flow chart of the processing being modeled. For the purposes of this simulation, I will be focus on a tennis tournament simulation based on game winnings. The flow chart of the process goes like this:

Using the code

Here is an example where we have two players, each with varying skills and talent. Assume they each have played 10 tournaments and we have gotten the minimum and maximum number of wins (assuming winning the tournament has the same number of wins as all the other tournaments). We can also compute the mean win for all the 10 tournaments they have played. Collapse | Copy Code

double[] double[] double[] double[] double[]

VarWins = new double[]{player1,player2}]; minWins = new double[]{11,12}; meanWins = new double[]{12.5,13}; maxWins = new double[]{17,15}; results = MonteCarlo.simulate(2, new double[] { 11, 12 }, new double[] { 12.5, 13 }, new double[] { 17, 15 });

The main simulation is done using the following triangular distribution: Collapse | Copy Code

public static double triangular(double Min,double Mode,double Max) { // Declarations double R=0.0; // Initialise Random r = new Random(); R = r.NextDouble(); //between 0.0 and 1.0 gaussian // Triangular if ( R == (( Mode - Min) / ( Max - Min))) { return Mode; } else if ( R < (( Mode - Min) / ( Max - Min))) { return Min + Math.Sqrt( R * ( Max - Min) * ( Mode - Min)); else { return Max - Math.Sqrt((1 - R) * ( Max - Min) * ( Max - Mode)); } }

As seen in the output, player 1 is the most likely to win. However, player 2 is not far behind. So, with a little more practice and training, player 2 manages to reduce the mean number of wins on the next 10 tournaments. Simulating this again will show that he has greatly improved his performance and is most likely to win over player 1. Output: Collapse | Copy Code

Simulation Results Wins5251 Wins4749 Player1 win probability 52.51 - Player2 win probability47.49

Output after training:

Collapse | Copy Code

Simulation Results Wins4119 Wins5881 Player1 win probability 41.19 - Player2 win probability58.81:

There are some instances where it is the other way around. As I have mentioned ealier, by the Law of Large Numbers, if you increase the number of simulated tournaments, you will notice that either player will always be the winner because the expected result will always come closer because of the large number of simulations. More examples on the Monte Carlo algorithm can be found here.

Application Monte Carlo methods are often used in Physics to determine neutron trajectories or to simulate atomic clusters. In the 1950s, Monte Carlo methods were used to study nuclear cascades. When a particle collides with a nucleus, a nuclear cascade is produced. The particle’s path after the collision determines whether or not a new particle, a pion, is created. Monte Carlo methods keep track of the colliding particle’s path after the collision and determine the probability of the production of pions. Monte Carlo algorithms are useful here because it is difficult and time-consuming to simulate the collisions. By running the algorithm enough times, the results are usually statistically significant, and we are able to successfully calculate the probability of the production of nuclear cascades. Chemistry is another field where Monte Carlo methods come in handy. Monte Carlo methods are often used to evaluate integrals. While integrals can be evaluated by many other mathematical methods, it is, sometimes too complex for other methods to evaluate an integral as its dimensionality grows. Not being able to evaluate an integral with growing dimensions is known as the ‘curse of dimensionality’. However, Monte Carlo methods outdo this curse because they work increasingly better with more complex integrals that contain high dimensions. Evaluation of multidimensional integrals is useful in calculating the volume of molecules. Since molecules do not have defined boundaries, volume calculations involve multidimensional integration. It is difficult to calculate the volume without the use of Monte Carlo methods because of the complexity involved and also because definite boundaries of the molecules are undefined. Volumes of molecules are often defined by a range of electron densities, so Monte Carlo methods use random number generation and a probability distribution function to evaluate the volumes of molecules. One of the most popular applications of the Monte Carlo algorithms is in the field of finance. Monte Carlo methods aid the analysis of financial instruments, portfolios, and assets. When a financial instrument or asset is being analyzed to label its value, many complex equations, the values of which may be uncertain,

are used to reach a final answer. Since Monte Carlo methods work well with highly complex equations, their use becomes vital in the calculation of uncertain values, which then in turn help analyze the final value of the instrument or asset in question. A specific ‘Monte Carlo Option Model’ is used to evaluate future prices of options. Many uncertain values affect the final value of these financial options; Monte Carlo methods use random number generation to lay the various price paths and then calculate a final option value.

/* Code is in Monte Carlo.zip */