布朗运动¶
| 日期 | 2017-03-03（最后修改），2010-01-17（创建） |
|---|
布朗运动是一个随机过程。布朗运动方程的一种形式是
$X(0) = X_0$
$X(t + dt) = X(t) + N(0, (delta)^2 dt; t, t+dt)$
其中 $N(a, b; t_1, t_2)$ 是一个均值为 a，方差为 b 的正态分布随机变量。参数 $t_1$ 和 $t_2$ 明确了N在不同时间间隔上的统计独立性；也就是说，如果 $[t_1, t_2)$ 和 $[t_3, t_4)$ 是不相交的间隔，那么 $N(a, b; t_1, t_2)$ 和 $N(a, b; t_3, t_4)$ 是独立的。
计算实际上非常简单。一个打印n步布朗运动的朴素实现可能如下所示
In [1]
from scipy.stats import norm
# Process parameters
delta = 0.25
dt = 0.1
# Initial condition.
x = 0.0
# Number of iterations to compute.
n = 20
# Iterate to compute the steps of the Brownian motion.
for k in range(n):
x = x + norm.rvs(scale=delta**2*dt)
print x
上面的代码可以很容易地修改为将迭代保存到数组中，而不是打印它们。
上面代码的问题是它很慢。如果我们想要计算大量的迭代，我们可以做得更好。关键是注意到计算是来自正态分布的样本的累积和。可以通过首先使用 scipy.stats.norm.rvs() 一次性生成来自正态分布的所有样本，然后使用 numpy 的cumsum函数来形成累积和，来实现一个快速版本。
以下函数使用这个想法来实现brownian()函数。该函数允许初始条件为数组（或任何可以转换为数组的东西）。x0的每个元素都被视为布朗运动的初始条件。
In [3]
"""
brownian() implements one dimensional Brownian motion (i.e. the Wiener process).
"""
# File: brownian.py
from math import sqrt
from scipy.stats import norm
import numpy as np
def brownian(x0, n, dt, delta, out=None):
"""
Generate an instance of Brownian motion (i.e. the Wiener process):
X(t) = X(0) + N(0, delta**2 * t; 0, t)
where N(a,b; t0, t1) is a normally distributed random variable with mean a and
variance b. The parameters t0 and t1 make explicit the statistical
independence of N on different time intervals; that is, if [t0, t1) and
[t2, t3) are disjoint intervals, then N(a, b; t0, t1) and N(a, b; t2, t3)
are independent.
Written as an iteration scheme,
X(t + dt) = X(t) + N(0, delta**2 * dt; t, t+dt)
If `x0` is an array (or array-like), each value in `x0` is treated as
an initial condition, and the value returned is a numpy array with one
more dimension than `x0`.
Arguments
---------
x0 : float or numpy array (or something that can be converted to a numpy array
using numpy.asarray(x0)).
The initial condition(s) (i.e. position(s)) of the Brownian motion.
n : int
The number of steps to take.
dt : float
The time step.
delta : float
delta determines the "speed" of the Brownian motion. 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.
out : numpy array or None
If `out` is not None, it specifies the array in which to put the
result. If `out` is None, a new numpy array is created and returned.
Returns
-------
A numpy array of floats with shape `x0.shape + (n,)`.
Note that the initial value `x0` is not included in the returned array.
"""
x0 = np.asarray(x0)
# For each element of x0, generate a sample of n numbers from a
# normal distribution.
r = norm.rvs(size=x0.shape + (n,), scale=delta*sqrt(dt))
# If `out` was not given, create an output array.
if out is None:
out = np.empty(r.shape)
# This computes the Brownian motion by forming the cumulative sum of
# the random samples.
np.cumsum(r, axis=-1, out=out)
# Add the initial condition.
out += np.expand_dims(x0, axis=-1)
return out
示例¶
这是一个使用该函数和 matplotlib 的 pylab 模块来绘制多个布朗运动实现的脚本。
在 [5] 中
import numpy
from pylab import plot, show, grid, xlabel, ylabel
# The Wiener process parameter.
delta = 2
# Total time.
T = 10.0
# Number of steps.
N = 500
# Time step size
dt = T/N
# Number of realizations to generate.
m = 20
# Create an empty array to store the realizations.
x = numpy.empty((m,N+1))
# Initial values of x.
x[:, 0] = 50
brownian(x[:,0], N, dt, delta, out=x[:,1:])
t = numpy.linspace(0.0, N*dt, N+1)
for k in range(m):
plot(t, x[k])
xlabel('t', fontsize=16)
ylabel('x', fontsize=16)
grid(True)
show()
在 [6] 中
import numpy
from pylab import plot, show, grid, axis, xlabel, ylabel, title
# The Wiener process parameter.
delta = 0.25
# Total time.
T = 10.0
# Number of steps.
N = 500
# Time step size
dt = T/N
# Initial values of x.
x = numpy.empty((2,N+1))
x[:, 0] = 0.0
brownian(x[:,0], N, dt, delta, out=x[:,1:])
# Plot the 2D trajectory.
plot(x[0],x[1])
# Mark the start and end points.
plot(x[0,0],x[1,0], 'go')
plot(x[0,-1], x[1,-1], 'ro')
# More plot decorations.
title('2D Brownian Motion')
xlabel('x', fontsize=16)
ylabel('y', fontsize=16)
axis('equal')
grid(True)
show()
章节作者：WarrenWeckesser，Warren Weckesser