Overdamped

import numpy as np
import folie as fl
import matplotlib.pyplot as plt

Let’s first simulate some trajectories with various timestep

model_simu = fl.models.OrnsteinUhlenbeck()
model_simu.coefficients = np.array([0.1, 1.2, 2.0])
data_dts = []
list_dts = [1e-3, 5e-3, 1e-2, 5e-2]
for dt in list_dts:
    simulator = fl.Simulator(fl.simulations.ExactStepper(model_simu), dt)
    data_dts.append(simulator.run(5000, np.zeros((25,)), 25))

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[2], line 7
      5 for dt in list_dts:
      6     simulator = fl.Simulator(fl.ExactDensity(model_simu), dt)
----> 7     data_dts.append(simulator.run(5000, np.zeros((25,)), 25))

File ~/Projets/folie/folie/simulations/__init__.py:34, in Simulator.run(self, nsteps, x0, ntrajs, save_every, **kwargs)
     32     x = self.transition.run_step(x, self.dt, dW[:, n])
     33     if n % save_every == 0:
---> 34         x_val[:, n // save_every, :] = x
     35 data = Trajectories(dt=self.dt * save_every)
     36 for i in range(ntrajs):

ValueError: could not broadcast input array from shape (25,25) into shape (25,1)

We can then run the estimation for various likelihood at all timesteps

fig, axs = plt.subplots(1, len(model_simu.coefficients))

for name, transitioncls in zip(
    ["Exact", "Euler", "Elerian", "Kessler", "Drozdov"],
    [
        fl.ExactDensity,
        fl.EulerDensity,
        fl.ElerianDensity,
        fl.KesslerDensity,
        fl.DrozdovDensity,
    ],
):
    model = fl.models.OrnsteinUhlenbeck()
    estimator = fl.LikelihoodEstimator(transitioncls(model))
    coeffs_vals = np.empty((len(data_dts), len(model.coefficients)))
    for n, data in enumerate(data_dts):
        res = estimator.fit_fetch(data)
        coeffs_vals[n, :] = res.coefficients
    for n in range(len(axs)):
        axs[n].plot(list_dts, np.abs(coeffs_vals[:, n] - model_simu.coefficients[n]), "-+", label=name)
    print(coeffs_vals)
for n in range(len(axs)):
    axs[n].legend()
    axs[n].set_yscale("log")
    axs[n].grid()
    axs[n].set_xlabel(r"$\Delta t$")
    axs[n].set_ylabel(r"$|c-c_{real}|$")
plt.show()