在 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 Cookbook 编写的代码：

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()