绘制 gBM 问题
Plotting gBM issues
Geometric Brownian motion (gBM)是一个随机过程，可以认为是标准布朗运动的扩展。
我正在尝试编写一个函数来模拟 gBM 的不同路径（ntraj 路径），然后在列表 tcheck 中指定的特定点绘制直方图。一旦它绘制了这些图，该函数就意味着每次在图上叠加一个对数正态分布。
输出应该是这样的
除了 gBM 而不是标准的布朗运动过程。到目前为止，我有一个函数可以生成 gBM 的多个路径，
def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
'''
DOCSTRING:
1D geomwtric brownian motion
INPUTS:
ntraj = "number of trajectories"
n = "length of a trajectory"
T = "last time point, i.e final tradjectory t = {0,1,...,T}"
sigma= volatility
mu= percentage drift
'''
np.random.seed(52323)
S_0 = 0
# Discretize, dt = time step = $t_{j+1}- t_{j}$
dt = T/(n)
sqrtdt = np.sqrt(dt)
# Container for different colors for each trajectory
colors = plt.cm.jet(np.linspace(0,1,nTraj))
# Container for trajectories
xtraj=np.zeros(n+1,float)
ztraj=np.zeros(n+1,float)
trange=np.linspace(start = 0,stop = T ,num = n+1)
# Simulation
# Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1)
for j in range(nTraj):
# Loop over time
for i in range(n):
xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu
# Loop again over time in order to make geometric drift
ztraj = S_0 * np.exp(xtraj) # ztraj[z+1]= ztraj[0]+ np.exp(xtraj[z])
plt.plot(trange , xtraj,'b-',alpha=0.2, color=colors[j], lw=3.0,label="$\sigma$={}, $\mu$={:.5f}".format(sigma,mu))
plt.title("1D Geometric Brownian Motion:\n nTraj={}, T={},\n $\Delta t$={:.3f}, $\sigma$={}, $\mu$={:.3f}".format(nTraj,T,dt,sigma,mu))
plt.xlabel(r'$t$')
plt.ylabel(r'$Z_t$');
oneDGeometricBM(nTraj=5,n=10**3,T=10.0,sigma=0.8,mu=1.1)
我已经看到很多关于如何绘制 gBM 的多条路径的问题的答案，但我对如何查看特定时间的直方图然后查看分布感兴趣。以下是我到目前为止的功能。它不起作用，但我无法弄清楚我做错了什么。我还添加了我得到的输出。
import yfinance as yf
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.stats import norm, lognorm
ntraj = 10000
S_0 =0
sigma=1
mu=1
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 1.0
'''
ntraj = 10**4
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 5.0 # limits
'''
n=int(tfinal/dt)
sqrtdt = np.sqrt(dt)
x=np.zeros(shape=[ntraj,n+1], dtype=float)
z=np.zeros(shape=[ntraj,n+1], dtype=float)
zrange=np.arange(start=-xv, stop=xv, step=dt)
# Calculate the number of the bins
binval = math.ceil(np.sqrt(ntraj))
# Nested for loop to create Drifted BM
for i in range(n):
for j in range(ntraj):
x[j,i+1]=x[j,i]+ sqrtdt*np.random.randn()
#Nested loop to create gBM
for j0 in range(ntraj):
for i0 in range(n+1):
z[j0,i0] = 0 + np.exp(x[j0,i0])
# Loop to plot the distribution of gBM tradjectories at different times
for i1 in range(n):
# Compute histogram at every tsample , sample at time t
t=(i1+1)*dt
if t in tcheck:
# Plot histogram on sample
plt.hist(z[:,i1],bins=30,density=False,alpha=0.6,label=['t ={}'.format(t)] )
# Superimpose each samples mean
xbar = np.average(z[:,i1])
plt.axvline(xbar, color='RED', linestyle='dashed', linewidth=2)
# Plot theoretic distribution { N(0, sqrt[t] ) }
#plt.plot(xrange,norm.pdf(xrange,0.0,np.sqrt(t)),'k--')
总结一下我的问题。我正在尝试模拟 gBM 的多个轨迹，将我的结果存储在一个数组中，然后遍历该数组并使用 matplotlib 在特定点上绘制直方图，然后最后在我的直方图上叠加对数正态分布。
编辑 1：
如果可能的话，我需要在 GBM 和 Cauchy 上叠加对数正态分布。我的问题是，当我编辑@Paul Harris 的更正时，我得到了，
def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
'''
DOCSTRING:
INPUTS:
ntraj = "number of trajectories"
n = "length of a trajectory"
T = "last time point, i.e final tradjectory t = {0,1,...,T}"
sigma= volatility
mu= percentage drift
'''
np.random.seed(52323)
S0 = 10
# Discretize, dt = time step = $t_{j+1}- t_{j}$
dt = T/(n)
sqrtdt = np.sqrt(dt)
# Container for different colors for each trajectory
colors = plt.cm.jet(np.linspace(0,1,nTraj))
# Container for trajectories
xtraj=np.zeros(n+1,float)
ztraj=np.zeros(n+1,float)
trange=np.linspace(start = 0,stop = T ,num = n+1)
out = []
# Simulation
# Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1)
for j in range(nTraj):
# Loop over time
for i in range(n):
xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu
# Loop again over time in order to make geometric drift
ztraj = S0 * np.exp(xtraj)
# Return gBM
return ztraj
# Plotting
fig, ax = plt.subplots(ncols=2, figsize=plt.figaspect(1./2))
colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]
sigma=0.8
mu=1.1
for c, T in zip(colors, T):
ztraj = oneDGeometricBM(nTraj=5,n=10**4,T=T,sigma=0.8,mu=1.1)
# Plot Emperical Values
xrange = range(0,80,1)
ax[0].hist(ztraj, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(0, 95))
# Plot the theoretical values
theoretic_mean = math.exp(mu * T + 0.5 * sigma**2 * T)
theoretic_var = math.exp(2* mu * T + sigma**2 * T)* (math.exp(sigma**2 * T) - 1)
ax[0].plot(xrange,lognorm.pdf(xrange, theoretic_mean , theoretic_var ),'k--')
# Plot the differences between consecutive elements of gBM (an array)
diff = np.ediff1d(ztraj)
ax[1].hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-5, 5))
ax[0].set_xlabel('z')
ax[0].set_ylabel('$p(z,T)$')
ax[0].set_title('Histogram of ztraj positions')
ax[1].set_xlabel('dz')
ax[1].set_ylabel('$p(dz,T)$')
ax[1].set_title('Histogram of d(ztraj) positions\nbetween time steps')
ax[0].legend()
fig.tight_layout()
所以总而言之，我需要叠加每个时间点的分布，即 gBM 的理论分布，即对数正态分布。
所以我已经看过你的问题了。我已经编辑了你的函数以停止绘图和 return xtraj 我认为这是你的布朗运动：
def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
'''
DOCSTRING:
1D geomwtric brownian motion
INPUTS:
ntraj = "number of trajectories"
n = "length of a trajectory"
T = "last time point, i.e final tradjectory t = {0,1,...,T}"
sigma= volatility
mu= percentage drift
'''
np.random.seed(52323)
S_0 = 10
# Discretize, dt = time step = $t_{j+1}- t_{j}$
dt = T/(n)
sqrtdt = np.sqrt(dt)
# Container for different colors for each trajectory
colors = plt.cm.jet(np.linspace(0,1,nTraj))
# Container for trajectories
xtraj=np.zeros(n+1,float)
ztraj=np.zeros(n+1,float)
trange=np.linspace(start = 0,stop = T ,num = n+1)
out = []
# Simulation
# Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1)
for j in range(nTraj):
# Loop over time
for i in range(n):
xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu
# Loop again over time in order to make geometric drift
ztraj = S_0 * np.exp(xtraj) # ztraj[z+1]= ztraj[0]+ np.exp(xtraj[z])
return ztraj
每个时间步长的位移就是数组 xtraj 中的差异：dx = np.ediff1d(oneDGeometricBM(...))，因此我们计算这些值的直方图：
fig, ax = plt.subplots()
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=10.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, label='T=10', density=True)
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=1.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, color='k', label='T=1', density=True)
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=5.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, color='r', label='T=5', density=True)
ax.set_xlabel('x')
ax.set_ylabel('$p(x,T)$')
ax.legend()
我使用了 3 个不同的 T 值，如示例中所示。标准化直方图，使 y 轴现在代表概率 p(x,T)，即。所有 p*x = 1 的总和，我们使用 density=True 参数。
编辑
我已将 oneDGeometricBM 函数编辑为 return ztraj = S0*np.exp(xtraj)。您的初始 S0 值为 0，因此我将其设置为非零值。
您可以将 ztraj 差异绘制为：
fig, ax = plt.subplots()
colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]
for c, T in zip(colors, T):
ztraj = oneDGeometricBM(nTraj=5,n=10**3,T=T,sigma=0.8,mu=1.1)
diff = np.ediff1d(ztraj)
ax.hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-10, 10))
ax.set_xlabel('x')
ax.set_ylabel('$p(x,T)$')
ax.legend()
编辑2
通过仔细查看您生成的直方图，我认为您的建模是正确的，只是应该调整图的 xrange，因为 ztraj 对于大 T 变大，您可以限制直方图使用 range 参数。所以我为三个单独的 T 绘制了 ztraj 和 d(ztraj)。 ztraj 确实似乎近似服从对数正态分布，而 ztraj 的差异似乎近似服从洛伦兹分布（必须检查那个理论，也许是高斯分布）。重现代码：
fig, ax = plt.subplots(ncols=2, figsize=plt.figaspect(1./2))
colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]
for c, T in zip(colors, T):
ztraj = oneDGeometricBM(nTraj=5,n=10**4,T=T,sigma=0.8,mu=1.1)
ax[0].hist(ztraj, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(0, 95))
diff = np.ediff1d(ztraj)
ax[1].hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-5, 5))
ax[0].set_xlabel('z')
ax[0].set_ylabel('$p(z,T)$')
ax[0].set_title('Histogram of ztraj positions')
ax[1].set_xlabel('dz')
ax[1].set_ylabel('$p(dz,T)$')
ax[1].set_title('Histogram of d(ztraj) positions\nbetween time steps')
ax[0].legend()
fig.tight_layout()
这是您的数据和绘图，但限制了直方图 range=(0, 10)：
编辑3
我已经包含了适合对数正态分布的代码，并在您的原始图中显示了它们。我们将 lognormal function 定义为：
from scipy.optimize import curve_fit
def lognorm(x, x0, A, sigma):
return A * np.exp(-(np.log(x)-x0)**2 / (2*sigma**2))
然后在最终循环中使用直方图中的值和分箱进行拟合：
# Loop to plot the distribution of gBM tradjectories at different times
for i1 in range(n):
# Compute histogram at every tsample , sample at time t
t=(i1+1)*dt
if t in tcheck:
# Plot histogram on sample
v, b, patches = plt.hist(z[:,i1],bins=200,density=False,alpha=0.6,label=['t ={}'.format(t)], range=(0, 10) )
# second term is bin centre locations rather than bin edges
popt, pcov = curve_fit(lognorm, b[:-1] + np.ediff1d(b), v, p0=(0.1, 300, 0.3))
# make colors match their original data but no transparency
plt.plot(b, lognorm(b, *popt), color=patches[0].get_facecolor()[:3])
print(f'tcheck: {t} with parameters: {popt}')
输出：
tcheck: 0.5 with parameters: [ -0.42334757 358.38545736 0.6748076 ]
tcheck: 1.0 with parameters: [ -0.90719967 321.03944864 0.96137893]
tcheck: 4.0 with parameters: [ -3.66426932 721.41708932 1.86376987]
编辑4
生成上述输出的全部代码为：
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.stats import norm, lognorm
from scipy.optimize import curve_fit
def lognorm(x, x0, A, sigma):
return A * np.exp(-(np.log(x)-x0)**2 / (2*sigma**2))
ntraj = 10000
S_0 =0
sigma=1
mu=1
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 1.0
'''
ntraj = 10**4
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 5.0 # limits
'''
n=int(tfinal/dt)
sqrtdt = np.sqrt(dt)
x=np.zeros(shape=[ntraj,n+1], dtype=float)
z=np.zeros(shape=[ntraj,n+1], dtype=float)
zrange=np.arange(start=-xv, stop=xv, step=dt)
# Calculate the number of the bins
binval = math.ceil(np.sqrt(ntraj))
# Nested for loop to create Drifted BM
for i in range(n):
for j in range(ntraj):
x[j,i+1]=x[j,i]+ sqrtdt*np.random.randn()
#Nested loop to create gBM
for j0 in range(ntraj):
for i0 in range(n+1):
z[j0,i0] = 0 + np.exp(x[j0,i0])
# Loop to plot the distribution of gBM tradjectories at different times
for i1 in range(n):
# Compute histogram at every tsample , sample at time t
t=(i1+1)*dt
if t in tcheck:
# Plot histogram on sample
v, b, patches = plt.hist(z[:,i1],bins=200,density=True,alpha=0.6,label=['t ={}'.format(t)], range=(0, 10))
popt, pcov = curve_fit(lognorm, b[:-1] + np.ediff1d(b), v, p0=(0.1, 300, 0.3))
# make colors match their original data but no transparency
plt.plot(b, lognorm(b, *popt), color=patches[0].get_facecolor()[:3])
print(f'tcheck: {t} for parameters: {popt}')
Geometric Brownian motion (gBM)是一个随机过程，可以认为是标准布朗运动的扩展。
我正在尝试编写一个函数来模拟 gBM 的不同路径（ntraj 路径），然后在列表 tcheck 中指定的特定点绘制直方图。一旦它绘制了这些图，该函数就意味着每次在图上叠加一个对数正态分布。
输出应该是这样的
除了 gBM 而不是标准的布朗运动过程。到目前为止，我有一个函数可以生成 gBM 的多个路径，
def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
'''
DOCSTRING:
1D geomwtric brownian motion
INPUTS:
ntraj = "number of trajectories"
n = "length of a trajectory"
T = "last time point, i.e final tradjectory t = {0,1,...,T}"
sigma= volatility
mu= percentage drift
'''
np.random.seed(52323)
S_0 = 0
# Discretize, dt = time step = $t_{j+1}- t_{j}$
dt = T/(n)
sqrtdt = np.sqrt(dt)
# Container for different colors for each trajectory
colors = plt.cm.jet(np.linspace(0,1,nTraj))
# Container for trajectories
xtraj=np.zeros(n+1,float)
ztraj=np.zeros(n+1,float)
trange=np.linspace(start = 0,stop = T ,num = n+1)
# Simulation
# Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1)
for j in range(nTraj):
# Loop over time
for i in range(n):
xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu
# Loop again over time in order to make geometric drift
ztraj = S_0 * np.exp(xtraj) # ztraj[z+1]= ztraj[0]+ np.exp(xtraj[z])
plt.plot(trange , xtraj,'b-',alpha=0.2, color=colors[j], lw=3.0,label="$\sigma$={}, $\mu$={:.5f}".format(sigma,mu))
plt.title("1D Geometric Brownian Motion:\n nTraj={}, T={},\n $\Delta t$={:.3f}, $\sigma$={}, $\mu$={:.3f}".format(nTraj,T,dt,sigma,mu))
plt.xlabel(r'$t$')
plt.ylabel(r'$Z_t$');
oneDGeometricBM(nTraj=5,n=10**3,T=10.0,sigma=0.8,mu=1.1)
我已经看到很多关于如何绘制 gBM 的多条路径的问题的答案，但我对如何查看特定时间的直方图然后查看分布感兴趣。以下是我到目前为止的功能。它不起作用，但我无法弄清楚我做错了什么。我还添加了我得到的输出。
import yfinance as yf
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.stats import norm, lognorm
ntraj = 10000
S_0 =0
sigma=1
mu=1
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 1.0
'''
ntraj = 10**4
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 5.0 # limits
'''
n=int(tfinal/dt)
sqrtdt = np.sqrt(dt)
x=np.zeros(shape=[ntraj,n+1], dtype=float)
z=np.zeros(shape=[ntraj,n+1], dtype=float)
zrange=np.arange(start=-xv, stop=xv, step=dt)
# Calculate the number of the bins
binval = math.ceil(np.sqrt(ntraj))
# Nested for loop to create Drifted BM
for i in range(n):
for j in range(ntraj):
x[j,i+1]=x[j,i]+ sqrtdt*np.random.randn()
#Nested loop to create gBM
for j0 in range(ntraj):
for i0 in range(n+1):
z[j0,i0] = 0 + np.exp(x[j0,i0])
# Loop to plot the distribution of gBM tradjectories at different times
for i1 in range(n):
# Compute histogram at every tsample , sample at time t
t=(i1+1)*dt
if t in tcheck:
# Plot histogram on sample
plt.hist(z[:,i1],bins=30,density=False,alpha=0.6,label=['t ={}'.format(t)] )
# Superimpose each samples mean
xbar = np.average(z[:,i1])
plt.axvline(xbar, color='RED', linestyle='dashed', linewidth=2)
# Plot theoretic distribution { N(0, sqrt[t] ) }
#plt.plot(xrange,norm.pdf(xrange,0.0,np.sqrt(t)),'k--')
总结一下我的问题。我正在尝试模拟 gBM 的多个轨迹，将我的结果存储在一个数组中，然后遍历该数组并使用 matplotlib 在特定点上绘制直方图，然后最后在我的直方图上叠加对数正态分布。
编辑 1：
如果可能的话，我需要在 GBM 和 Cauchy 上叠加对数正态分布。我的问题是，当我编辑@Paul Harris 的更正时，我得到了，
def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
'''
DOCSTRING:
INPUTS:
ntraj = "number of trajectories"
n = "length of a trajectory"
T = "last time point, i.e final tradjectory t = {0,1,...,T}"
sigma= volatility
mu= percentage drift
'''
np.random.seed(52323)
S0 = 10
# Discretize, dt = time step = $t_{j+1}- t_{j}$
dt = T/(n)
sqrtdt = np.sqrt(dt)
# Container for different colors for each trajectory
colors = plt.cm.jet(np.linspace(0,1,nTraj))
# Container for trajectories
xtraj=np.zeros(n+1,float)
ztraj=np.zeros(n+1,float)
trange=np.linspace(start = 0,stop = T ,num = n+1)
out = []
# Simulation
# Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1)
for j in range(nTraj):
# Loop over time
for i in range(n):
xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu
# Loop again over time in order to make geometric drift
ztraj = S0 * np.exp(xtraj)
# Return gBM
return ztraj
# Plotting
fig, ax = plt.subplots(ncols=2, figsize=plt.figaspect(1./2))
colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]
sigma=0.8
mu=1.1
for c, T in zip(colors, T):
ztraj = oneDGeometricBM(nTraj=5,n=10**4,T=T,sigma=0.8,mu=1.1)
# Plot Emperical Values
xrange = range(0,80,1)
ax[0].hist(ztraj, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(0, 95))
# Plot the theoretical values
theoretic_mean = math.exp(mu * T + 0.5 * sigma**2 * T)
theoretic_var = math.exp(2* mu * T + sigma**2 * T)* (math.exp(sigma**2 * T) - 1)
ax[0].plot(xrange,lognorm.pdf(xrange, theoretic_mean , theoretic_var ),'k--')
# Plot the differences between consecutive elements of gBM (an array)
diff = np.ediff1d(ztraj)
ax[1].hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-5, 5))
ax[0].set_xlabel('z')
ax[0].set_ylabel('$p(z,T)$')
ax[0].set_title('Histogram of ztraj positions')
ax[1].set_xlabel('dz')
ax[1].set_ylabel('$p(dz,T)$')
ax[1].set_title('Histogram of d(ztraj) positions\nbetween time steps')
ax[0].legend()
fig.tight_layout()
所以总而言之，我需要叠加每个时间点的分布，即 gBM 的理论分布，即对数正态分布。
所以我已经看过你的问题了。我已经编辑了你的函数以停止绘图和 return xtraj 我认为这是你的布朗运动：
def oneDGeometricBM(nTraj=100,n=100,T=1.0,sigma=1,mu=0):
'''
DOCSTRING:
1D geomwtric brownian motion
INPUTS:
ntraj = "number of trajectories"
n = "length of a trajectory"
T = "last time point, i.e final tradjectory t = {0,1,...,T}"
sigma= volatility
mu= percentage drift
'''
np.random.seed(52323)
S_0 = 10
# Discretize, dt = time step = $t_{j+1}- t_{j}$
dt = T/(n)
sqrtdt = np.sqrt(dt)
# Container for different colors for each trajectory
colors = plt.cm.jet(np.linspace(0,1,nTraj))
# Container for trajectories
xtraj=np.zeros(n+1,float)
ztraj=np.zeros(n+1,float)
trange=np.linspace(start = 0,stop = T ,num = n+1)
out = []
# Simulation
# Random Variable $X_{n}$ is distributed np.sqrt(dt)* N(mu=0,sigma=1)
for j in range(nTraj):
# Loop over time
for i in range(n):
xtraj[i+1]=xtraj[i]+ sqrtdt * np.random.randn() + dt*mu
# Loop again over time in order to make geometric drift
ztraj = S_0 * np.exp(xtraj) # ztraj[z+1]= ztraj[0]+ np.exp(xtraj[z])
return ztraj
每个时间步长的位移就是数组 xtraj 中的差异：dx = np.ediff1d(oneDGeometricBM(...))，因此我们计算这些值的直方图：
fig, ax = plt.subplots()
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=10.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, label='T=10', density=True)
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=1.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, color='k', label='T=1', density=True)
ax.hist(np.ediff1d(oneDGeometricBM(nTraj=5,n=10**3,T=5.0,sigma=0.8,mu=1.1)), bins=50, alpha=0.5, color='r', label='T=5', density=True)
ax.set_xlabel('x')
ax.set_ylabel('$p(x,T)$')
ax.legend()
我使用了 3 个不同的 T 值，如示例中所示。标准化直方图，使 y 轴现在代表概率 p(x,T)，即。所有 p*x = 1 的总和，我们使用 density=True 参数。
编辑
我已将 oneDGeometricBM 函数编辑为 return ztraj = S0*np.exp(xtraj)。您的初始 S0 值为 0，因此我将其设置为非零值。
您可以将 ztraj 差异绘制为：
fig, ax = plt.subplots()
colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]
for c, T in zip(colors, T):
ztraj = oneDGeometricBM(nTraj=5,n=10**3,T=T,sigma=0.8,mu=1.1)
diff = np.ediff1d(ztraj)
ax.hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-10, 10))
ax.set_xlabel('x')
ax.set_ylabel('$p(x,T)$')
ax.legend()
编辑2
通过仔细查看您生成的直方图，我认为您的建模是正确的，只是应该调整图的 xrange，因为 ztraj 对于大 T 变大，您可以限制直方图使用 range 参数。所以我为三个单独的 T 绘制了 ztraj 和 d(ztraj)。 ztraj 确实似乎近似服从对数正态分布，而 ztraj 的差异似乎近似服从洛伦兹分布（必须检查那个理论，也许是高斯分布）。重现代码：
fig, ax = plt.subplots(ncols=2, figsize=plt.figaspect(1./2))
colors = ['k', 'r', 'b']
T = [1.0, 2.0, 5.0]
for c, T in zip(colors, T):
ztraj = oneDGeometricBM(nTraj=5,n=10**4,T=T,sigma=0.8,mu=1.1)
ax[0].hist(ztraj, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(0, 95))
diff = np.ediff1d(ztraj)
ax[1].hist(diff, bins=100, alpha=0.5, label=f'T={T}', density=True, color=c, range=(-5, 5))
ax[0].set_xlabel('z')
ax[0].set_ylabel('$p(z,T)$')
ax[0].set_title('Histogram of ztraj positions')
ax[1].set_xlabel('dz')
ax[1].set_ylabel('$p(dz,T)$')
ax[1].set_title('Histogram of d(ztraj) positions\nbetween time steps')
ax[0].legend()
fig.tight_layout()
这是您的数据和绘图，但限制了直方图 range=(0, 10)：
编辑3
我已经包含了适合对数正态分布的代码，并在您的原始图中显示了它们。我们将 lognormal function 定义为：
from scipy.optimize import curve_fit
def lognorm(x, x0, A, sigma):
return A * np.exp(-(np.log(x)-x0)**2 / (2*sigma**2))
然后在最终循环中使用直方图中的值和分箱进行拟合：
# Loop to plot the distribution of gBM tradjectories at different times
for i1 in range(n):
# Compute histogram at every tsample , sample at time t
t=(i1+1)*dt
if t in tcheck:
# Plot histogram on sample
v, b, patches = plt.hist(z[:,i1],bins=200,density=False,alpha=0.6,label=['t ={}'.format(t)], range=(0, 10) )
# second term is bin centre locations rather than bin edges
popt, pcov = curve_fit(lognorm, b[:-1] + np.ediff1d(b), v, p0=(0.1, 300, 0.3))
# make colors match their original data but no transparency
plt.plot(b, lognorm(b, *popt), color=patches[0].get_facecolor()[:3])
print(f'tcheck: {t} with parameters: {popt}')
输出：
tcheck: 0.5 with parameters: [ -0.42334757 358.38545736 0.6748076 ]
tcheck: 1.0 with parameters: [ -0.90719967 321.03944864 0.96137893]
tcheck: 4.0 with parameters: [ -3.66426932 721.41708932 1.86376987]
编辑4
生成上述输出的全部代码为：
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import math
from scipy.stats import norm, lognorm
from scipy.optimize import curve_fit
def lognorm(x, x0, A, sigma):
return A * np.exp(-(np.log(x)-x0)**2 / (2*sigma**2))
ntraj = 10000
S_0 =0
sigma=1
mu=1
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 1.0
'''
ntraj = 10**4
tfinal = 4.0
tcheck = [0.5, 1.0, 4.0]
dt = 0.01
xv = 5.0 # limits
'''
n=int(tfinal/dt)
sqrtdt = np.sqrt(dt)
x=np.zeros(shape=[ntraj,n+1], dtype=float)
z=np.zeros(shape=[ntraj,n+1], dtype=float)
zrange=np.arange(start=-xv, stop=xv, step=dt)
# Calculate the number of the bins
binval = math.ceil(np.sqrt(ntraj))
# Nested for loop to create Drifted BM
for i in range(n):
for j in range(ntraj):
x[j,i+1]=x[j,i]+ sqrtdt*np.random.randn()
#Nested loop to create gBM
for j0 in range(ntraj):
for i0 in range(n+1):
z[j0,i0] = 0 + np.exp(x[j0,i0])
# Loop to plot the distribution of gBM tradjectories at different times
for i1 in range(n):
# Compute histogram at every tsample , sample at time t
t=(i1+1)*dt
if t in tcheck:
# Plot histogram on sample
v, b, patches = plt.hist(z[:,i1],bins=200,density=True,alpha=0.6,label=['t ={}'.format(t)], range=(0, 10))
popt, pcov = curve_fit(lognorm, b[:-1] + np.ediff1d(b), v, p0=(0.1, 300, 0.3))
# make colors match their original data but no transparency
plt.plot(b, lognorm(b, *popt), color=patches[0].get_facecolor()[:3])
print(f'tcheck: {t} for parameters: {popt}')