Types of Simulations

• Progression of time
  • Time is constant – static simulation
  • Time varies – dynamic simulation
    • Time increases continuously – continuous simulation
    • Time increases in discrete time steps – discrete event simulation

Types of Simulations

• Simulation outcomes
  • No randomness, outcome is always identical – deterministic simulation
  • Randomness/probability may affect the outcome – non-deterministic simulation
    • What If simulations - randomized initial state, predictable updates
    • Stochastic simulations - randomized updates

Types of Simulations

• Stochastic simulation involves using random values, combined with probabilities, iteratively
  • Necessarily, such a simulation will yield different results each time it is executed
  • However, certain trends should repeat themselves across multiple simulations

Why Stochastic Simulation?

• Say we want to model the relationship between predator and prey
  • Presumably, we know the death rates and birth rates for each, as well as the rate of consumption of prey by predators
  • We could perform a deterministic simulation based on these probabilities over a population

pred = pred  - (pred * deathRate)
                            

• An example that works this way is the Lotka-Volterra model

Why Stochastic Simulation?

• A rate such as deathRate is not constant
  • If the value was determined empirically, it is likely an average over values that change over time
  • If the change in value can be represented by a differential equation, then we can include it as a dependent variable in our simulation
  • What if the phenomenon is not well understood?
• A stochastic simulation, in contrast, can simulate fluctuations over time

Random Numbers

• Computers are not very good at creating random numbers
  • Nearly all functionality in a computer is deterministic, even pseudorandom number generators (PRNGs)
• To get truly random values, you often require an external device to obtain entropy (randomness) from the environment
  • e.g. A device that measures quantum oscillations to generate random data
  • e.g. A device that measures fluctuations in biometric data
• NumPy has a pseudorandom number generator that is random enough for our purposes

Random Numbers: Uniform Distribution

• Often, we want our random numbers to come from a uniform probability distribution
  • i.e. any number is equally likely to be chosen
• A uniform probability distribution looks like this:

Random Numbers: Gaussian Distribution

• Sometimes, we want our random numbers to be chosen from another, non-uniform, distribution
  • Such distributions are characterized by a probability which varies depending upon the value
  • e.g. grades are normally distributed along a Gaussian (or normal) distribution
• A Gaussian probability distribution looks like this:

Random Numbers: Triangular Distribution

• A triangular distribution is used when little is known about the distribution of values
  • We don't have a set of actual values (e.g. prediction)
  • It uses the minimum, maximum, and most likely (often a mode) value
  • The shape doesn't describe the actual distribution exactly, but it is a best estimate, based on limited data
• A triangular probability distribution looks like this:

Random Numbers: Triangular Distribution

• There are plenty of other distributions:

Monte Carlo Simulations

• For example, imagine there is a 30% chance that a predator will eat a prey in its vicinity
  • One could generate a random number between 0.0 and 1.0
  • If the random number is between 0.0 and 0.3, then the predator eats the prey
  • If the random number is between 0.3 an 1.0, then the prey escapes
• A stochastic simulation that follows this method repeatedly is called a Monte Carlo simulation

What If Simulations

• In contrast to Monte Carlo simulations is a what if simulation
  • Randomize the initial conditions
  • Use deterministic methods to update state variables
• What if simulations have limited use compared to Monte Carlo technique

Monte Carlo Simulations

1. Define a distribution of possible inputs for each input random variable
2. Generate inputs randomly from these distributions
3. Perform a deterministic computation using that set of inputs
4. Aggregate the results of the individual computations into the final result

Monte Carlo Simulations

1. Define a distribution of possible inputs for each input random variable
   • Use a known distribution
   • Estimate the distribution, based on empirical evidence
2. Generate inputs randomly from these distributions
3. Perform a deterministic computation using that set of inputs
4. Aggregate the results of the individual computations into the final result

Monte Carlo Simulations

1. Define a distribution of possible inputs for each input random variable
2. Generate inputs randomly from these distributions
   • For each dependent variable, generate a vector of inputs by repeatedly sampling the variable's distribution
3. Perform a deterministic computation using that set of inputs
4. Aggregate the results of the individual computations into the final result

Monte Carlo Simulations

1. Define a distribution of possible inputs for each input random variable
2. Generate inputs randomly from these distributions
3. Perform a deterministic computation using that set of inputs
   • Iteratively, choose one of the randomized input values from the vector, and applying equations to update state variables based on these input values
   • Continue until all randomized input values from the vector have been exhausted
4. Aggregate the results of the individual computations into the final result

Monte Carlo Simulations

1. Define a distribution of possible inputs for each input random variable
2. Generate inputs randomly from these distributions
3. Perform a deterministic computation using that set of inputs
4. Aggregate the results of the individual computations into the final result
   • Collect all values for the state variables
   • Save, analyze, and/or visualize the resulting data

Example: Rolling Doubles

• Let's try to figure out how likely it is to roll doubles using dice
  • This problem is relatively easy to solve using probability
    • 6 possible ways to achieve doubles, out of 62=36 possible roll combinations = 1/6 = 0.166666…
  • …but let's use stochastic simulation
• Let's roll a large number of pairs of dice, and see how many of them were doubles

Example: Rolling Doubles

• Here is an implementation:

def probabilityOfDoubles(iterations):
   numRolls = 0
   numDoubles = 0
   for i in range(iterations):
      die1 = int(np.random.rand() * 6 + 1)
      die2 = int(np.random.rand() * 6 + 1)
      numRolls = numRolls + 1

      if die1 == die2:
         numDoubles = numDoubles + 1

   return float(numDoubles) / float(numRolls)
                        

Aside: Handling Probabilities

• If there is an event which occurs with probability p:
  • Generate a random number r between 0.0 and 1.0
    • The random number should generally be taken from a uniform distribution
    • If r ≤ p, then the event occurred
    • Otherwise, the event did not occur

Example: Atomic Decay

• Some unstable elements decay, probabilistically, at a predictable rate
  • Let the probability of an atom decaying in some time period, Δt, be P
  • After one time period, Δt, the probability of the atom intact is (1-P)
  • After two time periods, 2*Δt, the probability of the atom intact is (1-P)²
  • …
  • After k time periods, k*Δt, the probability of the atom intact is (1-P)^k

Example: Atomic Decay

• Basic strategy:

Given numAtoms, decayProb, and timeSteps:
For each of timeSteps:
    For each of numAtoms:
        Determine if this atom will decay (decayProb)
        If yes, then decrease numAtoms
                        

Example: Atomic Decay (Stochastic)

• To simulate each atom:

def simulateHalfLife(numAtoms, decayProb, timeSteps):
   initialAtoms = numAtoms
   for t in xrange(timeSteps):
      for m in xrange(numAtoms):
         if np.random.rand() < decayProb:
            numAtoms -= 1
   return initialAtoms, numAtoms
                        

Example: Atomic Decay (Deterministic)

• For comparison, here is a deterministic function, which plots the number of atoms over time:

def plotDecay(numAtoms, decayProb, timeSteps):
   numRemaining = [numAtoms]
   for t in range(timeSteps):
      count = numAtoms*((1-decayProb)**t)
      numRemaining.append(count)
   plt.plot(numRemaining)
   plt.xlabel("Time Steps")
   plt.ylabel("Atoms")
   plt.grid(True)
   plt.show()

                        

Example: Estimating π

• One clever way to estimate π proposed by Georges-Louis LeClerc:
  • Imagine a (radius 1) circle inside a (width 2) square

\begin{equation}\label{01} \begin{split} \frac{area_{circle}}{area_{square}} &= \frac{π\times r^2}{w^2} = \frac{π\times1^2}{2^2} = \frac{π}{4} \end{split} \end{equation}

Example: Estimating π

• Suppose we start dropping objects randomly into the square
  • What is the probability of dropping the object into each shape?
  • The probabilities are proportional to the area of their shapes:

\begin{equation}\label{02} \begin{split} prob_{square} &= 1.0\\ prob_{circle} &= \frac{π}{4}\\ π &\simeq 4\times\frac{prob_{circle}}{prob_{square}}\\ \end{split} \end{equation}

Example: Estimating π

• Here is an example run:

\begin{equation}\label{03} \begin{split} π_{estimated} &= 4\times\frac{156}{200} = 3.12\\ \end{split} \end{equation}

Example: Estimating π

• Here is the implementation:

import math

def estimatePi(iterations):
	numInsideCircle = 0
	for i in range(iterations):
		if dropObject():
			numInsideCircle += 1
	return 4.0 * float(numInsideCircle) / float(iterations)

def dropObject():
	x = np.random.random()
	y = np.random.random()
	if math.sqrt(x**2 + y**2) < 1.0:
		return True
	else:
	return False
                        

Example: Estimating Earnings

• Earnings due to sales involves multiple variables, with uncertainty for each
• What if we don't know the distribution for each of these variables?
  • We can use a triangular distribution

Example: Estimating Earnings

• We could use a triangular distribution to estimate our earnings directly
  • We'd need to guess at minimum, maximum, and most likely earnings
  • This merely gives us a single figure, with limited variation
  • Plus, earnings is also more tricky to estimate, due to its dependence on the other variables
• A Monte Carlo method based on triangular distribution for the input variables is a better option

Example: Estimating Earnings

• Generate triangular distributions for:
  • unitPrice
  • unitSales
  • fixedCosts
  • variableCosts
• Calculate the earnings by taking random numbers from these distributions and plugging them into the formula:
  \begin{equation}\label{04} \begin{split} earnings_{estimated} = unitCost\times unitSales - varCosts - fixedCosts\\ \end{split} \end{equation}
• Plot the distribution of earnings, to give us a good idea of what to expect

Random Walks

• A random walk is named after the analogy typically used to describe it:
  • Imagine a drunken bar patron late at night
  • At each time step, they take a step in a random direction
  • Where do they end up?

Random Walks

• Random walk simulations can be used to describe many phenomena:
  • Brownian motion of particles suspended in fluids
  • Motion of bacteria
  • Effect of random mutations on a genome
  • Proportions of nucleotide bases required during DNA replication
  • Motion of entities within an interactive simulation/game
  • Motion of users within a mobile network
  • Motion of animals in a eco-simulation

Random Walks

• Named after 19^th century botanist Robert Brown
  • Trying to explain the apparently random motion of particulate matter suspended in a fluid
    • e.g. pollen suspended in water
  • Was the first to recognize that the random motion was not caused by living organisms
    • This was the prevailing theory at the time
    • Though, he could not explain the cause of the phenomena

Random Walks

• General case – a cellular automaton
  • An n-dimensional space, divided into discrete cells like a grid
    • The cell within that space represents the current state
    • For the simulation to be visually useful, we need to restrict n to [1, 2, 3]
  • A set of transition rules
    • Describe how to change state, represented as a position change in the n-dimensional space
    • In a random walk, the transition rules are selected randomly from the set of currently applicable choices