【机器学习】穷理至极，观微知著：微积分的哲思之旅与算法之道

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import dblquad

# -------------------- 高维期望值部分 --------------------

# 定义联合概率密度函数 f_X(x1, x2) = 6x1x2 for 0 <= x1 <=1 and 0 <= x2 <=1
def joint_pdf(x2, x1):
    return 6 * x1 * x2 if 0 <= x1 <=1 and 0 <= x2 <=1 else 0

# 计算 E[X1] = ∫∫ x1 * f_X(x1, x2) dx2 dx1
E_X1, error1 = dblquad(lambda x2, x1: x1 * joint_pdf(x2, x1), 0, 1, lambda x1: 0, lambda x1:1)
print(f"二维随机变量 X1 的期望值 E[X1]: {E_X1:.4f}, 估计误差: {error1:.4e}")

# 计算 E[X2] = ∫∫ x2 * f_X(x1, x2) dx2 dx1
E_X2, error2 = dblquad(lambda x2, x1: x2 * joint_pdf(x2, x1), 0, 1, lambda x1: 0, lambda x1:1)
print(f"二维随机变量 X2 的期望值 E[X2]: {E_X2:.4f}, 估计误差: {error2:.4e}")

# 可视化联合概率密度函数
x1 = np.linspace(0, 1, 100)
x2 = np.linspace(0, 1, 100)
X1, X2 = np.meshgrid(x1, x2)
F = 6 * X1 * X2

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X1, X2, F, cmap='plasma', alpha=0.7)
ax.set_title('联合概率密度函数 $f_X(x_1, x_2) = 6x_1x_2$', fontsize=14)
ax.set_xlabel('$x_1$', fontsize=12)
ax.set_ylabel('$x_2$', fontsize=12)
ax.set_zlabel('$f_X(x_1, x_2)$', fontsize=12)
plt.show()

二维随机变量 X1 的期望值 E[X1]: 0.5000, 估计误差: 1.3878e-17
二维随机变量 X2 的期望值 E[X2]: 0.5000, 估计误差: 1.3878e-17

import numpy as np
from scipy.special import gamma

# -------------------- 蒙特卡罗方法计算单位球体体积 --------------------

# 定义蒙特卡罗积分函数
def monte_carlo_sphere_volume(dimensions, samples=10**6):
    # 在 [-1, 1] 的超立方体中均匀采样点
    points = np.random.uniform(-1, 1, (samples, dimensions))
    # 计算采样点的欧几里得范数平方
    distances_squared = np.sum(points**2, axis=1)
    # 判断点是否落在单位球体内
    inside_sphere = distances_squared <= 1
    # 计算单位球体体积：球体点占比乘以超立方体体积
    volume = (2**dimensions) * np.mean(inside_sphere)
    return volume

# 计算四维单位球体的体积
dimensions = 4
samples = 10**6  # 设置采样数量
V4_monte_carlo = monte_carlo_sphere_volume(dimensions, samples)
print(f"蒙特卡罗方法计算的四维单位球体体积: {V4_monte_carlo:.6f}")

# 理论公式计算
V4_theoretical = (np.pi**(dimensions / 2)) / gamma((dimensions / 2) + 1)
print(f"四维单位球体的理论体积: {V4_theoretical:.6f}")

蒙特卡罗方法计算的四维单位球体体积: 4.937760
四维单位球体的理论体积: 4.934802

import numpy as np
import matplotlib.pyplot as plt

# -------------------- 蒙特卡罗积分部分 --------------------

# 定义维度（4维单位球体）
n = 4
# 定义样本数量（采样点数）
N = 1000000
# 生成 N 个 n 维均匀分布的样本点在 [-1,1]^n
samples = np.random.uniform(-1, 1, (N, n))
# 计算样本点到原点的距离的平方
dist_squared = np.sum(samples**2, axis=1)
# 统计在单位球体内的样本点数量
count_inside = np.sum(dist_squared <= 1)
# 估计单位球体积
volume_estimate = (2**n) * (count_inside / N)
print(f"四维单位球体的蒙特卡罗估计体积: {volume_estimate:.6f}")

