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

Question

我试图画出Ornstein-Uhlenbeck过程的时间演化图，这是一个随机过程，然后找到每个时间步的概率分布。我能够为流程的1000实现绘制图表。每个实现都有一个1000时间步长，时间步长的宽度为.001。我使用了一个1000 x 1000数组来存储数据。每行保存每个实现的值。和按列的i-th列对应于i-th时间步上1000实现的值

现在我想把每个时间步的bin结果放在一起，然后绘制每个时间步对应的概率分布。我对做这件事感到很困惑

我从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个bin，所以您的代码可能应该读取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()