stochrare.dynamics.diffusion¶

Simulating diffusion processes in arbitrary dimensions¶

This module defines the DiffusionProcess class, representing generic diffusion processes with arbitrary drift and diffusion coefficients, in arbitrary dimension.

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 ConstantDiffusionProcess class, for which the diffusion term is constant and proportional to the identity matrix, the OrnsteinUhlenbeck class representing the particular case of the Ornstein-Uhlenbeck process, and the Wiener class corresponding to Brownian motion. These classes form a hierarchy deriving from the base class, DiffusionProcess.

class stochrare.dynamics.diffusion.DiffusionProcess(vecfield, sigma, **kwargs)¶

Bases: object

Generic class for diffusion processes in arbitrary dimensions.

It corresponds to the family of SDEs \(dx_t = F(x_t, t)dt + \sigma(x_t, t)dW_t\), where \(F\) is a time-dependent \(N\)-dimensional vector field and \(W\) the \(M\)-dimensional Wiener process. The diffusion matrix sigma has size NxM.

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

update(xn, tn, **kwargs)¶

Return the next sample for the time-discretized process.

Keyword Arguments:
Parameters:	xn (ndarray) – A n-dimensional vector (in \(\mathbb{R}^n\)). tn (float) – The current time.
	dt (float) – The time step. dw (ndarray) – 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:	ndarray

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\), for a fixed time step \(\Delta t\), where \(\Delta W_n\) is a random vector distributed according to the standard normal distribution [1] [2].

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]	(1, 2) G. Maruyama, “Continuous Markov processes and stochastic equations”, Rend. Circ. Mat. Palermo 4, 48-90 (1955).

[2]	(1, 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 (ndarray) – The initial position (in \(\mathbb{R}^n\)). 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). 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_generator(x0, t0, nsteps, **kwargs)¶

Integrate the SDE with given initial condition, generator version.

Keyword Arguments:
Parameters:	x0 (ndarray) – The initial position (in \(\mathbb{R}^n\)). t0 (float) – The initial time. nsteps (int) – The number of samples to generate.
	dt (float) – The time step, forwarded to the `update()` routine (default 0.1, unless overridden by a subclass). observable (function with two arguments) – Time-dependent observable \(O(x, t)\) to compute (default \(O(x, t)=x\))
Yields:	t, y (ndarray, ndarray) – Time-discrete sample path (or observable) for the stochastic process with initial conditions (t0, x0). The array t contains the time discretization and y=O(x, t) the value of the observable (it may be the stochastic process itself) at these instants.

sample_mean(x0, t0, nsteps, nsamples, **kwargs)¶

Compute the sample mean of a time dependent observable, conditioned on initial conditions.

Keyword Arguments:
Parameters:	x0 (ndarray) – The initial position (in \(\mathbb{R}^n\)). t0 (float) – The initial time. nsteps (int) – The number of samples in each sample path. nsamples (int) – The number of sample paths in the ensemble.
	dt (float) – The time step, forwarded to the `update()` routine (default 0.1, unless overridden by a subclass). observable (function with two arguments) – Time-dependent observable \(O(x, t)\) to compute (default \(O(x, t)=x\))
Yields:	t, y (ndarray, ndarray) – Time-discrete ensemble mean for the observable, conditioned on the initial conditions (t0, x0). The array t contains the time discretization and \(y=\mathbb{E}[O(x, t)]\) the value of the sample mean of the observable (it may be the stochastic process itself) at these instants.

class stochrare.dynamics.diffusion.ConstantDiffusionProcess(vecfield, Damp, dim, **kwargs)¶

Bases: stochrare.dynamics.diffusion.DiffusionProcess

Diffusion processes, in arbitrary dimensions, with constant diffusion coefficient.

It corresponds to the family of SDEs \(dx_t = F(x_t, t)dt + \sigma dW_t\), where \(F\) is a time-dependent \(N\)-dimensional vector field and \(W\) the \(N\)-dimensional Wiener process. The diffusion coefficient \(\sigma\) is independent of the stochastic process (additive noise) and time, and we further assume that it is proportional to the identity matrix: all the components of the noise are independent.

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

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 (ndarray) – A n-dimensional vector (in \(\mathbb{R}^n\)). tn (float) – The current time.
	dt (float) – The time step. dw (ndarray) – 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:	ndarray

See also

DiffusionProcess.update(): for details about the Euler-Maruyama method.

Notes

This is the same as the DiffusionProcess.update() method from the parent class DiffusionProcess, except that a matrix product is no longer necessary.

class stochrare.dynamics.diffusion.OrnsteinUhlenbeck(mu, theta, D, dim, **kwargs)¶

Bases: stochrare.dynamics.diffusion.ConstantDiffusionProcess

The Ornstein-Uhlenbeck process, in arbitrary dimensions.

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

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

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 [3] . 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

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

potential(x)¶

Compute the potential from which the force derives.

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

Notes

Not all diffusion processes derive from a potential, but the Ornstein Uhlenbeck does. It is a gradient system, with a quadratic potential: \(dx_t = -\nabla V(x_t)dt + \sqrt{2D} dW_t\), with \(V(x) = \theta(\mu-x)^2/2\).

class stochrare.dynamics.diffusion.Wiener(dim, D=1, **kwargs)¶

Bases: stochrare.dynamics.diffusion.OrnsteinUhlenbeck

The Wiener process, in arbitrary dimensions.

Parameters:	dim (int) – The dimension of the system. 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 point where we want to compute the potential.
Returns:	V – The potential from which the force derives, at the given point.
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

`ConstantDiffusionProcess`(vecfield, Damp, …)	Diffusion processes, in arbitrary dimensions, with constant diffusion coefficient.
`DiffusionProcess`(vecfield, sigma, **kwargs)	Generic class for diffusion processes in arbitrary dimensions.
`OrnsteinUhlenbeck`(mu, theta, D, dim, **kwargs)	The Ornstein-Uhlenbeck process, in arbitrary dimensions.
`Wiener`(dim[, D])	The Wiener process, in arbitrary dimensions.