Stochastic simulation

A probability distribution is used to describe the potential outcome of a random variable.

Limits the outcomes where the variable can only take on discrete values.[4]

A random variable X is Bernoulli-distributed with parameter p if it has two possible outcomes usually encoded 1 (success or default) or 0 (failure or survival)[5] where the probabilities of success and failure are ${\displaystyle P(X=1)=p}$ and ${\displaystyle P(X=0)=1-p}$ where ${\displaystyle 0\leq p\leq 1}$.

To produce a random variable X with a Bernoulli distribution from a U(0,1) uniform distribution made by a random number generator, we define

${\displaystyle X={\begin{cases}1,&{\text{if }}0\leq U<p\\0,&{\text{if }}1\geq U\geq p\end{cases}}}$

such that the probability for${\displaystyle P(X=1)=P(0\leq U<p)=p}$and ${\displaystyle P(X=0)=P(1\geq U\geq p)=1-p}$.[2]

Define

X = 1 if head comes up and
X = 0 if tail comes up

For a fair coin, both realizations are equally likely. We can generate realizations of this random variable X from a ${\displaystyle U(1,0)}$ uniform distribution provided by a random number generator (RNG) by having ${\displaystyle X=1}$ if the RNG outputs a value between 0 and 0.5 and${\displaystyle X=0}$ if the RNG outputs a value between 0.5 and 1.

P (X = 1) = P(0 ≤ U < 1/2) = 1/2

P (X = 0) = P(1 ≥ U ≥ 1/2) = 1/2

Of course, the two outcomes may not be equally likely (e.g. success of medical treatment).[6]

A binomial distributed random variable Y with parameters n and p is obtained as the sum of n independent and identically Bernoulli-distributed random variables X₁, X₂, ..., X_n[4]

Example: A coin is tossed three times. Find the probability of getting exactly two heads. This problem can be solved by looking at the sample space. There are three ways to get two heads.

HHH, HHT, HTH, THH, TTH, THT, HTT, TTT

The answer is 3/8 (= 0.375).[7]

A poisson process is a process where events occur randomly in an interval of time or space.[2][8] The probability distribution for poisson processes with constant rate λ per time interval is given by the following equation.[4]

${\displaystyle P(k{\text{ events in interval}})={\frac {\lambda ^{k}e^{-\lambda }}{k!}}}$

Defining ${\displaystyle N(t)}$ as the number of events that occur in the time interval ${\displaystyle t}$

${\displaystyle P(N(t)=k)={\frac {(t\lambda )^{k}}{k!}}e^{-t\lambda }}$

It can be shown that inter-arrival times for events is exponentially distributed with a cumulative distribution function(CDF) of ${\displaystyle F(t)=1-e^{{-t}{\lambda }}}$. The inverse of the exponential CDF is given by

${\displaystyle t=-{\frac {1}{\lambda }}ln(u)}$

where ${\displaystyle u}$ is an ${\displaystyle U(0,1)}$ uniformly distributed random variable.[2]

Simulating a Poisson process with a constant rate ${\displaystyle \lambda }$ for the number of events ${\displaystyle N}$ that occur in interval ${\displaystyle [t_{start},t_{end}]}$ can be carried out with the following algorithm.[9]

1. Begin with ${\displaystyle N=0}$ and ${\displaystyle t=t_{start}}$
2. Generate random variable ${\displaystyle u}$ from ${\displaystyle U(0,1)}$ uniform distribution
3. Update the time with ${\displaystyle t=t+(-1/\lambda )ln(u)}$
4. If ${\displaystyle t>t_{end}}$, then stop. Else continue to step 5.
5. ${\displaystyle N=N+1}$
6. Continue to step 2

Published by Dan Gillespie in 1977, and is a linear search on the cumulative array. See Gillespie algorithm.

Gillespie’s Stochastic Simulation Algorithm (SSA) is essentially an exact procedure for numerically simulating the time evolution of a well-stirred chemically reacting system by taking proper account of the randomness inherent in such a system.[10]

It is rigorously based on the same microphysical premise that underlies the chemical master equation and gives a more realistic representation of a system’s evolution than the deterministic reaction rate equation (RRE) represented mathematically by ODEs.[10]

As with the chemical master equation, the SSA converges, in the limit of large numbers of reactants, to the same solution as the law of mass action.

Published 2000 by Gibson and Bruck.[11] This is an improvement over the first reaction method where the unused reaction times are reused. To make the sampling of reactions more efficient, an indexed priority queue is used to store the reaction times. On the other hand, to make the recomputation of propensities more efficient, a dependency graph is used. This dependency graph tells which reaction propensities to update after a particular reaction has fired.

Published 2004[12] and 2005. These methods sort the cumulative array to reduce the average search depth of the algorithm. The former runs a presimulation to estimate the firing frequency of reactions, whereas the latter sorts the cumulative array on-the-fly.

