Построение множественных реализаций одного процесса c в Python

Question

Я пытаюсь построить график временной эволюции для процесса Орнштейна-Уленбека, который представляет собой процесс c случайности, а затем найти распределение вероятностей на каждом временном шаге. Я могу построить график для 1000 реализации процесса. Каждая реализация имеет временной шаг 1000 с шириной временного шага .001. Я использовал массив 1000 x 1000 для хранения данных. Каждая строка содержит значение каждой реализации. А по столбцам i-th столбцы соответствуют значению на временном шаге i-th для 1000 реализаций.

Теперь мне нужно bin результатов на каждом временном шаге вместе, а затем построить распределение вероятностей, соответствующее каждому временному шагу. Я совершенно запутался в этом ( Я пробовал модифицировать код из I Python Cookbook, где они не хранят каждую реализацию в памяти ).

Код, который я сделал из I Python Поваренная книга:

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)

График для 1000 реализаций процесса Орнштейна-Уленбека:

Распределение, сгенерированное из приведенного выше кода:

Я действительно потерялся с присвоением значения bin и построение гистограммы с ее помощью. Я хочу знать, верен ли мой код для построения распределений, соответствующих каждому временному шагу, используя bin. Если нет, скажите, пожалуйста, какие изменения мне нужно внести в мой код.

Answer 1

Последний для l oop должен повторяться по 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()