stochrare.dynamics.diffusion1d¶

Simulating 1D diffusion processes¶

This module defines the DiffusionProcess1D class, representing diffusion processes with arbitrary drift and diffusion coefficients in 1D.

This class can be subclassed for specific diffusion processes for which methods can be specialized, both to simplify the code (e.g. directly enter analytical formulae when they are available) and for performance. As an exemple of this mechanism, we also provide in this module the ConstantDiffusionProcess1D class, for which the diffusion term is constant (additive noise), the OrnsteinUhlenbeck1D class representing the particular case of the Ornstein-Uhlenbeck process, and the Wiener1D class corresponding to Brownian motion. These classes form a hierarchy deriving from the base class, DiffusionProcess1D.

class stochrare.dynamics.diffusion1d.DiffusionProcess1D(vecfield, sigma, **kwargs)¶

Bases: object

Generic class for 1D diffusion processes.

It corresponds to the family of 1D SDEs \(dx_t = F(x_t, t)dt + \sigma(x_t, t)dW_t\), where \(F\) is a time-dependent vector field and \(W\) the Wiener process.

Parameters:	vecfield (function with two arguments) – The vector field \(F(x, t)\). sigma (function with two arguments) – The diffusion coefficient \(\sigma(x, t)\).

potential(X, t)¶

Compute the potential from which the force derives.

Parameters:	X (ndarray) – The points where we want to compute the potential.
Returns:	V – The potential from which the force derives, at the given points.
Return type:	ndarray

Notes

We integrate the vector field to obtain the value of the underlying potential at the input points. Caveat: This works only for 1D dynamics.

update(xn, tn, **kwargs)¶

Return the next sample for the time-discretized process.

Keyword Arguments:
Parameters:	xn (float) – The current position. tn (float) – The current time.
	dt (float) – The time step (default 0.1 if not overriden by a subclass). dw (float) – The brownian increment if precomputed. By default, it is generated on the fly from a Gaussian distribution with variance \(dt\).
Returns:	x – The position at time tn+dt.
Return type:	float

Notes

This method uses the Euler-Maruyama method [1] [2]: \(x_{n+1} = x_n + F(x_n, t_n)\Delta t + \sigma(x_n, t_n) \Delta W_n\).

It is the straightforward generalization to SDEs of the Euler method for ODEs.

The Euler-Maruyama method has strong order 0.5 and weak order 1.

References

[1]	G. Maruyama, Continuous Markov processes and stochastic equations, Rend. Circ. Mat. Palermo 4, 48-90 (1955).

[2]	P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer (1992).

trajectory(x0, t0, **kwargs)¶

Integrate the SDE with given initial condition.

Keyword Arguments:
Parameters:	x0 (float) – The initial position. t0 (float) – The initial time.
	dt (float) – The time step, forwarded to the `update()` routine (default 0.1, unless overridden by a subclass). T (float) – The time duration of the trajectory (default 10). brownian_path ((ndarray, ndarray)) – A precomputed Brownian path with respect to which we integrate the SDE. If not provided (default behavior), one will be computed one the fly. deltat (float) – The time step for the Brownian path, when generated on the fly (default: dt). finite (bool) – Filter finite values before returning trajectory (default False).
Returns:	t, x – Time-discrete sample path for the stochastic process with initial conditions (t0, x0). The array t contains the time discretization and x the value of the sample path at these instants.
Return type:	ndarray, ndarray

trajectory_conditional(x0, t0, pred, **kwargs)¶

Compute sample path satisfying arbitrary condition.

Keyword Arguments:
Parameters:	x0 (float) – The initial position. t0 (float) – The initial time. pred (function with two arguments) – The predicate to select trajectories.
	dt (float) – The time step, forwarded to the `update()` routine (default 0.1, unless overridden by a subclass). T (float) – The time duration of the trajectory (default 10). finite (bool) – Filter finite values before returning trajectory (default False).
Returns:	t, x – Time-discrete sample path for the stochastic process with initial conditions (t0, x0). The array t contains the time discretization and x the value of the sample path at these instants.
Return type:	ndarray, ndarray

empirical_vector(x0, t0, nsamples, *args, **kwargs)¶

