Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
Thanks, guys. Yeah, I realized the problem w/ the uniform-increment-variable-direction approach this morning: physically, it ignores the fact that the particles hitting the particle being tracked are going to have a distribution of momentum, not all the same, just varying in direction. But I don't quite understand Warren's observation: "the 'angles' that describe the position undergo a random walk [actually, it would seem that they don't, since they too fail the varying-as-white-noise test], so the particle tends to move in the same direction over short intervals"--is this just another way of saying that, since I was varying the angles by -1, 0, or 1 unit each time, the simulation is susceptible to "unnaturally" long strings of -1, 0, or 1 increments? Thanks again,
DG
On Wed, Sep 25, 2013 at 1:41 PM, David Goldsmith d.l.goldsmith@gmail.comwrote:
Thanks, guys. Yeah, I realized the problem w/ the uniform-increment-variable-direction approach this morning: physically, it ignores the fact that the particles hitting the particle being tracked are going to have a distribution of momentum, not all the same, just varying in direction. But I don't quite understand Warren's observation: "the 'angles' that describe the position undergo a random walk [actually, it would seem that they don't, since they too fail the varying-as-white-noise test], so the particle tends to move in the same direction over short intervals"--is this just another way of saying that, since I was varying the angles by -1, 0, or 1 unit each time, the simulation is susceptible to "unnaturally" long strings of -1, 0, or 1 increments? Thanks again,
Note: I was interpreting your code as the discretization of a stochastic process, and I was experimenting with values of `incr` that were small, e.g. `incr = 0.01`.
This code
t = 2*np.pi*incr*(R.randint(3, size=(N,))-1) t[0] = 0 t = t.cumsum()
makes `t` a (discrete) random walk. At each time step, t either remains the same, or changes by +/- 2*np.pi*incr. If `incr` is small, then `t[1]` is a small step from `t[0]`. Similarly, `p[1]` will be close to `p[0]`. So the particle "remembers" its direction. A particle undergoing Brownian motion does not have this memory.
Warren
DG
NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On 25 September 2013 19:41, David Goldsmith d.l.goldsmith@gmail.com wrote:
"the 'angles' that describe the position undergo a random walk [actually, it would seem that they don't, since they too fail the varying-as-white-noise test], so the particle tends to move in the same direction over short intervals"--is this just another way of saying that, since I was varying the angles by -1, 0, or 1 unit each time, the simulation is susceptible to "unnaturally" long strings of -1, 0, or 1 increments?
In the 1D case, the white noise has a gaussian probability distribution of being positive or negative. Translated to the Wiener process, it means you would have to sum normally distributed values. When you go 3D you can do the same thing, taking a random displacement from a N(0,1) and two random angles.
The issue here is that the polar angles cannot be taken uniformly, but instead they have to be distributed proportionally to the jacobian. As you have it now, your particle will tend to move towards the poles. If you want to visualize it: take a sphere and imagine dots spaced evenly at angles (intersection of meridians and parallels, for example): they are much more dense at the poles.
The simplest way is to do it in cartesian coordinates: take x, y, and z independently from N(0,1). If you want to generate only one normal number per step, consider the jacobian in the angles.
David.
On 26 September 2013 10:02, Daπid davidmenhur@gmail.com wrote:
The simplest way is to do it in cartesian coordinates: take x, y, and z independently from N(0,1). If you want to generate only one normal number per step, consider the jacobian in the angles.
Actually, this is wrong, as it would allow displacements (at 1 sigma) of 1 along the axis, but up to sqrt(3) along diagonals. What you actually want is a multivariate normal distribution with covariance proportional to the identity (uncorrelation between axis and isotropy).
See formulae here: http://en.wikipedia.org/wiki/Multivariate_normal_distribution
On Fri, Sep 27, 2013 at 9:00 AM, Daπid davidmenhur@gmail.com wrote:
On 26 September 2013 10:02, Daπid davidmenhur@gmail.com wrote:
The simplest way is to do it in cartesian coordinates: take x, y, and z
independently from N(0,1). If you want to generate only one normal number per step, consider the jacobian in the angles.
Actually, this is wrong, as it would allow displacements (at 1 sigma) of
1 along the axis, but up to sqrt(3) along diagonals. What you actually want is a multivariate normal distribution with covariance proportional to the identity (uncorrelation between axis and isotropy).
No, you were right the first time. Sampling 3 independent N(0,1) variates is equivalent to an isotropic 3D multivariate normal. This is a special property of the normal distribution because of the behavior of exp(-x**2). The multivariate normal PDF can be decomposed into a product of univariate normals.
exp(-(x**2 + y**2 + z**2)) = exp(-x**2) * exp(-y**2) * exp(-z**2)
-- Robert Kern
participants (4)
-
David Goldsmith -
Daπid -
Robert Kern -
Warren Weckesser