Additive noise implementation

Support
,
2

Hi @pietro. Your implementation using xi is correct, and looking at the difference between time steps is one way of checking that it works correctly (even though I don’t think you can flatten it like this, since this introduces a big jump between the last value of one trajectory, and the first value of another trajectory). But your calculation is slightly wrong, I think: you seem to expect that the standard deviation for each step is the same if normalized by dt, i.e. if one step has a standard deviation of 0.1 for a time step of 0.1ms, it should have a standard deviation of 0.2 for a time step of twice that size. But the sum of two steps with a standard variation of 0.1 each would have a total standard deviation of 0.141… not 0.2 – it’s their variance that adds up.
An alternative way to reason about (which I find more intuitive) noise is that we are modeling a Wiener process, and after a certain time, we should have the same probability distribution of values, i.e. if we simulate a trajectory multiple times and look at the value at times t, then the distribution of values at time t should be independent of the dt we used for the simulation. This is the case with the model above, see this plot (on the right in black is the histogram of the final values of the trajectory, in red is a fit of a normal distribution):
Figure_1

Printing the values confirms that everything works as expected:

dt=100. us, std dev: 0.315 (expected: 0.316)
dt=0.5 ms, std dev: 0.317 (expected: 0.316)
dt=1. ms, std dev: 0.316 (expected: 0.316)
dt=5. ms, std dev: 0.317 (expected: 0.316)

For the full code, see below. Hope this helped!

Code 
import brian2 as br
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import norm
T = 100*br.ms
sigma = 1
sigma_int = sigma * br.Hz ** 0.5
dt_values = [0.1, 0.5, 1.0, 5.0] * br.ms
results = []
for dt in dt_values:
    br.defaultclock.dt = dt 
    eqs='''
    dphi/dt = sigma_int * xi : 1
    '''
    neurons = br.NeuronGroup(10000,eqs,method='heun')
    
    noise = br.StateMonitor(neurons,'phi',record=True)
    br.run(T)
    # a detail, but the monitor records at the start of the time step,
    # so we need to record an extra step to actually have a final value
    # at the same time point, not at t-dt
    noise.record_single_timestep()
    results.append((noise.t, noise.phi[:]))

fig, axs = plt.subplots(4, 2,
                        sharex="col", sharey=True,
                        gridspec_kw={"width_ratios": [3, 1]},
                        layout="constrained")
for idx, (dt, (t, result)) in enumerate(zip(dt_values, results)):
    axs[idx, 0].set_title(f'dt = {dt}')
    axs[idx, 0].plot(t/br.ms, result[::100].T, alpha=0.1, color='black')  # plot a subset of trajectories
    axs[idx, 1].hist(result[:, -1], bins=100, orientation="horizontal", histtype="step", color="black", density=True)
    # Fit normal distribution
    mu, std = norm.fit(result[:, -1])
    print(f"dt={dt}, std dev: {std:.3f} (expected: {np.sqrt(sigma*float(T)):.3f})")
    axs[idx, 1].plot(norm.pdf(np.linspace(-1, 1, 100), mu, std), np.linspace(-1, 1, 100), color="darkred", alpha=0.5)
    axs[idx, 1].axis("off")
axs[-1, 0].set_xlabel("Time (ms)")
plt.show()