StochRare documentation

Installation

Dependencies

• Python 3.6 or above
• numpy
• scipy
• matplotlib
• Optional: jupyter

stochrare only requires a running Python 3 install and standard scientific packages. Tutorial notebooks require Jupyter.

Basic Install

The easiest way to install stochrare is from the Python Package Index, using the Python Package Installer:

pip install stochrare

Development Mode Install

If you intend to edit the code, you should install stochrare in development mode.

First clone the repository:

git clone https://github.com/cbherbert/stochrare

Then go to the directory and install the package in development mode:

cd stochrare
pip install -e .

Tutorials

To illustrate the functionnalities and use cases of stochrare, we provide Jupyter notebooks, which can either be viewed by clicking the links below or run interactively, either by downloading your local copy of the repository, or on Binder, without installing anything: .

1D Diffusion

This tutorial illustrates some basic features of the stochrare package with simple 1D diffusion processes.

First, let us import standard packages: numpy and matplotlib.

import numpy as np
import matplotlib.pyplot as plt

Monte-Carlo simulation of 1D diffusions

One of the very common tasks stochrare can be used for is sampling trajectories of a stochastic process, or estimating directly some statistical properties of the process. We have tried to design the package so that such tasks can be easily fulfilled with an intuitive interface. The generic class for 1D diffusion processes is stochrare.dynamics.diffusion1d.DiffusionProcess1D. As we shall see later in this tutorial, any 1D diffusion process can be represented as members of this class, simply by providing the drift and diffusion functions to the constructor. But the package also ships with predefined standard processes, like the Wiener process. Let us import this process:

from stochrare.dynamics.diffusion1d import Wiener1D

Sample paths

First, we seed the random number generator with a fixed value, so that the results are always the same when we run the whole notebook. Of course, if you run a given cell multiple times, you will get different realizations of the stochastic processes.

np.random.seed(seed=100)

The basic tool to generate sample paths is the trajectory method. It integrates numerically stochastic differential equations of the form:

\[dX_t = F(X_t, t) + \sigma(X_t, t) dW_t.\]

Let us start by plotting sample paths of the Wiener process \(W\) (also known as Brownian motion) with a one-liner:

Wiener1D().trajectoryplot(*[Wiener1D().trajectory(0.,0.,T=100) for _ in range(100)]);

Or, in a different style:

from stochrare.io.plot import trajectory_plot1d
def ensemble_plot1d(*args, **kwargs):
    fig, ax = trajectory_plot1d(*((t, x, {'color': 'skyblue', 'alpha': 0.1}) for t, x in args), **kwargs)
    ax.plot(*np.array(args).mean(axis=0), color='steelblue', lw=2)
    return fig, ax
with plt.style.context(('dark_background')):
    ensemble_plot1d(*[Wiener1D().trajectory(0., 0., T=1, dt=0.01) for _ in range(200)]);

We have used two convenience methods that we provide for making quick plots: trajectory_plot1d, which is a function from the stochrare.io.plot submodule, and trajectoryplot, a method of DiffusionProcess1D objects which just provides an interface to trajectory_plot1d for quick access, and which can be overriden by subclasses to systematically include specific details in trajectory plots.

Observables and PDFs

We now estimate some statistical properties of stochastic processes.

Let us first check that for the Brownian motion the mean square displacement increases linearly with time:

ensemble = np.array([Wiener1D().trajectory(0., 0., T=10, dt=0.01) for _ in range(1000)])
time = np.average(ensemble[:, 0, :], axis=0)
ax = plt.axes(xlabel='t', ylabel='r')
ax.plot(time, np.average(ensemble[:, 1, :]**2, axis=0), label=r'$\langle \Delta x^2 \rangle$')
ax.plot(time, 2*time, label=r'$2Dt$')
ax.legend()
plt.show()

Now we estimate the probability for the stochastic process to take value \(x\) at time \(t\), knowing the initial condition \(X_{t_0}=x_0\), i.e. the transition probability \(p(x, t | x_0, t_0)\). For this, the empirical_vector method should be used: it simulates on the fly an ensemble of sample paths and returns the histogram of the values at each desired time.

Again, we work with the Wiener process, for which an analytical solution is known:

\[p(x, t | 0, 0) = \frac{1}{\sqrt{4\pi D t}}e^{-x^2/4Dt}.\]

This solution is hard-coded in the Wiener1D class, in the _fpthsol method.

Both solutions are represented using the pdf_plot1d function, another tool for quick plots dedicated to probability distributions. In the figure below, the Monte-Carlo estimate is the solid line and the theoretical result is the dotted line.

from stochrare.io.plot import pdf_plot1d
pdf = list(Wiener1D().empirical_vector(0, 0, 10000, 1, 5, bins=20))
fig, ax, lines = pdf_plot1d(*((0.5*(xx[1:]+xx[:-1]), pp, {'label': f't={t}'}) for t, pp, xx in pdf));
pdf_plot1d(*((0.5*(xx[1:]+xx[:-1]), Wiener1D()._fpthsol(0.5*(xx[1:]+xx[:-1]), t),
              {'ls': 'dotted', 'color': l.get_color()}) for l, (t, _, xx) in zip(lines, pdf)),
           fig=fig, ax=ax);

Numerical convergence

The sample paths computed with the trajectory method are discrete approximations of the sample paths of the stochastic differential equation, using the Euler-Maruyama method, which consists in computing a sequence of random numbers \(X_n\) defined by:

\[X_{n+1} = X_n + F(X_n, t_n)\Delta t + \sigma(X_n, t_n) \Delta W_n,\]

where \(t_n\) is the sample time, \(\Delta t = t_{n+1}-t_n\) the time step of the method and \(\Delta W_n = W_{t_{n+1}}-W_{t_n}\) a Gaussian random variable with zero mean and variance \(\Delta t\).

Let us now illustrate the numerical convergence of the Euler-Maruyama method. For this, we precompute the Brownian path with respect to which we integrate the SDE, and we vary the time step used for the Euler-Maruyama method.

We do this for a stochastic process which can be analytically solved:

\[dX_t = 2X_t dt + X_t dW_t,\]

and compare the numerical approximation to the analytical solution, \(X_t = X_0 e^{3t/2+W_t}\). This stochastic process is easily constructed through the constructor of the generic class DiffusionProcess1D:

from stochrare.dynamics.diffusion1d import DiffusionProcess1D
model = DiffusionProcess1D(lambda x, t: 2*x, lambda x, t: x)

dt_brownian = 0.001
brownian_path = Wiener1D(D=0.5, deterministic=True).trajectory(0., 0., T=1, dt=dt_brownian)
model.trajectoryplot((brownian_path[0], np.exp(1.5*brownian_path[0]+brownian_path[1])),
                     model.trajectory(1., 0., T=1, dt=4*dt_brownian, brownian_path=brownian_path),
                     model.trajectory(1., 0., T=1, dt=16*dt_brownian, brownian_path=brownian_path),
                    labels=('Exact', r'$\Delta t= 4 \delta t$', r'$\Delta t = 16 \delta t$'));

Slightly more precisely, let us try to illustrate that the Euler-Maruyama has strong order of convergence \(1/2\):

\[\mathbb{E}|X_n-X_{t_n}| \leq C \Delta t^{1/2}.\]

Below, we plot the error on the final point as a function of the time step used in the Euler-Maruyama method, and we compare to a linear regression and to a reference line with slope \(1/2\).

import scipy.stats

ensemble_size = 1000
dt_brownian = 0.002
dtarr = np.array((1, 2, 4, 8, 16))*dt_brownian
model = DiffusionProcess1D(lambda x, t : 2*x, lambda x, t: x)
trajs = np.array([[np.abs(model.trajectory(1., 0., T=1, dt=dt,
                                           brownian_path=brownian_path)[1][-1]-np.exp(1.5*brownian_path[0][-1]+brownian_path[1][-1]))
                   for dt in dtarr]
                  for brownian_path in (Wiener1D(D=0.5).trajectory(0., 0., T=1, dt=dt_brownian) for _ in range(ensemble_size))])
error = trajs.mean(axis=0)

