Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
On 26 September 2013 10:02, Daπid davidmen...@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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
On Fri, Sep 27, 2013 at 9:00 AM, Daπid davidmen...@gmail.com wrote: On 26 September 2013 10:02, Daπid davidmen...@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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
On 25 September 2013 19:41, David Goldsmith d.l.goldsm...@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. ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
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. 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
Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
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
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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
On Wed, Sep 25, 2013 at 1:41 PM, David Goldsmith d.l.goldsm...@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 ___ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
