function custom_sde

deeptime.data.custom_sde(dim: int, rhs: Callable, sigma: ndarray, h: float, n_steps: int)

This function allows the definition of custom stochastic differential equations (SDEs) of the form

\[\mathrm{d}X_t = F(X_t) \mathrm{d}t + \sigma(t, X_t)\mathrm{d}W_t,\]

where the right-hand side \(F\) should map an dim-dimensional array-like object to an dim-dimensional array-like object. The prefactor in front of the Wiener process \(W_t\) is assumed to be constant with respect to time and state, i.e.,

\[\sigma(t, X_t) = \sigma \in\mathbb{R}^{\mathrm{dim}\times\mathrm{dim}}.\]

(Source code, png, hires.png, pdf)

../../_images/plot_custom_sde.png
Parameters:

  • dim (int, positive and less or equal to 5) – The dimension of the SDE’s state vector \(X_t\). Must be less or equal to 5.

  • rhs (Callable) – The right-hand side function \(F(X_t)\). It must map a dim-dimensional array like object to a dim-dimensional array or list.

  • sigma ((dim, dim) ndarray) – The sigma parameter.

  • h (float) – Step size for the Euler-Maruyama integrator.

  • n_steps (int) – Number of integration steps per evaluation / recording of the state.

Returns:

system – The system.

Return type:

CustomSystem

Examples

First, some imports.

>>> import deeptime as dt
>>> import numpy as np

Then, we can define the right-hand side. Here, we choose the force of an harmonic spherical inclusion potential.

>>> def harmonic_sphere_force(x, radius=.5, k=1.):
...     dist_to_origin = np.linalg.norm(x)
...     dist_to_sphere = dist_to_origin - radius
...     if dist_to_sphere > 0:
...         return -k * dist_to_sphere * np.array(x) / dist_to_origin
...     else:
...         return [0., 0.]

This we can use as right-hand side to define our SDE with \(\sigma = \mathrm{diag}(1, 1)\).

>>> sde = dt.data.custom_sde(dim=2, rhs=lambda x: harmonic_sphere_force(x, radius=.5, k=1),
...                          sigma=np.diag([1., 1.]), h=1e-3, n_steps=1)

Here, h is the step-size of the (Euler-Maruyama) integrator and n_steps refers to the number of integration steps for each evaluation. Given the SDE instance, we can generate trajectories via

>>> trajectory = sde.trajectory(x0=[[0., 0.]], length=10, seed=55)
>>> assert trajectory.shape == (10, 2)

or propagate (in this case 300) sample points by n_steps:

>>> propagated_samples = sde(np.random.normal(scale=.1, size=(300, 2)))
>>> assert propagated_samples.shape == (300, 2)