ax = plt.axes(xlim=(0.001, 0.1), xlabel=r'$\Delta t$', ylabel=r'$\mathbb{E}|\hat{X}_T-X_T|$', xscale='log', yscale='log')
ax.scatter(dtarr, error)
slope, intercept, _, _, _ = scipy.stats.linregress(np.log10(dtarr), np.log10(error))
ax.plot(dtarr, 10**intercept*dtarr**slope, label=r'$\Delta t^{'+format(slope, '.2f')+'}$', color='C1')
ax.plot(dtarr, error[0]*(dtarr/dtarr[0])**0.5, label=r'$\Delta t^{1/2}$', color='C2')
ax.legend();

Numerical solution of the Fokker-Planck equation

Above we have considered stochastic processes fron the standpoint of stochastic differential equations. We now turn to an alternative point of view, the probability distribution. Markov processes are fully determined by the transition probabilities \(p(x',t'|x,t)\), which satisfy the Fokker-Planck equations:

\[\frac{\partial P}{\partial t} = - \frac{\partial}{\partial x} [FP] + \frac{1}{2}\frac{\partial^2}{\partial x^2}[\sigma^2 P].\]

Some properties of stochastic processes are much more conveniently addressed in the Fokker-Planck framework: rather than requiring the simulation of many sample paths, probability distributions can be computed directly by solving a partial differential equation.

The package includes a basic finite-difference solver for the 1D Fokker-Planck equation. We illustrate its use below.

Let us first import the submodule:

import stochrare.fokkerplanck as fp

All we need to do is create an object representing the Fokker-Planck equation (an instance of the FokkerPlanck1D class) and then use its fpintegrate method:

fpe = fp.FokkerPlanck1D(lambda x, t: 0, lambda x, t: 1)
t, X, P = fpe.fpintegrate(0, 1, dt=0.001, npts=400, bounds=(-20., 20.), P0=fpe.dirac1d(0, np.linspace(-20, 20, 400)),
                          bc=('absorbing', 'absorbing'))
pdf_plot1d((np.array(X), np.array(P)), legend=False);

Note that we have provided as arguments to fpintegrate the domain on which to solve the equation (bounds), the space resolution (npts), the time step (dt), as well as the boundary conditions (bc; here they are both absorbing boundary conditions) and an initial distribution (P0).

FokkerPlanck1D objects offer another method to compute easily the probability distribution at different times: fpintegrate_generator is a generator yielding the pdf at times given as arguments (see API documentation for more information). It relies on fpintegrate under the hood. For convenience, ConstantDiffusionProcess1D objects offer a pdfplot method which wraps the FokkerPlanck1D.fpintegrate_generator method. As an example, the solution of the heat equation at different times can be obtained as a one-liner:

Wiener1D().pdfplot(1, 5, 10, dt=0.001, npts=400, bounds=(-20.0, 20.0), t0=0.0,
                   P0=fpe.dirac1d(0, np.linspace(-20, 20, 400)), bc=('absorbing', 'absorbing'), th=True);

The theoretical solution is shown as a dotted line with the same color.

Effect of Boundary Conditions

Above we have compared the solution of the Fokker-Planck equation with absorbing boundary conditions with the theoretical solution on an infinite domain. Discrepancies should start appearing when we reach the boundaries.

Wiener1D().pdfplot(1, 10, 50, dt=0.01, npts=100, bounds=(-20.0, 20.0), t0=0.0,
                   P0=fpe.dirac1d(0, np.linspace(-20, 20, 100)), bc=('absorbing', 'absorbing'), th=True);

Below, we solve the equation with reflecting boundary conditions on both sides. In that case, we expect the system to reach a stationary state where the probability distribution is uniform on the interval.

Wiener1D().pdfplot(1, 10, 50, 200, dt=0.01, npts=100, bounds=(-20.0, 20.0), t0=0.0,
                   P0=fpe.dirac1d(0, np.linspace(-20, 20, 100)), bc=('reflecting', 'reflecting'));

Comments on numerical aspects

With the default Euler time-stepping scheme, the time step \(dt\) should be much smaller than \(dx^2/D\), with \(dx\) the spatial resolution and \(D\) the diffusivity. For instance, if we decrease the diffusivity by a factor 10, we can afford to multiply the time step by the same factor (of course the probability distribution will spread 10 times slower).

