Different noise realizations

[pabloriera]

Hi! Great library.

I wonder how can I get repeated simulations with differente noise realizations

I was tring this:

import numpy as np
import sdeint

M = 10

b = 1.0
tau = 0.8
tspan = np.linspace(0.0, 5.0, 5001)
x0 = 0.1*ones(M)

def f(x, t):
    return -x/tau

def g(x, t):
    return b*np.ones(x.shape)

result = sdeint.itoSRI2(f, [g]*M, x0, tspan)

But the outputs are all the same.

Thanks in advance.

[pabloriera]

Or even this

import numpy as np
import sdeint

M = 10

b = 1.0
tau = 0.8
tspan = np.linspace(0.0, 5.0, 5001)
x0 = 0.1*ones(M)

def f(x, t):
    return -x/tau

def g(x, t):
    return b*np.ones((M,1))

res = sdeint.itoSRI2(f, g, x0, tspan)

[pabloriera]

I got it, you have to create the diagonal matrix. It would be nice if there was a shortcut to create an ensemble of simulations without using an external loop.

I have a library for GPU ODE simulations, but it only has RK4. It has the option for a noisy external force but it is not a proper SDE. It runs very fast though. It would be nice to add an SDE solver, which one do you recommend for general pourpose SDEs?.

https://github.com/pabloriera/gpuODE

[mattja]

Would be great to have a GPU implementation. Yes, a ODE RK4 integrator will give the wrong answer if you just put noise in. (Simulated noise intensity will vary as the time-step varies!) The Euler-Maruyama algorithm would be a good starting place, it's very simple.

[mattja]

If you just want to do repeated simulations easily, have a look at mattja/nsim. It lets you write this:

import nsim

class Blah(nsim.ItoModel):
    b = 1.0
    tau = 0.8
    y0 = 0.1

    def f(self, y, t):
        return -y/self.tau

    def G(self, x, t):
        return self.b

sims = nsim.RepeatedSim(Blah, T=5.0, dt=0.001, repeat=10)
results = sims.timeseries
results.plot()

It will compute them all in parallel if you have a lot of CPU cores. But right now the parallel is not efficient if you are on a single machine. It is really designed for a large networked computing cluster where the RAM is not shared.

[pabloriera]

Regarding the RK4 and the noise intensity, I was normalizing the noise with the square root of the time-step (similar to Euler-Maruyama). I don't know if is accurately but gives a more or less stable noise intensity when the time-step varies. Maybe this only is valid for linear systems though.

Great, I will try nsim. I have a machine with shared RAM but I don't need long runs.

Maybe in the near future we could see to add GPU to sdeint. My module is not so clean and tidy but it has the heavier work already done.

Thanks