Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Extract noise? #10

Open
laurito2 opened this issue Apr 4, 2017 · 2 comments
Open

Extract noise? #10

laurito2 opened this issue Apr 4, 2017 · 2 comments
Labels

Comments

@laurito2
Copy link

laurito2 commented Apr 4, 2017

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
Copy link

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 mattja added the question label Apr 6, 2017
@mattja
Copy link
Owner

mattja commented Apr 7, 2017

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 i_definition_trans
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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

3 participants