On Wed, Sep 25, 2013 at 12:51 PM, Warren Weckesser < warren.weckesser@gmail.com> wrote:

On Wed, Sep 25, 2013 at 9:36 AM, Neal Becker <ndbecker2@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion