On Wed, Sep 25, 2013 at 12:51 PM, Warren Weckesser < warren.weckes...@gmail.com> wrote:
> > On Wed, Sep 25, 2013 at 9:36 AM, Neal Becker <ndbeck...@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 >> > >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion