[Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?

[Warren Weckesser]

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] = 0
>> >     r = r.cumsum()
>> >     t = 2*np.pi*incr*(R.randint(3, size=(N,))-1)
>> >     t[0] = 0
>> >     t = t.cumsum()
>> >     p = np.pi*incr*(R.randint(3, size=(N,))-1)
>> >     p[0] = 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
>> of
>> 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).

Warren



>     The array of values generated by this function simulate a Wiener
> process.
>
>     Arguments
>     ---------
>     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
> variable
>         of the position at time t, X(t), has a normal distribution whose
> mean
>         is the position at time t=0 and whose variance is delta**2*t.
>
>     Returns
>     -------
>     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.0
>
>     # This computes the random walk by forming the cumulative sum of
>     # the random sample.
>     x = r.cumsum(axis=0)
>     x += x0
>
>     return x
>
>
>
>
> Warren
>
>
> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20130925/1a58b834/attachment.html>