Montecarlo Simulations¶
evolve the state of a system given a probability of transition from each state
The master equation example: I have a system where I want to keep track of the number of molecules in an environment.
the number of molecules can change due to two different processes:
- a particle from the outside can come in with a probability that is constant (or can be approximated as such)
- a particle can escape the system, with a probability that is constant for each particle (so the total probability is proportional to the number of particles)
import numpy as np
import pylab as plt
import random as rn
def removal(state):
return state*0.1
def increase(state):
return 1
transitions = [removal, increase]
transitions_names = ['removal', 'increase']
state = 5
rates = [f(state) for f in transitions]
rates
the Stochastic Simulation Algorithm¶
the total rate $\lambda_{tot}$ is the sum of the rates
the chance of each one happening is the proportion compared to the total rate
the time for which the system keeps being in the same state is distributed with an exponential time $\tau = \lambda_{tot}^{-1}$
total_rate = sum(rates)
total_rate
time = np.random.exponential(1/total_rate)
time
choosing the transition¶
we can use the same sampling method as before!
normalized_rates = np.cumsum(rates)
normalized_rates /= normalized_rates[-1]
normalized_rates
p = np.random.rand()
idx = np.searchsorted(normalized_rates, p)
idx
event = transitions_names[idx]
event
With python 3.6 onward, we can actually replace this with a simple call to random.choices, that handles weighted extractions
rn.choices(transitions_names, weights=rates)
event = rn.choices(transitions_names, weights=rates)[0]
print(event)
updating the status¶
given the chosen event we can update the current status and then start all over
if event == 'increase':
state += 1
elif event == 'removal':
state -= 1
else:
raise ValueError("transition not recognized")
the complete cycle¶
the complete simulation requires to keep track of which transition happened, how much time the system has spent in each status and so on.
one other thing we might want to keep track of is the total time ellapsed in the simulation.
When should we stop our simulation?
ideally weshould stop after a certain amount of "simulated time" has passed
import typing
from enum import Enum
class Transition(Enum):
INCREASE = 'increase'
DECREASE = 'removal'
class Observation(typing.NamedTuple):
state: typing.Any
time_of_observation: float
time_of_residency: float
transition: Transition
Observation(5, 1.0, 0.33, Transition.DECREASE)
transitions = [removal, increase]
transitions_names = [Transition.DECREASE, Transition.INCREASE]
observed_states = []
state = 5
total_time = 0.0
time_limit = 5.0
while total_time < time_limit:
rates = [f(state) for f in transitions]
total_rate = sum(rates)
time = np.random.exponential(1/total_rate)
event = rn.choices(transitions_names, weights=rates)[0]
observation = Observation(state, total_time, time, event)
print(observation)
observed_states.append(observation)
total_time += time
if event == Transition.INCREASE:
state += 1
elif event == Transition.DECREASE:
state -= 1
else:
raise ValueError("transition not recognized")
let's make it a function so we have a more general way of managing it
def simulation(starting_state, time_limit):
observed_states = []
state = starting_state
total_time = 0.0
while total_time < time_limit:
rates = [f(state) for f in transitions]
total_rate = sum(rates)
time = np.random.exponential(1/total_rate)
event = rn.choices(transitions_names, weights=rates)[0]
observation = Observation(state, total_time, time, event)
observed_states.append(observation)
total_time += time
if event == Transition.INCREASE:
state += 1
elif event == Transition.DECREASE:
state -= 1
else:
raise ValueError("transition not recognized")
return observed_states
time_limit = 100
result = simulation(starting_state=5, time_limit=time_limit)
len(result)
for res in result[-5:]:
print(res)
One thing we can observe is that the last state overshoot the time limit a bit for any simulation we could do.
The solution to this is to trim the time the system stays in the state so that the final overall time is exactly equal to the given limit.
Being processes with constant probability in time this does not deform the probabilities.
for now we can ignore it, but will come back later to this, because it is also related to problem of absorbing states and fixed time events.
final probability distribution¶
to estimate the probability of observing our system in a given state, we have to sum up the time the system has spent in each state.
BEWARE: the probability depends on the time of residency, not on the number of times that state has appeared!
the expected distribution for this process is a poisson distribution with average 10
from collections import Counter
distribution = Counter()
for observation in result:
state = observation.state
residency_time = observation.time_of_residency
distribution[state] += residency_time/time_limit
print(distribution)
import scipy.stats as st
fig, ax = plt.subplots()
ax.bar(distribution.keys(), distribution.values())
values = np.arange(20)
pmf = st.poisson(10).pmf(values)
ax.bar(values, pmf, alpha=0.5)
let's make it run for longer...
def generate_distribution(observation_sequence):
distribution = Counter()
for observation in observation_sequence:
state = observation.state
residency_time = observation.time_of_residency
distribution[state] += residency_time
total_time_observed = sum(distribution.values())
for state in distribution:
distribution[state] /= total_time_observed
return distribution
time_limit = 1_000
result = simulation(starting_state=5, time_limit=time_limit)
distribution = generate_distribution(result)
fig, ax = plt.subplots()
ax.bar(distribution.keys(), distribution.values())
values = np.arange(20)
pmf = st.poisson(10).pmf(values)
ax.bar(values, pmf, alpha=0.5)
warn-in period¶
for monodimensional systems it's relatively easy to converge to the right distribution, but the starting point can have a substancial effect!
time_limit = 1_000
result = simulation(starting_state=50, time_limit=time_limit)
distribution = generate_distribution(result)
fig, ax = plt.subplots()
ax.bar(distribution.keys(), distribution.values())
values = np.arange(20)
pmf = st.poisson(10).pmf(values)
ax.bar(values, pmf, alpha=0.5)
usually for these kind of simulations the first period of time gets removed to avoid problem of convergence
multiple simulations¶
Estimating the convergence to the final distribution can be quite hard with a single simulation, so a good rule of thumb is to run at least 3 or 4 identical simulations to assess how much the convergence has been reached
time_limit = 1_000
result_1 = simulation(starting_state=50, time_limit=time_limit)
result_2 = simulation(starting_state=50, time_limit=time_limit)
result_3 = simulation(starting_state=50, time_limit=time_limit)
def plot_observations(observation_sequence, ax=None):
if ax is None:
ax = plt.gca()
values = [obs.state for obs in observation_sequence]
times = [obs.time_of_observation for obs in observation_sequence]
ax.plot(times, values, linestyle='steps-post')
fig, ax = plt.subplots()
plot_observations(result_1)
plot_observations(result_2)
plot_observations(result_3)
Absorbing states¶
an absorbing state is one where the transition rates for all the possible transitions is equal to zero.
When this happens, the expected residency time becomes infinite (the system is stuck there)
let's rewrite our functions so that 0 is an absorbing state
class Transition(Enum):
INCREASE = 'increase'
DECREASE = 'removal'
ABSORPION = 'absorbed'
def removal_absorbing(state):
return state*0.5
def increase_absorbing(state):
return 1 if state>0 else 0
transitions = [removal_absorbing, increase_absorbing]
transitions_names = [Transition.DECREASE, Transition.INCREASE]
state = 0
rates = [f(state) for f in transitions]
rates
total_rate = sum(rates)
total_rate
time = np.random.exponential(1/total_rate)
time
if total_rate>0:
time = np.random.exponential(1/total_rate)
else:
time = np.inf
time
at this point we should not be advancing anymore!
def simulation(starting_state, time_limit, transitions, transitions_names):
observed_states = []
state = starting_state
total_time = 0.0
while total_time < time_limit:
rates = [f(state) for f in transitions]
total_rate = sum(rates)
if total_rate>0:
time = np.random.exponential(1/total_rate)
event = rn.choices(transitions_names, weights=rates)[0]
else:
time = np.inf
event = Transition.ABSORPION
observation = Observation(state, total_time, time, event)
observed_states.append(observation)
total_time += time
if event == Transition.INCREASE:
state += 1
elif event == Transition.DECREASE:
state -= 1
elif event == Transition.ABSORPION:
pass
else:
raise ValueError("transition not recognized")
return observed_states
result = simulation(5, time_limit=100,
transitions=[removal_absorbing, increase_absorbing],
transitions_names=[Transition.DECREASE, Transition.INCREASE])
for observation in result:
print(observation)
print()
fixed time events¶
One possibility that wemight be interested in could be to register events happening at fixed point in time (or space, if we are describing the movement of an object).
the advantage of working with markovian transitions is that there is no problem interrupting at a fixed point in time and then keep going from there as nothign happened!
class Transition(Enum):
INCREASE = 'increase'
DECREASE = 'removal'
ABSORPION = 'absorbed'
RINGTONE = 'ringtone'
def removal_absorbing(state):
return state*0.5
def increase_absorbing(state):
return 1 if state>0 else 0
transitions = [removal_absorbing, increase_absorbing]
transitions_names = [Transition.DECREASE, Transition.INCREASE]
state = 5
rates = [f(state) for f in transitions]
rates
total_rate = sum(rates)
total_rate
proposed_residency_time = np.random.exponential(1/total_rate)
proposed_residency_time
current_time = 0.1
time_ringtone = 0.3
if current_time<time_ringtone and current_time + proposed_residency_time > time_ringtone:
proposed_residency_time = time_ringtone - current_time
event = Transition.RINGTONE
else:
pass
event
note that we have to decide how to interact with the absorbing states! my personal preference is to give priority to the timed events, as they might change the parameters of the simulation and change and absorbing point into a non absorbing one
def simulation(starting_state, time_limit, transitions, transitions_names):
observed_states = []
state = starting_state
total_time = 0.0
while total_time < time_limit:
rates = [f(state) for f in transitions]
total_rate = sum(rates)
if total_rate>0:
time = np.random.exponential(1/total_rate)
event = rn.choices(transitions_names, weights=rates)[0]
else:
time = np.inf
event = Transition.ABSORPION
# the fixed events overrides the others if they would happen before the chosen time
time_ringtone = 1.0
if total_time < time_ringtone and total_time + time > time_ringtone:
time = time_ringtone - total_time
event = Transition.RINGTONE
observation = Observation(state, total_time, time, event)
observed_states.append(observation)
total_time += time
if event == Transition.INCREASE:
state += 1
elif event == Transition.DECREASE:
state -= 1
elif event == Transition.ABSORPION or event == Transition.RINGTONE:
pass
else:
raise ValueError("transition not recognized")
return observed_states
result = simulation(5, time_limit=100,
transitions=[removal_absorbing, increase_absorbing],
transitions_names=[Transition.DECREASE, Transition.INCREASE])
for observation in result[:5]:
print(observation)
print()
Continuous time random walk¶
these simulations are the basics for the generalization to the continuous time random walk simulation.
These are used for example to simulate the behavior of high energy particles in matter to study the image resolution properties of certain instrumental settings, or as input for deep neural models for superresolution and so on.
In particle simulation for example the events are going to be related to a particle, that can generate new particles on its way (for example by ionization).
The fixed time events are replaced by geometrical considerations, but the underlying concept is the same.
Exercise¶
there are still many things that can be done to implement this simulator:
- allow the user to specify the effect of the various transitions
- return the simulation not as a list, but as a dedicated objct that can :
- be sliced over time:
my_simulation[0.5:100.0] - have a plot function dedicated
- automatically estimate the distribution
- be sliced over time:
- always include multiple runs of the simulation, an djoin this with point 2