Published in 2006. This is a binary search on the cumulative array, thus reducing the worst-case time complexity of reaction sampling to O (log M).

Published in 2009, 2010, and 2011 (Ramaswamy 2009, 2010, 2011). Use factored-out, partial reaction propensities to reduce the computational cost to scale with the number of species in the network, rather than the (larger) number of reactions. Four variants exist:

• PDM, the partial-propensity direct method. Has a computational cost that scales linearly with the number of different species in the reaction network, independent of the coupling class of the network (Ramaswamy 2009).
• SPDM, the sorting partial-propensity direct method. Uses dynamic bubble sort to reduce the pre-factor of the computational cost in multi-scale reaction networks where the reaction rates span several orders of magnitude (Ramaswamy 2009).
• PSSA-CR, the partial-propensity SSA with composition-rejection sampling. Reduces the computational cost to constant time (i.e., independent of network size) for weakly coupled networks (Ramaswamy 2010) using composition-rejection sampling (Slepoy 2008).
• dPDM, the delay partial-propensity direct method. Extends PDM to reaction networks that incur time delays (Ramaswamy 2011) by providing a partial-propensity variant of the delay-SSA method (Bratsun 2005, Cai 2007).

The use of partial-propensity methods is limited to elementary chemical reactions, i.e., reactions with at most two different reactants. Every non-elementary chemical reaction can be equivalently decomposed into a set of elementary ones, at the expense of a linear (in the order of the reaction) increase in network size.

A general drawback of stochastic simulations is that for big systems, too many events happen which cannot all be taken into account in a simulation. The following methods can dramatically improve simulation speed by some approximations.

Since the SSA method keeps track of each transition, it would be impractical to implement for certain applications due to high time complexity. Gillespie proposed an approximation procedure, the tau-leaping method which decreases computational time with minimal loss of accuracy.[13] Instead of taking incremental steps in time, keeping track of X(t) at each time step as in the SSA method, the tau-leaping method leaps from one subinterval to the next, approximating how many transitions take place during a given subinterval. It is assumed that the value of the leap, τ, is small enough that there is no significant change in the value of the transition rates along the subinterval [t, t + τ]. This condition is known as the leap condition. The tau-leaping method thus has the advantage of simulating many transitions in one leap while not losing significant accuracy, resulting in a speed up in computational time.[14]

This method approximates reversible processes (which includes random walk/diffusion processes) by taking only net rates of the opposing events of a reversible process into account. The main advantage of this method is that it can be implemented with a simple if-statement replacing the previous transition rates of the model with new, effective rates. The model with the replaced transition rates can thus be solved, for instance, with the conventional SSA.[15]

The random variable X is said to be normally distributed with parameters μ and σ, abbreviated by X ∈ N (μ, σ²), if the density of the random variable is given by the formula [4] ${\displaystyle f_{X}(x)={\frac {1}{\sqrt {2\pi \sigma ^{2}}}}e^{-{\frac {(x-\mu )^{2}}{2\sigma ^{2}}}}.}$ x ∈ R.[4]

Many things actually are normally distributed, or very close to it. For example, height and intelligence are approximately normally distributed; measurement errors also often have a normal distribution.[19]

Exponential distribution describes the time between events in a Poisson process, i.e. a process in which events occur continuously and independently at a constant average rate.

The exponential distribution is popular, for example, in queuing theory when we want to model the time we have to wait until a certain event takes place. Examples include the time until the next client enters the store, the time until a certain company defaults or the time until some machine has a defect.[4]

Student's t-distribution are used in finance as probabilistic models of assets returns. The density function of the t-distribution is given by the following equation:[4] ${\displaystyle f(t)={\frac {\Gamma ({\frac {\nu +1}{2}})}{{\sqrt {\nu \pi }}\,\Gamma ({\frac {\nu }{2}})}}\left(1+{\frac {t^{2}}{\nu }}\right)^{-{\frac {\nu +1}{2}}},\!}$

where ${\displaystyle \nu }$ is the number of degrees of freedom and ${\displaystyle \Gamma }$ is the gamma function.

For large values of n, the t-distribution doesn't significantly differ from a standard normal distribution. Usually, for values n > 30, the t-distribution is considered as equal to the standard normal distribution.

• Generalized extreme value distribution

There are wide possibilities for use of Monte Carlo Method:[1]

• Statistic experiment using generation of random variables (e.g. dice)
• sampling method
• Mathematics (e.g. numerical integration, multiple integrals)
• Reliability Engineering
• Project Management (SixSigma)
• Experimental particle physics
• Simulations
• Risk Measurement/Risk Management (e.g. Portfolio value estimation)
• Economics (e.g. finding the best fitting demand curve)
• Process Simulation
• Operations Research