免费网站建设案例,四川省建设领域信用系统网站,wordpress上帝模式,东营会计信息网官网首页目录
Python实例题
题目
处理复杂概率分布的方法
自定义概率分布#xff1a;使用scipy.stats创建自定义分布
使用混合分布#xff1a;将多个简单分布组合成复杂分布
接受 - 拒绝采样#xff1a;从复杂分布中生成样本
MCMC 方法#xff1a;使用马尔可夫链蒙特卡洛采样…目录
Python实例题
题目
处理复杂概率分布的方法
自定义概率分布使用scipy.stats创建自定义分布
使用混合分布将多个简单分布组合成复杂分布
接受 - 拒绝采样从复杂分布中生成样本
MCMC 方法使用马尔可夫链蒙特卡洛采样
Python 计算常微分方程
使用 scipy.integrate.odeint最常用的 ODE 求解器
使用 scipy.integrate.solve_ivp更现代的 ODE 求解器
求解系统 ODE多个相互关联的微分方程
高阶 ODE 转换为一阶系统处理二阶及以上的 ODE
边界值问题 (BVP)使用 scipy.integrate.solve_bvp
处理复杂概率分布与 ODE 结合的随机微分方程 Python实例题
题目
Python计算常微分方程
处理复杂概率分布的方法
在计算随机过程时处理复杂概率分布是常见的挑战。以下是几种有效的方法 自定义概率分布使用scipy.stats创建自定义分布
from scipy import stats
import numpy as np
import matplotlib.pyplot as plt# 定义自定义概率密度函数
def custom_pdf(x):return np.exp(-x**2/2) * (1 0.5 * np.sin(2*x)) / np.sqrt(2*np.pi)# 创建自定义分布
class CustomDistribution(stats.rv_continuous):def _pdf(self, x):return custom_pdf(x)# 实例化分布
custom_dist CustomDistribution(namecustom)# 生成随机样本
samples custom_dist.rvs(size1000)# 可视化
plt.figure(figsize(10, 6))
plt.hist(samples, bins50, densityTrue, alpha0.6, label样本分布)
x np.linspace(-4, 4, 1000)
plt.plot(x, custom_pdf(x), r-, label理论PDF)
plt.legend()
plt.title(自定义概率分布)
plt.show()使用混合分布将多个简单分布组合成复杂分布
# 混合高斯分布
def mixture_gaussian(x, weights, means, stds):result 0for w, mu, sigma in zip(weights, means, stds):result w * stats.norm.pdf(x, mu, sigma)return result# 参数
weights [0.3, 0.5, 0.2]
means [-2, 0, 3]
stds [0.8, 1.2, 0.7]# 生成随机样本
n_samples 1000
component np.random.choice(len(weights), sizen_samples, pweights)
samples np.zeros(n_samples)for i in range(n_samples):samples[i] np.random.normal(means[component[i]], stds[component[i]])# 可视化
plt.figure(figsize(10, 6))
plt.hist(samples, bins50, densityTrue, alpha0.6, label样本分布)
x np.linspace(-6, 6, 1000)
plt.plot(x, mixture_gaussian(x, weights, means, stds), r-, label理论PDF)
plt.legend()
plt.title(混合高斯分布)
plt.show()接受 - 拒绝采样从复杂分布中生成样本
def complex_pdf(x):# 复杂分布的PDFreturn 0.3 * np.exp(-0.2 * x**2) 0.7 * np.exp(-0.2 * (x-10)**2)def sampling_complex(n_samples, proposal_diststats.norm, c1.5):samples []while len(samples) n_samples:# 从提议分布中采样x proposal_dist.rvs()# 计算接受概率accept_prob complex_pdf(x) / (c * proposal_dist.pdf(x))# 决定是否接受if np.random.random() accept_prob:samples.append(x)return np.array(samples)# 生成样本
samples sampling_complex(1000)# 可视化
plt.figure(figsize(10, 6))
plt.hist(samples, bins50, densityTrue, alpha0.6, label样本分布)
x np.linspace(-5, 15, 1000)
plt.plot(x, complex_pdf(x), r-, label理论PDF)
plt.legend()
plt.title(接受-拒绝采样)
plt.show()MCMC 方法使用马尔可夫链蒙特卡洛采样
def metropolis_hastings(target_pdf, n_samples, proposal_std1.0, x00):samples [x0]current_x x0for i in range(n_samples - 1):# 提议新样本proposed_x current_x np.random.normal(0, proposal_std)# 计算接受概率ratio target_pdf(proposed_x) / target_pdf(current_x)accept_prob min(1, ratio)# 决定是否接受if np.random.random() accept_prob:current_x proposed_xsamples.append(current_x)return np.array(samples)# 使用Metropolis-Hastings从复杂分布采样
samples metropolis_hastings(complex_pdf, 5000)
samples samples[1000:] # 丢弃burn-in期# 可视化
plt.figure(figsize(10, 6))
plt.hist(samples, bins50, densityTrue, alpha0.6, label样本分布)
x np.linspace(-5, 15, 1000)
plt.plot(x, complex_pdf(x), r-, label理论PDF)
plt.legend()
plt.title(MCMC采样)
plt.show()Python 计算常微分方程
Python 提供了多种方法来求解常微分方程 (ODEs)。下面介绍主要的方法和库 使用 scipy.integrate.odeint最常用的 ODE 求解器
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt# 定义微分方程 dy/dt f(y, t)
def model(y, t, k):dydt -k * yreturn dydt# 初始条件
y0 5# 时间点
t np.linspace(0, 20, 100)# 求解ODE
k 0.1
y odeint(model, y0, t, args(k,))# 可视化结果
plt.figure(figsize(10, 6))
plt.plot(t, y, b-, linewidth2)
plt.xlabel(时间)
plt.ylabel(y(t))
plt.title(一阶常微分方程 dy/dt -ky 的解)
plt.grid(True)
plt.show()使用 scipy.integrate.solve_ivp更现代的 ODE 求解器
from scipy.integrate import solve_ivp# 定义微分方程 dy/dt f(t, y)
def model(t, y, k):return -k * y# 时间跨度
t_span (0, 20)
t_eval np.linspace(0, 20, 100)# 初始条件
y0 [5]# 求解ODE
k 0.1
sol solve_ivp(lambda t, y: model(t, y, k), t_span, y0, t_evalt_eval)# 可视化结果
plt.figure(figsize(10, 6))
plt.plot(sol.t, sol.y[0], b-, linewidth2)
plt.xlabel(时间)
plt.ylabel(y(t))
plt.title(使用solve_ivp求解一阶常微分方程)
plt.grid(True)
plt.show()求解系统 ODE多个相互关联的微分方程
# 定义Lotka-Volterra捕食者-猎物模型
def lotka_volterra(t, z, a, b, c, d):x, y zdxdt a * x - b * x * ydydt -c * y d * x * yreturn [dxdt, dydt]# 参数
a 1.0
b 0.1
c 1.0
d 0.02# 初始条件
z0 [10, 5] # 初始猎物数量和捕食者数量# 时间跨度
t_span (0, 15)
t_eval np.linspace(0, 15, 1000)# 求解ODE
sol solve_ivp(lambda t, z: lotka_volterra(t, z, a, b, c, d), t_span, z0, t_evalt_eval, methodRK45)# 可视化结果
plt.figure(figsize(12, 6))
plt.plot(sol.t, sol.y[0], b-, label猎物数量)
plt.plot(sol.t, sol.y[1], r-, label捕食者数量)
plt.xlabel(时间)
plt.ylabel(种群数量)
plt.title(Lotka-Volterra捕食者-猎物模型)
plt.legend()
plt.grid(True)
plt.show()高阶 ODE 转换为一阶系统处理二阶及以上的 ODE
# 阻尼振荡器方程: y 2*zeta*omega_n*y omega_n^2*y 0
def oscillator(t, z, omega_n, zeta):y, v z # y是位置v是速度dydt vdvdt -omega_n**2 * y - 2 * zeta * omega_n * vreturn [dydt, dvdt]# 参数
omega_n 2.0 # 自然频率
zeta 0.1 # 阻尼比# 初始条件
z0 [1.0, 0.0] # 初始位置和速度# 时间跨度
t_span (0, 10)
t_eval np.linspace(0, 10, 1000)# 求解ODE
sol solve_ivp(lambda t, z: oscillator(t, z, omega_n, zeta), t_span, z0, t_evalt_eval)# 可视化结果
plt.figure(figsize(12, 6))
plt.plot(sol.t, sol.y[0], b-, label位置)
plt.plot(sol.t, sol.y[1], r-, label速度)
plt.xlabel(时间)
plt.ylabel(状态)
plt.title(阻尼振荡器)
plt.legend()
plt.grid(True)
plt.show()边界值问题 (BVP)使用 scipy.integrate.solve_bvp
from scipy.integrate import solve_bvp# 定义边界值问题 y y 0
def bvp_ode(x, y):return np.vstack((y[1], -y[0]))# 定义边界条件 y(0) 0, y(pi/2) 1
def bvp_bc(ya, yb):return np.array([ya[0], yb[0] - 1])# 初始网格和猜测
x np.linspace(0, np.pi/2, 5)
y np.zeros((2, x.size))
y[0] x # 初始猜测# 求解BVP
sol solve_bvp(bvp_ode, bvp_bc, x, y)# 可视化结果
x_plot np.linspace(0, np.pi/2, 100)
y_plot sol.sol(x_plot)[0]plt.figure(figsize(10, 6))
plt.plot(x_plot, y_plot, b-, linewidth2)
plt.xlabel(x)
plt.ylabel(y(x))
plt.title(边界值问题 y\\ y 0 的解)
plt.grid(True)
plt.show()处理复杂概率分布与 ODE 结合的随机微分方程
当需要求解随机微分方程 (SDE) 时上述技术可以结合使用
# 几何布朗运动的随机微分方程 dS mu*S*dt sigma*S*dW
def geometric_brownian_motion(S0, mu, sigma, T, N, M):模拟几何布朗运动参数:S0: 初始价格mu: 漂移率sigma: 波动率T: 总时间N: 时间步数M: 模拟路径数dt T / Nt np.linspace(0, T, N1)# 生成随机增量dW np.random.normal(0, np.sqrt(dt), (M, N))# 计算价格路径S np.zeros((M, N1))S[:, 0] S0for i in range(N):# Euler-Maruyama方法S[:, i1] S[:, i] * (1 mu * dt sigma * dW[:, i])return t, S# 参数
S0 100 # 初始价格
mu 0.05 # 漂移率
sigma 0.2 # 波动率
T 1.0 # 总时间
N 252 # 交易日数
M 10 # 模拟路径数# 模拟
t, S geometric_brownian_motion(S0, mu, sigma, T, N, M)# 可视化
plt.figure(figsize(12, 6))
for i in range(M):plt.plot(t, S[i], labelf路径 {i1})plt.xlabel(时间)
plt.ylabel(价格)
plt.title(几何布朗运动模拟)
plt.legend()
plt.grid(True)
plt.show()
