Extract noise?

[laurito2]

I am using sdeint as it is shown in the 2nd example. So multi-dimensional independent driving Wiener processes. Is it possible to save the noise for every dimension/timestep?

[superconvergence]

It is possible to use directly the underlying integrator (sdeint.itoSRI2) which offers an option to explicitly specify Wiener increments (argument dW). You can generate the increments using e.g. numpy.random.randn and from those it is straightforward to calculate the noise.

[mattja]

Yes, superconvergence is right. You pre-generate the noise, then pass it in to itoSRI2() or itoEuler(). You can use the function sdeint.deltaW(n, m, h) to pre-generate the noise increments. (It just calls numpy.random.normal internally)

But please note that if you want to save enough information so you can replay the same realization of the Wiener process, then storing the increments dW is not enough.

You also need to store a realization of the repeated integrals I This is because, even if you know that the Wiener process went from point A to point B duing the short time h, there are many different paths it could have taken from A to B during that short time. The repeated integrals will be different depending on which path was taken between time steps in a specific realization.

These repeated integrals I can be pre-generated and stored too, by calling the function sdeint.Ikpw(dW, h) with your pre-generated noise increments. (or similarly Jkpw in the case of a Stratonovich system).

So if you want to be able to replay the noise process exactly, the example code becomes:

import numpy as np
import sdeint

A = np.array([[-0.5, -2.0],
              [ 2.0, -1.0]])

B = np.diag([0.5, 0.5]) # diagonal, so independent driving Wiener processes

x0 = np.array([3.0, 3.0])

def f(x, t):
    return A.dot(x)

def G(x, t):
    return B

tspan = np.linspace(0.0, 10.0, 10001)
dW = sdeint.deltaW(10000, 2, 0.001)   # pre-generate Wiener increments
_, I = sdeint.Ikpw(dW, 0.001)  # pre-generate repeated integrals

result = sdeint.itoSRI2(f, G, x0, tspan, dW=dW, I=I)

Then you can save dW and I and re-use them. You don't need to worry about I if you are using the simple algorithms itoEuler() or stratHeun(). For those, saving the noise increments dW is enough. But any higher-order algorithm will use a realization of repeated integrals I. Hope this helps.