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

[Warren Weckesser]

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.

    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/bd92ad85/attachment.html>