在 Python 中绘制随机过程的多个实现

Question

我正在尝试绘制 Ornstein-Uhlenbeck 过程的时间演化图，这是一个随机过程，然后找到每个时间步长的概率分布。我能够为 1000 过程的实现绘制图表。每个实现都有一个 1000 时间步长，时间步长为 .001。我使用 1000 x 1000 数组来存储数据。每行都包含每个实现的值。并且按列 i-th 列对应于 i-th 时间步的值 1000 实现。

现在我想要 bin 每个时间步的结果在一起，然后绘制每个时间步对应的概率分布。我对这样做感到很困惑 (I tried modifying code from IPython Cookbook, where they don't store each realizations in the memory)。

我根据 IPython 食谱编写的代码：

import numpy as np
import matplotlib.pyplot as plt

sigma = 1.  # Standard deviation.
mu = 10.  # Mean.
tau = .05  # Time constant.

dt = .001  # Time step.
T = 1.  # Total time.
n = int(T / dt)  # Number of time steps.
ntrails = 1000 # Number of Realizations.
t = np.linspace(0., T, n)  # Vector of times.

sigmabis = sigma * np.sqrt(2. / tau)
sqrtdt = np.sqrt(dt)

x = np.zeros((ntrails,n))  # Vector containing all successive values of our process

for j in range (ntrails):  # Euler Method
    for i in range(n - 1):     
        x[j,i + 1] = x[j,i] + dt * (-(x[j,i] - mu) / tau) + sigmabis * sqrtdt * np.random.randn()

for k in range(ntrails): #plotting 1000 realizations
    plt.plot(t, x[k])

# Time averaging of each time stamp using bin 

# Really lost from this point onwrds.
bins = np.linspace(-2., 15., 100)
fig, ax = plt.subplots(1, 1, figsize=(12, 4))
for i in range(ntrails):
    hist, _ = np.histogram(x[:,[i]], bins=bins)
    ax.plot(hist)

Ornstein-Uhlenbeck 过程的 1000 个实现图：

从上面的代码生成的分布：

我真的迷失了分配 bin 值并使用它绘制直方图。我想知道我的代码使用 bin 绘制对应于每个时间步长的分布是否正确。如果不是请告诉我我需要对我的代码进行哪些修改。

Answer 1

最后一个 for 循环应该遍历 n，而不是 ntrails（这里恰好是相同的值），否则代码和绘图看起来是正确的（除了一些小问题，例如因为这需要 101 次休息才能获得 100 个箱子，所以您的代码可能应该读作 bins = np.linspace(-2., 15., 101)）。

不过你的情节还可以改进一下。一个好的指导原则是尽可能少地使用墨水来传达您要表达的观点。您总是试图绘制所有数据，这最终会掩盖您的绘图。此外，您可以从更多地关注颜色中受益。颜色应该具有意义，或者根本不使用。

以下是我的建议：

import numpy as np
import matplotlib.pyplot as plt

import matplotlib as mpl
mpl.rcParams['axes.spines.top'] = False
mpl.rcParams['axes.spines.right'] = False

sigma = 1.  # Standard deviation.
mu = 10.  # Mean.
tau = .05  # Time constant.

dt = .001  # Time step.
T = 1  # Total time.
n = int(T / dt)  # Number of time steps.
ntrails = 10000 # Number of Realizations.
t = np.linspace(0., T, n)  # Vector of times.

sigmabis = sigma * np.sqrt(2. / tau)
sqrtdt = np.sqrt(dt)

x = np.zeros((ntrails,n))  # Vector containing all successive values of our process

for j in range(ntrails):  # Euler Method
    for i in range(n - 1):
        x[j,i + 1] = x[j,i] + dt * (-(x[j,i] - mu) / tau) + sigmabis * sqrtdt * np.random.randn()

fig, ax = plt.subplots()
for k in range(200): # plotting fewer realizations shows the distribution better in this case
    ax.plot(t, x[k], color='k', alpha=0.02)

# Really lost from this point onwards.
bins = np.linspace(-2., 15., 101) # you need 101 breaks to get 100 bins
fig, ax = plt.subplots(1, 1, figsize=(12, 4))
# plotting a smaller selection of time points spaced out using a log scale prevents
# the asymptotic approach to the mean from dominating the plot
for i in np.logspace(0, np.log10(n)-1, 21):
    hist, _ = np.histogram(x[:,[int(i)]], bins=bins)
    ax.plot(hist, color=plt.cm.plasma(i/20))
plt.show()

在 Python 中绘制随机过程的多个实现

Plotting Multiple Realizations of a Stochastic Process in Python

python

numpy

matplotlib

histogram

scipy