<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Sep 25, 2013 at 12:51 PM, Warren Weckesser <span dir="ltr"><<a href="mailto:warren.weckesser@gmail.com" target="_blank">warren.weckesser@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><div class="gmail_quote"><div><div class="h5">On Wed, Sep 25, 2013 at 9:36 AM, Neal Becker <span dir="ltr"><<a href="mailto:ndbecker2@gmail.com" target="_blank">ndbecker2@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><div>David Goldsmith wrote:<br>
<br>
> Is this a valid algorithm for generating a 3D Wiener process? (When I<br>
> graph the results, they certainly look like potential Brownian motion<br>
> tracks.)<br>
><br>
> def Wiener3D(incr, N):<br>
> r = incr*(R.randint(3, size=(N,))-1)<br>
> r[0] = 0<br>
> r = r.cumsum()<br>
> t = 2*np.pi*incr*(R.randint(3, size=(N,))-1)<br>
> t[0] = 0<br>
> t = t.cumsum()<br>
> p = np.pi*incr*(R.randint(3, size=(N,))-1)<br>
> p[0] = 0<br>
> p = p.cumsum()<br>
> x = r*np.cos(t)*np.sin(p)<br>
> y = r*np.sin(t)*np.sin(p)<br>
> z = r*np.cos(p)<br>
> return np.array((x,y,z)).T<br>
><br>
> Thanks!<br>
><br>
> DG<br>
<br>
</div></div>Not the kind of Wiener process I learned of. This would be the integral of<br>
white noise. Here you have used:<br>
<br>
1. discrete increments<br>
2. spherical coordinates<br>
<br></blockquote><div><br><br></div></div></div><div>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.<br>
<br></div><div>To simulate a Wiener process (i.e. Brownian motion) in 3D, you can simply evolve each coordinate independently as a 1D process.<br><br></div><div>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.<br>
<span style="font-family:courier new,monospace"><br><br>import numpy as np<br><br><br>def wiener(x0, n, dt, delta):<br> """Generate an n-dimensional random walk.<br><br></span></div></div></div></div></blockquote>
<div><br><br></div><div>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).<br>
<br>Warren<br><br> <br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div><span style="font-family:courier new,monospace"> The array of values generated by this function simulate a Wiener process.<br>
<br> Arguments<br> ---------<br> x0 : float or array<br> The starting point of the random walk.<br> n : int<br> The number of steps to take.<br> dt : float<br> The time step.<br> delta : float<br>
delta determines the "speed" of the random walk. The random variable<br> of the position at time t, X(t), has a normal distribution whose mean<br> is the position at time t=0 and whose variance is delta**2*t.<br>
<br> Returns<br> -------<br> x : numpy array<br> The shape of `x` is (n+1,) + x0.shape.<br> The first element in the array is x0.<br> """<br> x0 = np.asfarray(x0)<br> shp = (n+1,) + x0.shape<br>
<br> # Generate a sample numbers from a normal distribution.<br> r = np.random.normal(size=shp, scale=delta*np.sqrt(dt))<br><br> # Replace the first element with 0.0, so that x0 + r.cumsum() results<br> # in the first element being x0.<br>
r[0] = 0.0<br><br> # This computes the random walk by forming the cumulative sum of<br> # the random sample. <br> x = r.cumsum(axis=0)<br> x += x0<br><br> return x<span class="HOEnZb"><font color="#888888"><br>
</font></span></span><span class="HOEnZb"><font color="#888888"><br><br><br><br></font></span></div><span class="HOEnZb"><font color="#888888"><div>Warren<br>
<br><br></div></font></span><div class="im"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org" target="_blank">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
</blockquote></div></div><br></div></div>
</blockquote></div><br></div></div>