[Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
warren.weckesser at gmail.com
Wed Sep 25 13:39:03 EDT 2013
On Wed, Sep 25, 2013 at 12:51 PM, Warren Weckesser <
warren.weckesser at gmail.com> wrote:
> On Wed, Sep 25, 2013 at 9:36 AM, Neal Becker <ndbecker2 at gmail.com> wrote:
>> David Goldsmith wrote:
>> > Is this a valid algorithm for generating a 3D Wiener process? (When I
>> > graph the results, they certainly look like potential Brownian motion
>> > tracks.)
>> > def Wiener3D(incr, N):
>> > r = incr*(R.randint(3, size=(N,))-1)
>> > r = 0
>> > r = r.cumsum()
>> > t = 2*np.pi*incr*(R.randint(3, size=(N,))-1)
>> > t = 0
>> > t = t.cumsum()
>> > p = np.pi*incr*(R.randint(3, size=(N,))-1)
>> > p = 0
>> > p = p.cumsum()
>> > x = r*np.cos(t)*np.sin(p)
>> > y = r*np.sin(t)*np.sin(p)
>> > z = r*np.cos(p)
>> > return np.array((x,y,z)).T
>> > Thanks!
>> > DG
>> Not the kind of Wiener process I learned of. This would be the integral
>> white noise. Here you have used:
>> 1. discrete increments
>> 2. spherical coordinates
> I agree with Neal: that is not a Wiener process. In your process, the
> *angles* that describe the position undergo a random walk, so the particle
> tends to move in the same direction over short intervals.
> To simulate a Wiener process (i.e. Brownian motion) in 3D, you can simply
> evolve each coordinate independently as a 1D process.
> Here's a simple function to generate a sample from a Wiener process. The
> dimension is determined by the shape of the starting point x0.
> import numpy as np
> def wiener(x0, n, dt, delta):
> """Generate an n-dimensional random walk.
Whoops--that's a misleading docstring. The `n` in "an n-dimensional
random walk" is not the same `n` that is the second argument of the
function (which is the number of steps to compute).
> The array of values generated by this function simulate a Wiener
> x0 : float or array
> The starting point of the random walk.
> n : int
> The number of steps to take.
> dt : float
> The time step.
> delta : float
> delta determines the "speed" of the random walk. The random
> of the position at time t, X(t), has a normal distribution whose
> is the position at time t=0 and whose variance is delta**2*t.
> x : numpy array
> The shape of `x` is (n+1,) + x0.shape.
> The first element in the array is x0.
> x0 = np.asfarray(x0)
> shp = (n+1,) + x0.shape
> # Generate a sample numbers from a normal distribution.
> r = np.random.normal(size=shp, scale=delta*np.sqrt(dt))
> # Replace the first element with 0.0, so that x0 + r.cumsum() results
> # in the first element being x0.
> r = 0.0
> # This computes the random walk by forming the cumulative sum of
> # the random sample.
> x = r.cumsum(axis=0)
> x += x0
> return x
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion