Advanced Monte Carlo Methods:Direct SimulationProf. Dr. Michael MascagniE-mail: [email protected], http://www.cs.fsu.edu/~mascagniSeminar für Angewandte MathematikEidgenössische Technische Hochschule, CH-8044 Zürich SwitzerlandandDepartment of Computer Science and School of Computational ScienceFlorida State University, Tallahassee, Florida, 32306-4120 USADirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 2 of 16 Direct Simulation Probabilistic Problem or Model Directly simulate the physical random processes of the original problem Replace complicated “blocks” with randomized outcomes (neutronics) Direct the simulation’s random outcomes with random numbers Examples Controlling floodwater and construction of new dams on the Niles The quantity of water in the river varies randomly from season to season Use records of weather, rainfall, and water levels extending over many years Examine what may happen to the water if certain dams are built and certain possible policies of water control are exercised Evaluate artificial lake impact on people, agriculture, transportation, antiquitiesDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 3 of 16 Direct Simulation (Cont.) Examples (Cont.) Telephone networks Computer networks Growth of an insect population on the basis of certain assumed vital statistics of survival and reproduction Service rates (queuing theory) Operations research problems Actuarial simulations Insurance risk from tabulated date (mortality) Risk assessment of investment portfolios Computer games Roadway design simulation War gamingDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 4 of 16 Direct Simulation of the Lifetime of Comets A Long-period Comet Described as a sequence of elliptic orbits Sun at one focus of the orbit ellipse Energy of a comet is inversely proportional to the length of the semi-major axis of the ellipseDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 5 of 16 Direct Simulation of the Lifetime of Comets (Cont.) Behavior of the Comet Most of the time Moves at a great distance from the sun A relatively short time (instantaneous) Passes through the immediate vicinity of the sun and the planets At this instant, the gravitational field of the planet perturbs the cometary energy by a random component Successive energy perturbations (in suitable units of energy) may be taken as independent normally distributed random variables η1, η2, η3,…, ηn η1, η2, η3,…, ηncan be normalized to, for example, a standard normal distribution, N(0,1) Computing the perturbations exactly is dauntingDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 6 of 16 Direct Simulation of the Lifetime of Comets (Cont.) Comets under perturbation A comet, starting with an energy, G = –z0, has subsequent energies -z0, -z1=-z0+η1, -z2=-z1+η2, … The process continues until the first occasion on which z changes sign (negative = bound state; positive = free state) Once z changes sign, the comet departs on a hyperbolic orbit and is lost from the solar system Kepler’s Third Law the time taken to describe an orbit with energy –z is z-3/2 the total lifetime of the comet is zTis the first negative quantity in the sequence of z0, z1, …Direct SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 7 of 16 Direct Simulation of the Lifetime of Comets (Cont.) Problem is to determine distribution of G given z0 Very difficult problem in theoretical mechanics Easy to simulate using probabilistic model Seek CDF:Direct SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 8 of 16 Direct Simulation of the Lifetime of Comets (Cont.) Monte Carlo trick: use analytic/deterministic information when it is available Probability of escape in one orbit (available from N(0,1) tables: Know also:Direct SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 9 of 16 Direct Simulation of the Lifetime of Comets (Cont.) Now we need to estimate:N*samples (subsample) out of N have T > 1 Let 1-p*(g) be the proportion of values in the subsample with G>gDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 10 of 16 Direct Simulation of the Lifetime of Comets (Cont.) Using conditional probability we have:Direct SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 11 of 16 A Robotics Problem Take symmetric, concentric objects, O1and O2 in random orientations: what is the distribution of the smallest angle needed to rotate one into coincidence with the other? There is an analytic solution, but for simulation need to generate random orientations via random 3x3 orthogonal matricesDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 12 of 16 A Robotics Problem (Cont.) Let Q0= [x; y; z] where the vectors x, y, z are uniformly distributed on S2x = (x1, x2, x3)| x |-1 is uniform on S2 when x1, x2, x3 are i.i.d. N(0,1) random variables Similarly define y*=(y1, y2, y3)| y*|-1 Take y = (y*–Px)/(1 – P2)1/2 where P = x • y* (Gram-Schmidt procedure)Direct SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 13 of 16 A Robotics Problem (Cont.) Finally, let z = x x y which is orthogonal to the othersDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 14 of 16 Dimensional Analysis Consider the “Traveling Salesman” problem: given n towns to visit, no order, minimize the Hamiltonian path through the Euclidean weighted complete graph Let l be the length of the shortest path A total area of region containing cities n/A is the density of citiesDirect SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 15 of 16 Dimensional Analysis (Cont.) Assume l = Aa(n/A)b= nbA(a-b)Units l –length n – dimensionless (number) A –(length)2 Implies a-b = ½Direct SimulationProf. Dr. Michael Mascagni: Advanced Monte Carlo MethodsSlide 16 of 16 Dimensional Analysis (Cont.) Multiply the area by f while keeping the density constant, l is also multiplied by f fA replaces A fn replaces n fl replaces l l = nbA1/2implies fl = fbnbf1/2A1/2so b = ½ l = k