# 理论值计算
from scipy.special import gamma
theoretical_volume = (np.pi**(n/2)) / gamma(n/2 + 1)
print(f"四维单位球体的理论体积: {theoretical_volume:.6f}")

# -------------------- 可视化部分（二维投影） --------------------

if n >= 2:
    plt.figure(figsize=(8, 6))
    # 选择前两个维度进行可视化
    x, y = samples[:, 0], samples[:, 1]
    # 绘制位于单位球体内的点（蓝色）
    plt.scatter(x[dist_squared <= 1], y[dist_squared <= 1], color='blue', alpha=0.1, label='Inside Sphere')
    # 绘制位于单位球体外的点（红色）
    plt.scatter(x[dist_squared > 1], y[dist_squared > 1], color='red', alpha=0.1, label='Outside Sphere')
    # 绘制单位圆边界
    circle = plt.Circle((0, 0), 1, color='green', fill=False, linewidth=2, label='Unit Sphere Boundary')
    plt.gca().add_artist(circle)
    # 图形设置
    plt.title('蒙特卡罗方法：二维投影下的四维单位球体采样点分布', fontsize=14)
    plt.xlabel('$x_1$', fontsize=12)
    plt.ylabel('$x_2$', fontsize=12)
    plt.legend(fontsize=12)
    plt.grid(True)
    plt.axis('equal')
    plt.show()

四维单位球体的蒙特卡罗估计体积: 4.934521
四维单位球体的理论体积: 4.934802

import numpy as np
import matplotlib.pyplot as plt

# -------------------- 四阶龙格-库塔方法部分 --------------------

# 定义微分方程 dy/dx = y^2 - x
def f(x, y):
    return y**2 - x

# 四阶龙格-库塔方法
def runge_kutta_4(f, x0, y0, x_end, h):
    x_values = np.arange(x0, x_end + h, h)
    y_values = np.zeros(len(x_values))
    y_values[0] = y0
    
    for i in range(1, len(x_values)):
        x_n = x_values[i-1]
        y_n = y_values[i-1]
        
        k1 = h * f(x_n, y_n)
        k2 = h * f(x_n + h/2, y_n + k1/2)
        k3 = h * f(x_n + h/2, y_n + k2/2)
        k4 = h * f(x_n + h, y_n + k3)
        
        y_values[i] = y_n + (k1 + 2*k2 + 2*k3 + k4)/6
        
    return x_values, y_values

# 定义解析解（通过符号计算得出）
# 该微分方程的解析解为 y(x) = tan(x) - x
def analytical_solution(x):
    return np.tan(x) - x

# 定义参数
x0 = 0
y0 = 0.5
x_end = 2.0  # 注意tan(x)在x=pi/2存在奇点，选择x_end < pi/2
h = 0.01

# 计算数值解
x_numeric, y_numeric = runge_kutta_4(f, x0, y0, x_end, h)

# 计算解析解
x_analytic = np.linspace(x0, x_end, 1000)
y_analytic = analytical_solution(x_analytic)

# 绘制数值解与解析解
plt.figure(figsize=(10,6))
plt.plot(x_numeric, y_numeric, 'b-', label='数值解 (RK4)')
plt.plot(x_analytic, y_analytic, 'r--', label='解析解')
plt.title(r'非线性微分方程 $\frac{dy}{dx} = y^2 - x$ 的解', fontsize=14)
plt.xlabel('$x$', fontsize=12)
plt.ylabel('$y$', fontsize=12)
plt.legend(fontsize=12)
plt.grid(True)
plt.ylim(-10,10)  # 设置y轴范围，避免tan(x)的奇点影响
plt.show()

# 打印数值解与解析解在x=2处的值
y_num_at_2 = y_numeric[-1]
y_ana_at_2 = analytical_solution(2.0)
print(f"数值解 y(2) = {y_num_at_2:.4f}")
print(f"解析解 y(2) = {y_ana_at_2:.4f}")

数值解 y(2) = -0.9957
解析解 y(2) = -4.1850