Empirical vector at given times.

Keyword Arguments:
Parameters:	x0 (float) – Initial position. t0 (float) – Initial time. nsamples (int) – The size of the ensemble. args (variable length argument list*) – The times at which we want to estimate the empirical vector.
	**kwargs – Keyword arguments forwarded to `trajectory()` and to `numpy.histogram()`.
Yields:	t, pdf, bins (float, ndarray, ndarray) – The time and histogram of the stochastic process at that time.

Notes

This method computes the empirical vector, or in other words, the relative frequency of the stochastic process at different times, conditioned on the initial condition. At each time, the empirical vector is a random vector. It is an estimator of the transition probability \(p(x, t | x_0, t_0)\).

classmethod trajectoryplot(*args, **kwargs)¶

Plot 1D trajectories.

Keyword Arguments:
Parameters:	args (variable length argument list) – trajs* (tuple (t, x)) –
	fig (matplotlig.figure.Figure) – Figure object to use for the plot. Create one if not provided. ax (matplotlig.axes.Axes) – Axes object to use for the plot. Create one if not provided. **kwargs – Other keyword arguments forwarded to matplotlib.pyplot.axes.
Returns:	fig, ax – The figure.
Return type:	matplotlib.figure.Figure, matplotlib.axes.Axes

Notes

This is just an interface to the function stochrare.io.plot.trajectory_plot1d(). However, it may be overwritten in subclasses to systematically include elements to the plot which are specific to the stochastic process.

class stochrare.dynamics.diffusion1d.ConstantDiffusionProcess1D(vecfield, Damp, **kwargs)¶

Bases: stochrare.dynamics.diffusion1d.DiffusionProcess1D

Diffusion processes in 1D with constant diffusion coefficient (additive noise).

It corresponds to the family of SDEs \(dx_t = F(x_t, t)dt + \sigma dW_t\), where \(F\) is a time-dependent vector field and \(W\) the Wiener process. The diffusion coefficient \(\sigma\) is independent of space and time.

Parameters:	vecfield (function with two arguments) – The vector field \(F(x, t)\). Damp (float) – The amplitude of the noise.

Notes

The diffusion coefficient is given by \(\sigma=\sqrt{2\text{Damp}}\). This convention leads to simpler expressions, for instance for the Fokker-Planck equations.

update(xn, tn, **kwargs)¶

Return the next sample for the time-discretized process.

Keyword Arguments:
Parameters:	xn (float) – The current position. tn (float) – The current time.
	dt (float) – The time step (default 0.1 if not overriden by a subclass). dw (float) – The brownian increment if precomputed. By default, it is generated on the fly from a Gaussian distribution with variance \(dt\).
Returns:	x – The position at time tn+dt.
Return type:	float

Notes

This method uses the Euler-Maruyama method [3] [4]: \(x_{n+1} = x_n + F(x_n, t_n)\Delta t + \sqrt{2D} \Delta W_n\).

References

[3]	G. Maruyama, Continuous Markov processes and stochastic equations, Rend. Circ. Mat. Palermo 4, 48-90 (1955).

[4]	P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer (1992).

traj_cond_gen(x0, t0, tau, M, **kwargs)¶

Generate trajectories conditioned on the first-passage time tau at value M.

Keyword Arguments:
Parameters:	x0 (float) – Initial position. t0 (float) – Initial time. tau (float) – The value of the first passage time required. M (float) – The threshold for the first passage time.
	dt (float) – The integration timestep (default is self.default_dt). ttol (float) – The first-passage time tolerance (default is 1% of trajectory duration). num (int) – The number of trajectories generated (default is 10). interp (bool) – Interpolate to generate unifomly sampled trajectories. npts (int) – The number of points for interpolated trajectories (default (tau-t0)/dt).
Yields:	t, x (ndarray, ndarray) – Trajectories satisfying the condition on the first passage time.

pdfplot(*args, **kwargs)¶

Plot the pdf P(x,t) at various times.

Keyword Arguments:
Parameters:	args (variable length argument list) – The times at which to plot the PDF.
	t0 (float) – Initial time. potential (bool) – Plot potential on top of PDF. th (bool) – Plot theoretical solution, if it exists, on top of PDF.