Wiener1D(D=0.1).pdfplot(10, 100, dt=0.01, npts=400, bounds=(-20.0, 20.0), t0=0.0,
                        P0=fpe.dirac1d(0, np.linspace(-20, 20, 400)), bc=('absorbing', 'absorbing'), th=True);

With the implicit and Crank-Nicolson schemes, we can afford using much larger timesteps than with the explicit method, as illustrated below:

Wiener1D(D=1).pdfplot(1, 10, dt=0.005, npts=400, bounds=(-20.0, 20.0), t0=0.0,
                      P0=fpe.dirac1d(0, np.linspace(-20, 20, 400)), bc=('absorbing', 'absorbing'), method='explicit');

Wiener1D(D=1).pdfplot(1, 10, dt=0.05, npts=400, bounds=(-20.0, 20.0), t0=0.0,
                      P0=fpe.dirac1d(0, np.linspace(-20, 20, 400)), bc=('absorbing', 'absorbing'), method='implicit', th=True);

Wiener1D(D=1).pdfplot(1, 10, dt=0.025, npts=400, bounds=(-20.0, 20.0), t0=0.0,
                      P0=fpe.dirac1d(0, np.linspace(-20, 20, 400)), bc=('absorbing', 'absorbing'), method='cn', th=True);

More simulations

In the previous tutorial we only simulated simple 1D diffusion processes. Here, we show more examples of simulations with different dynamics, illustrating in particular the generic class stochrare.dynamics.diffusion.DiffusionProcess.

2D Diffusions: gradient dynamics

Let us start with basic diffusion processes in 2D: the Wiener process and the Ornstein-Uhlenbeck process.

import numpy as np
import matplotlib.pyplot as plt
from stochrare.dynamics.diffusion import Wiener, OrnsteinUhlenbeck

np.random.seed(seed=100)

ax = plt.axes(xlabel=r'$x$', ylabel=r'$y$')
for _ in range(4):
    t, x = Wiener(2).trajectory(np.array([0., 0.]), 0., dt=0.01)
    ax.plot(x[:, 0], x[:, 1])

time, dist2 = zip(*[(t,r) for t,r in Wiener(2).sample_mean(np.array([0.,0.]), 0., 1000, 1000, dt=0.01,
                                                           observable=lambda x, t: x[0]**2+x[1]**2)])
ax = plt.axes(xlabel=r'$t$', ylabel=r'$r$')
ax.grid()
ax.plot(time, dist2, label=r'$\langle \Delta x^2 \rangle$')
ax.plot(time, 4*np.array(time), label=r'$4Dt$')
ax.legend();

model = OrnsteinUhlenbeck(0,1, 0.1, 2)
xvec = np.linspace(-1., 1.)
yvec = np.linspace(-1., 1.)
potential = np.array([model.potential(np.array([x, y])) for x in xvec for y in yvec]).reshape(50, 50)
fig = plt.figure(figsize=(5, 5))
ax = plt.axes(xlabel=r'$x$', ylabel=r'$y$')
ax.axis('equal')
ax.contourf(xvec, yvec, potential, 30, cmap='copper')
t, x = model.trajectory(np.array([0., 0.]), 0., T=2, dt=0.001)
ax.plot(x[:, 0], x[:, 1], color='white');

Langevin equation

Now we simulate the Langevin dynamics:

\[\dot{x}=p/m, \quad \dot{p}=-V'(x)-\gamma p +\eta(t),\]

with a harmonic potential \(V(x)=x^2/2\) and \(\mathbb{E}[\eta(t)\eta(t')]=D\delta(t-t')\).

from stochrare.dynamics.diffusion import DiffusionProcess
gamma = 0
def langevin(gamma, D):
    return DiffusionProcess(lambda X, t: np.array([X[1], -X[0]-gamma*X[1]]),
                            lambda X, t: np.array([[0., 0.], [0., D]]))

Without friction and noise (\(\gamma=D=0\)), the system is Hamiltonian and \(H=x^2+p^2\) is conserved. Adding some friction \(\gamma >0\), the system relaxes towards equilibrium in a spiraling motion, because of inertia. With noise, we inject energy randomly in the system, and the stationary distribution spreads over a region of phase space centered on the origin.

xvec = np.linspace(-1., 1.)
pvec = np.linspace(-1., 1.)
ax = plt.axes(xlabel=r'$x$', ylabel=r'$p$')
t, x = langevin(0, 0).trajectory(np.array([1., 0.]), 0., T=10, dt=0.001)
ax.plot(x[:,0], x[:,1]);
t, x = langevin(0.5, 0).trajectory(np.array([1., 0.]), 0., T=10, dt=0.001)
ax.plot(x[:,0], x[:,1]);
t, x = langevin(0.5, 0.2).trajectory(np.array([1., 0.]), 0., T=20, dt=0.001)
ax.plot(x[:,0], x[:,1]);

First Passage Times: the Kramers problem

This tutorial illustrates the use of the stochrare package for first-passage time computations in the context of the Kramers problem [1] (diffusion in a double well potential). Here, we simply compute numerically the rate of transition between the two attractors.

References

• [1] Kramers, Physica, 7, 284-304 (1940)
• [2] Gardiner, Handbook of Stochastic Methods, Springer, chap. 9 and §5.5.
• [3] Risken, The Fokker-Planck equation, Springer, §5.10
• [4] Caroli, Caroli and Roulet, J. Stat. Phys., 21, 415-437 (1979)
• [5] Caroli, Caroli and Roulet, J. Stat. Phys., 26, 83-111 (1981)
• [6] Hanggi, Talkner, Borkovec, Rev. Mod. Phys., 62, 251-342 (1990)

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import stochrare as sr

Let us define the DoubleWell class, corresponding to a simple bistable system, by subclassing the generic 1D diffusion class ConstantDiffusionProcess1D:

class DoubleWell(sr.dynamics.diffusion1d.ConstantDiffusionProcess1D):
    """
    Double well potential model.
    """
    default_dt = 0.01

    def __init__(self, Damp, **kwargs):
        sr.dynamics.diffusion1d.ConstantDiffusionProcess1D.__init__(self, lambda x, t: -x*(x**2-1), Damp, **kwargs)

    def potential(self, X, t):
        """
        Return the value of the potential at the input points.
        """
        Y = X**2
        return Y*(Y-2.0)/4.0

The dynamics is given by the stochastic differential equation:

\[dX_t = -V'(X_t)dt+\sqrt{2D}dW_t,\]

with \(V(x)=x^4/4-x^2/2\).

Qualitative understanding

Deterministic dynamics

Let us first look at the dynamics for the deterministic system: we plot below the phase portrait of the system, superimposed upon the potential.

ax = plt.axes(xlabel='t', ylabel='x(t)')
X = np.linspace(-1.5, 1.5, num=100)
ax.contourf(np.linspace(0, 10), X, np.tile(DoubleWell(0).potential(X, 0), (50, 1)).T, 50, cmap='copper')
for x0 in X:
    ax.plot(*DoubleWell(0).trajectory(x0, 0, T=10), color='white')

\(x=-1\) and \(x=1\) are stable fixed points with basins of attraction \(]-\infty,0[\) and \(]0,+\infty[\), respectively, separated by an unstable fixed point at \(x=0\). Typical relaxation time is of order one but diverges as \(x_0\) goes to 0.

Noise-induced transitions between the two attractors

When the diffusion coefficient does not vanish, one should expect transitions between the two attractors. Let us illustrate this below with a moderate value of the noise.

DoubleWell.trajectoryplot(DoubleWell(0.1, deterministic=True).trajectory(-1, 0, T=1000, dt=0.01));

The stationary distribution of the system should therefore be bimodal. It is simply given by \(e^{-V(x)/D}\), up to a normalization factor. Let us compute this stationary distribution numerically by solving the Fokker-Planck equation with absorbing boundary conditions at a sufficient distance.

_, ax = DoubleWell(0.1).pdfplot(0.0, 10.0, 100.0, dt=0.0005, npts=600, bounds=(-3.0, 3.0),
                                P0=sr.fokkerplanck.FokkerPlanck1D.gaussian1d(-1, 0.1, np.linspace(-3, 3, 600)));
x = np.linspace(-3, 3)
V = DoubleWell(0.1).potential(x, 0)
ax.plot(x, np.exp(-V/0.1)/np.trapz(np.exp(-V/0.1), x=x), color='black', ls='dotted', lw=2);

/Users/corentin/codes/stochrare/stochrare/fokkerplanck.py:151: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if P0 == 'gauss':
/Users/corentin/codes/stochrare/stochrare/fokkerplanck.py:153: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if P0 == 'dirac':
/Users/corentin/codes/stochrare/stochrare/fokkerplanck.py:157: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison
  if P0 == 'uniform':

Particles initially concentrated in the left well tunnel through the potential barrier, and ultimately the relative probability of the two attractors is the same (because the quasi-potential is symmetric).

To study the statistical properties of the transitions between the two attractors, we introduce the first-passage time:

\[\tau_M = \inf \{ t>0, X_t > M | x_0=-1 \}.\]

It is a random variable entirely determined by the stochastic process \(X\). In stochrare, it can be represented by the FirstPassageProcess class:

tau = sr.firstpassage.FirstPassageProcess(DoubleWell(0.1))

Methods of the FirstPassageProcess class allow for sampling the random variable, estimating its mean, moments and PDF, using direct simulation or using the Fokker-Planck equation.

tau_samples = tau.escapetime_sample(-1, 0, 0.5, ntraj=10000)

tau.escapetime_pdfplot(tau.escapetime_pdf(tau_samples))

The same result can be obtained by solving the Fokker-Planck equation for the transition probability \(P(x, t | -1, 0)\) with a reflecting boundary condition on the left side of the domain (sufficiently far) and an absorbing condition at \(x=M\). Then \(G(x_0, t)=\int_{-\infty}^M P(x, t | x_0, 0) dx\), with \(x_0=-1\) is the probability that the particle is still in the domain after time \(t\), i.e. \(P[\tau_M > t]\). Taking the opposite of the derivative with respect to time, we obtain the PDF of the first-passage time.

This computation can be done using the method firstpassagetime_cdf:

tFP, pFP = tau.firstpassagetime_cdf(-1.0, 0.5, *np.arange(10.0, 201.0, 10.0), dt=0.0001,
                                    out='pdf', npts=200, bounds=(-4.0,0.5))

In fact, the function \(G\) defined above satisfies the backwards Fokker-Planck equation. It can therefore be computed directly using the adjoint Fokker-Planck operator:

tFPb, pFPb = tau.firstpassagetime_cdf_adjoint(-1.0, 0.5, *np.arange(10.0, 201.0, 10.0), dt=0.0001,
                                              out='pdf', npts=200, bounds=(-4.0,0.5))

When the noise is small enough, it is expected that the transition times are Poisson distributed. Their statistics is entirely determined by the transition rate (the parameter of the Poisson distribution), which is the inverse of the average time between two transitions.

Let us compare the PDF of the first-passage time computed in the three ways mentioned above with the Poisson distribution, estimating the parameter from the Monte-Carlo simulation.

ax = plt.axes(xlabel=r'$\tau_M$', ylabel=r'$p(\tau_M)$', yscale='log', ylim=(5e-4, 2e-2), xlim=(0, 200))
ax.scatter(*tau.escapetime_pdf(tau_samples), label='MC');
t = np.linspace(*ax.get_xlim())
ax.plot(tFP, pFP, label='FP',marker='+',markeredgewidth=1.5, color='C1')
ax.plot(tFPb, pFPb, label='FPadjoint',marker='x',markeredgewidth=1.5, color='C2')
ax.plot(t, np.exp(-t/np.mean(tau_samples))/np.mean(tau_samples), color='C3', label='Poisson');
ax.legend();

Computing the transition rate

Let us now study what determines the transition rate, or equivalently, the mean first-passage time. For 1D homogeneous processes, a theoretical result can be obtained analytically:

\[\mathbb{E}[\tau_M] = \frac{1}{D} \int_{-1}^{M} dx e^{V(x)/D} \int_{-\infty}^x e^{-V(y)/D}dy.\]

In the small noise limit \(D \to 0\), a saddle-point approximation yields the Eyring-Kramers formula:

\[\mathbb{E}[\tau_M] \approx \frac{2\pi}{\sqrt{|V''(0)V''(-1)|}}e^{\Delta V/D},\]

for \(M> 0\) and \(\Delta V=V(0)-V(-1)=1/4\) here.

mfpt_mc = np.array([sr.firstpassage.FirstPassageProcess(DoubleWell(D)).escapetime_avg(-1, 0, 0.5, ntraj=1000)
                    for D in 1./np.arange(1, 11)])

ax = plt.axes(xlabel=r'$1/D$', ylabel=r'$\mathbb{E}[\tau_M]$', yscale='log', ylim=(1, 100))
ax.scatter(np.arange(1, 11), mfpt_mc, label='MC');
invD = np.linspace(1, 10)
ax.plot(invD, np.array([sr.firstpassage.FirstPassageProcess(DoubleWell(D)).firstpassagetime_avg_theory(-1, 0.5)
                        for D in 1./invD])[:, 1, 0], label='Theory', color='C1')
ax.plot(invD, np.sqrt(2)*np.pi*np.exp(0.25*invD), label='Eyring-Kramers', color='C2')
ax.legend();

The code also allows for computing the mean first passage time using the Fokker-Planck equation or its adjoint. We will soon update this notebook to show this method.

One could also show how the prefactor of the first-passage time depends on the threshold \(M\) (above we chose \(M=0.5\)), for a fixed noise amplitude \(D\). For \(D\) small enough, there is a sharp transition around \(M=0\), and the first-passage time depends very little on \(M\) away from this boundary layer (and it is given by the above Eyring-Kramers formula). In general, the expression of the first-passage time is always given by the integral formula before the saddle-point approximation.

Finally, these transitions can also be characterized using the instanton formalism. We shall also illustrate this in a future version of the notebook.

Return times with rare event algorithms

This tutorial demonstrates some of the features of the stochrare package for rare event simulation.

As an example, we show how to compute return times using the block maximum method (a standard method applicable to any time series) and using a particular rare event algorithm, the Adaptive Multilevel Splitting algorithm.

Let us first import the modules:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import stochrare as sr
np.random.seed(seed=100)

As an illustration, we shall work with a very simple stochastic process, the Ornstein-Uhlenbeck process in 1D:

oup = sr.dynamics.diffusion1d.OrnsteinUhlenbeck1D(0, 1, 0.5)

We simulate a realization \(x(t)\) of this process with many samples:

%%time
reftraj = oup.trajectory(0., 0., T=10000)

CPU times: user 1.31 s, sys: 25 ms, total: 1.34 s
Wall time: 1.32 s

We are interested in rare events corresponding to extreme values of the process. For instance, we represent below occurrences when the signal reaches a certain threshold \(a\).

_, ax = sr.io.plot.trajectory_plot1d((reftraj[0][:10000], reftraj[1][:10000]));
ax.axhline(y=0.7, color='C3', ls='dashed');

Such events may be characterized by their probability of occurrence per unit time. Conversely, one may consider the typical time associated with the event, for instance by looking at the average time between two successive independent events. When the events are rare enough (i.e. the threshold \(a\) is sufficiently large), the two quantities are inverse of each other. The event follow Poisson statistics, and the parameter of the Poisson distribution is the inverse of the average time between two successive independent events. This time is called the return time of the event. It is a very useful metric to quantify rare events. In this notebook, we show how to compute the return time \(r(a)\) as a function of the amplitude \(a\) of the event using different methods.

This tutorial is inspired from the paper:

Computing return times or return periods with rare event algorithms, T. Lestang, F. Ragone, C.-E. Bréhier, C. Herbert and F. Bouchet, J. Stat. Mech. (2018).

Return times with the block maximum method

The idea of the block maximum method is to divide the trajectory into \(M\) blocks of duration \(\Delta T\) (so that the total length of the trajectory is \(T_d=M\Delta T\)) much larger than the correlation time of the timeseries \(\tau_c\) (to make sure that the events are independent). On each block, the maximum value is computed:

\[a_m = \max\{ x(t) \| (m-1)\Delta T \leq t \leq m\Delta T\}.\]

To the largest of these values, we associate a return time \(T_d\), to the second largest value we associate a return time \(T_d/2\), and so on.

This procedure is implemented in the stochrare.timeseries.blockmaximum routine. In addition to the time series itself, the block size \(\Delta T\) should be given to the routine. It should be chosen such that \(\tau_c \ll \Delta T \ll r(a)\).

%%time
aref, rref = zip(*sr.timeseries.blockmaximum(reftraj[1], 1000, mode='returntime', time=reftraj[0], modified=True))

CPU times: user 12 ms, sys: 5.31 ms, total: 17.3 ms
Wall time: 13.4 ms

ax = plt.axes(xlabel=r'$r(a)$', ylabel=r'$a$', xscale='log')
ax.plot(rref, aref);

Return times with the Adaptive Multilevel Splitting algorithm

As explained in Lestang et al. (2019), the method can be extended to non-equiprobable blocks like trajectories simulated by rare event algorithms.

As an example, we use the Trajectory Adaptive Multilevel Splitting algorithm, which is defined by a score function (below it is just the identity) and a fixed length for trajectories (equal to 5 below). The idea of the algorithm is to evolve an ensemble of trajectories through selection-mutation steps. At each iteration, the poorest performers of the ensemble (measured by the score function) are killed and replaced by a copy of a surviving trajectory resampled from the point where it reached the maximum score function level obtained by the killed trajectory.

We first define the TAMS object, which requires a dynamics, a score function and the duration for trajectories:

tams = sr.rare.ams.TAMS(oup, (lambda t, x: x), 5.)

Then, we run the algorithm. Here we use directly the returntimes method, which samples trajectories by running the algorithm and then computes the corresponding return times. The method takes as arguments the number of member of the initial ensemble (here 100), and the number of iterations for the algorithm (here 600).

%%time
aams, rams = tams.returntimes(100, 600)

CPU times: user 726 ms, sys: 6.14 ms, total: 732 ms
Wall time: 732 ms

Let us compare the solution to the one obtained with the Block Maximum method:

ax = plt.axes(xlabel=r'$r(a)$', ylabel=r'$a$', xscale='log')
ax.plot(rref, aref, label='Block Maximum');
ax.plot(rams, aams, label='TAMS');
ax.legend();

In this simple example, for a similar computational cost, the AMS algorithm allows to estimate return times 7 orders of magnitude larger than the Block Maximum method. This depends on the realization, and to properly characterize the performance of the AMS algorithm one would need to study the statistics over an ensemble of realizations.

Applications of the rare event algorithms

For more information about application of the algorithm shown here to sample efficiently rare events, we refer the reader to the following articles:

• 6. Cérou, A. Guyader, T. Lelièvre and D. Pommier, J. Chem. Phys. 134, 054108 (2011).
• 10. Rolland, F. Bouchet and E. Simonnet, J. Stat. Phys. 162, 277–311 (2015).
• 6. Ragone, J. Wouters and F. Bouchet, Proc. Nat. Acad. Sci. 115, 24-29 (2018).
• 20. Lestang, F. Ragone, C.-E. Bréhier, C. Herbert and F. Bouchet, J. Stat. Mech. (2018).
• 10. Rolland, Phys. Rev. E 97, 023109 (2018).
• 6. Bouchet, J. Rolland and E. Simonnet, Phys. Rev. Lett. 122, 074502 (2019).

In addition, a recent review of the AMS algorithm can be found here:

• 6. Cérou, A. Guyader, and M. Rousset, Chaos 29, 043108 (2019).

API Reference

stochrare is organized in several modules serving different purposes:

stochrare.dynamics	Sample stochastic processes
stochrare.fokkerplanck	Numerical solvers for the Fokker-Planck equations
stochrare.firstpassage	First-passage processes
stochrare.rare	Rare event algorithms
stochrare.io	Input/Output
stochrare.utils	Utilities