The Monte Carlo method

1. The Monte Carlo method

2. What is Monte Carlo? • A method to evaluate integrals • Popular in all areas of science, economics, sociology, etc. • Involves: • Random sampling of domain • Evaluation of function • Acceptance/rejection of points

3. Why Monte Carlo? • The domain of integration is very complex: • Find the integral of a complicated solvent-solute energy function outside the boundary of a protein. • Find the volume of an object from a tomographic data set. • The integrals are very high-dimensional: • A partition function or free energy for a protein • Analyzing metabolic networks • Predicting trends in the stock market

4. Basic Monte Carlo • Numerical integral: Divide the region [0,1] evenly into M slices • Equivalent to sampling a [0,1] uniform random number • Error typically decays as

5. Error analysis • Standard deviation: Error from the Monte Carlo Error from trapezoid rule High dimensional problems: statistical physics, many body problem For d>10, Monte Carlo method is the Unique workable quadrature.

6. Basic Monte Carlo Algorithm • Suppose we want to approximate in a high-dimensional space • For i = 1 to n • Pick a point xiat random • Accept or reject the point based on criterion • If accepted, then add f(xi) to total sum • Error estimates are “free” by calculating sums of squares • Error typically decays as

7. Example: the area of a circle • Sample points randomly from square surrounding circle of radius 5 • 10,000 sample points • Acceptance criterion: inside circle • Actual area: 78.54 • Calculated area: 78.66

8. Example: more complicated shapes

9. Example: multiple dimensions • What is the average of a variable for a N-dimensional probability distribution? • Two approaches: • Quadrature • Discretize each dimension into a set of n points • Possibly use adaptivity to guide discretization • For a reasonably smooth function, error decreases as n-N/2 • Monte Carlo • Sample m points from the space • Possibly weight sampling based on reference function • Error decreases as m-1/2

10. Problems: sampling tails of distributions • We want to • Integrate a sharply-peaked function • Use Monte Carlo with uniformly-distributed random numbers • What happens? • Very few points contribute to the integral (~9%) • Poor computational efficiency/convergence • Can we ignore the tails? NO! • Solution: use a different distribution

11. Improved sampling: change of variables • One way to improve sampling is to change variables: • New distribution is flatter • Uniform variates more useful • Advantages: • Simplicity • Very useful for generating distributions of non-uniform variates (coming up) • Disadvantages • Most useful for invertible functions

12. Change of variables: method • Given an integral • Transform variables • Choose variables to give (nearly) flat distribution • Integrate

13. Change of variables application: exponential • Given an integral • Transform variables • Choose variables to give (nearly) flat distribution • Integrate

14. Change of variables application: exponential • Before transformation: • 0.5% of points in domain contribute to integral • Slow convergence

15. Change of variables example: exponential • Before transformation: • 0.5% of points in domain contribute to integral • Slow convergence • After transformation: • All of the points in domain contribute • Rapid (exact) convergence • In practice: • Speed-up best when inversion is nearly exact

16. Importance sampling • Functions aren’t always invertible • Is there another way to improve sampling of “important” regions of the functions? • Find flat distributions • Bias sampling • Find a function that is almost proportional to the one of interest: • Rewrite your integral as a “weighted” integral:

17. Importance sampling example: a lumpy Gaussian • Our original integrand is • This is close to • Therefore: • Sample random numbers from (we’ll talk about how to do this later): • Evaluate the following integrand over the random number distribution:

18. Importance sampling example: a lumpy Gaussian • Convergence is pretty good (actual value 1.41377…)

19. Evolution of Monte Carlo methods so far… • Uniform points and original integrand… • …but this had very poor efficiency…. • Uniform points and transformed integrand… • …but this only worked for certain integrands…. • Non-uniform points and scaled integrand… • …but this is very cumbersome for complicated integrands… • Now, we try Markov chain approaches…

20. Markov chains • Properties • A sequence of randomly-chosen states • The probability of transitions between states is independent of history • The entire chain represents a stationary probability distribution • Examples • Random number generators • Brownian motion • Hidden Markov Models • Financial market

21. Markov chain Monte Carlo • Assembling the entire distribution for MC is usually hard: • Complicated energy landscapes • High-dimensional systems • Extraordinarily difficult normalization • Solution: • Build up distribution from Markov chain • Choose local transition probabilities which generate distribution of interest (i.e., ensure detailed balance) • Each random variable is chosen based on the previous variable in the chain • “Walk” along the Markov chain until convergence reached • Result: Normalization not required, calculations are local

26. Detailed balance: example • Markov chain with a Gaussian stationary probability distribution • Detailed balance satisfied

29. A B Application: stochastic transitions • This is a very simplistic version of kinetic Monte Carlo… • Each Monte Carlo step corresponds to a time interval • The probability of moving between states in that time interval is related to the rate constants • Simulate to give mean first passage times, transient populations, etc.

30. Metropolis Monte Carlo • Start with the detailed balance condition • Derive an “acceptance ratio” condition • Choose a particular acceptance ratio

31. Application: canonical ensemble • Our un-normalized probabilities look like Boltzmann factors: • Our acceptance ratio is therefore:

32. Algorithm: NVT Metropolis Monte Carlo Metropolis MC in a harmonic potential

33. Choose step size • [0,1] random number • Accepting rate ~50%

34. Advantages of Metropolis MC simulations • Does not require forces • Rapidly-changing energy functions • No differentiation required • Amenable to complex move sets: • Torsions • Rotamers • Tautomers • Etc.

35. Monte Carlo moves • Trial moves • Rigid body translation • Rigid body rotation • Internal conformational changes (soft vs. stiff modes) • Titration/electronic states • … • Questions: • How “big” a move should we take? • Move one particle or many?

36. Monte Carlo moves • How “big” a move should we take? • Smaller moves: better acceptance rate, slower sampling • Bigger moves: faster sampling, poorer acceptance rate • Amortize mean squared displacement with respect to CPU time • “Rules of thumb”: percent acceptance rate • Move one particle or many? • Possible to achieve more efficient sampling with correct multi-particle moves • One-particle moves must choose particles at random

37. Summary • A method for sampling • Integrals • Configuration space • Error space • Any distribution that can be represented by transitions can be sampled • Sampling can be accelerated using various tricks