action(*args)¶: Compute the action for all the trajectories given as arguments

class stochrare.dynamics.diffusion1d.OrnsteinUhlenbeck1D(mu, theta, D, **kwargs)¶

Bases: stochrare.dynamics.diffusion1d.ConstantDiffusionProcess1D

The 1D Ornstein-Uhlenbeck process.

It corresponds to the SDE \(dx_t = \theta(\mu-x_t)dt + \sqrt{2D} dW_t\), where \(\theta>0\) and \(\mu\) are arbitrary coefficients and \(D>0\) is the amplitude of the noise.

Parameters:	mu (float) – The expectation value. theta (float) – The inverse of the relaxation time. D (float) – The amplitude of the noise.

Notes

The Ornstein-Uhlenbeck process has been used to model many systems. It was initially introduced to describe the motion of a massive Brownian particle with friction [5] . It may also be seen as a diffusion process in a harmonic potential.

Because many of its properties can be computed analytically, it provides a useful toy model for developing new methods.

References

[5]	G. E. Uhlenbeck and L. S. Ornstein, “On the theory of Brownian Motion”. Phys. Rev. 36, 823–841 (1930).

update(xn, tn, **kwargs)¶

Return the next sample for the time-discretized process, using the Gillespie method.

Keyword Arguments:
Parameters:	xn (float) – The current position. tn (float) – The current time.
	dt (float) – The time step (default 0.1 if not overriden by a subclass). dw (float) – The brownian increment if precomputed. By default, it is generated on the fly from a standard Gaussian distribution. method (str) – The numerical method for integration: ‘gillespie’ (default) or ‘euler’.
Returns:	x – The position at time tn+dt.
Return type:	float

Notes

For the Ornstein-Uhlenbeck process, there is an exact method, the Gillespie algorithm [6]. This method is selected by default. If necessary, the Euler-Maruyama method can still be chosen using the method keyword argument.

References

[6]	D. T. Gillespie, Exact numerical simulation of the Ornstein-Uhlenbeck process and its integral, Phys. Rev. E 54, 2084 (1996).

mean_firstpassage_time(x0, a)¶

Return the mean first-passage time for the 1D Ornstein-Uhlenbeck process (exact formula).

Parameters:	x0 (float) – Initial position a (float) – Threshold
Returns:	tau – Mean first-passage time
Return type:	float

Notes

The first passage time is defined by \(\tau_a(x_0)=\inf \{t>0, X_t>a | X_0=x_0\}\). It is a random variable. Here, we compute only its expectation value, for which an analytical formula is known.

General methods for first-passage time conputations are avaiblable in the stochrare.firstpassage module.

class stochrare.dynamics.diffusion1d.Wiener1D(D=1, **kwargs)¶

Bases: stochrare.dynamics.diffusion1d.OrnsteinUhlenbeck1D

The 1D Wiener process.

Parameters:	D (float, optional) – The amplitude of the noise (default is 1).

Notes

The Wiener process is a central object in the theory or stochastic processes, both from a mathematical point of view and for its applications in different scientific fields. We refer to classical textbooks for more information about the Wiener process and Brownian motion.

classmethod potential(X)¶

Compute the potential from which the force derives.

Parameters:	X (ndarray) – The points where we want to compute the potential.
Returns:	V – The potential from which the force derives, at the given points.
Return type:	float

Notes

The Wiener Process is a trivial gradient system, with vanishing potential. It is useless (and potentially source of errors) to call the general potential routine, so we just return zero directly.

Classes

`ConstantDiffusionProcess1D`(vecfield, Damp, …)	Diffusion processes in 1D with constant diffusion coefficient (additive noise).
`DiffusionProcess1D`(vecfield, sigma, **kwargs)	Generic class for 1D diffusion processes.
`DrivenOrnsteinUhlenbeck1D`(mu, theta, D, A, …)	The 1D Ornstein-Uhlenbeck model driven by a periodic forcing:
`OrnsteinUhlenbeck1D`(mu, theta, D, **kwargs)	The 1D Ornstein-Uhlenbeck process.
`Wiener1D`([D])	The 1D Wiener process.