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 andim
-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)
- 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:
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 andn_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)