Python科学计算中的概率统计与随机模拟实战：NumPy随机数、SciPy统计分布与蒙特卡洛方法完全指南

在科学计算与数据分析领域，概率统计与随机模拟是不可或缺的核心工具。无论是金融风控中的风险评估、物理仿真中的粒子行为建模，还是机器学习中的贝叶斯推断，都离不开对随机现象的精确建模与高效模拟。Python 生态中，NumPy 提供了高性能的随机数生成引擎，SciPy 则封装了丰富的统计分布与假设检验工具。本文将从随机数生成基础出发，逐步深入到统计分布拟合、假设检验、蒙特卡洛模拟与方差缩减技术，帮助你系统掌握 Python 科学计算中的概率统计实战技能。

一、NumPy随机数生成引擎：从种子到分布

NumPy 的随机数模块经历了从

旧API到新一代

的重大演进。新API不仅性能更优，还采用了更现代的随机数生成算法。推荐始终使用

来创建生成器实例。

1.1 Generator vs legacy RandomState

旧版

使用 Mersenne Twister 算法，虽然质量不错但速度一般。新版

默认使用 PCG64 算法，在统计质量和生成速度上都有显著提升。来看一个直观的性能对比：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
import numpy as np

import time



# 旧版API

rng_legacy = np.random.RandomState(42)

start = time.perf_counter()

_ = rng_legacy.standard_normal(10_000_000)

legacy_time = time.perf_counter() - start



# 新版API

rng = np.random.default_rng(42)

start = time.perf_counter()

_ = rng.standard_normal(10_000_000)

new_time = time.perf_counter() - start



print(f"Legacy RandomState: {legacy_time:.3f}s")

print(f"New Generator:      {new_time:.3f}s")

print(f"Speedup: {legacy_time/new_time:.2f}x")

在大多数平台上，新版 Generator 生成标准正态分布随机数的速度约为旧版的 1.5-2 倍。更重要的是，

的API设计更加一致和直观，所有分布方法都遵循相同的命名规范。

1.2 常用分布的随机数生成

提供了数十种统计分布的随机数生成方法。以下是科学计算中最常用的几种：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
rng = np.random.default_rng(seed=42)



# 均匀分布 [0, 1)

uniform_samples = rng.random(1000)



# 正态分布 N(mu, sigma^2)

normal_samples = rng.normal(loc=0.0, scale=1.0, size=1000)



# 指数分布 Exp(lambda)

exp_samples = rng.exponential(scale=1.0, size=1000)  # scale = 1/lambda



# 泊松分布 Poi(lambda)

poisson_samples = rng.poisson(lam=5.0, size=1000)



# 二项分布 B(n, p)

binomial_samples = rng.binomial(n=10, p=0.3, size=1000)



# 多元正态分布

mean = [0, 0]

cov = [[1, 0.8], [0.8, 1]]

mvn_samples = rng.multivariate_normal(mean, cov, size=1000)

需要特别注意指数分布的参数：

参数是均值的倒数，即

。这是初学者最容易犯的错误之一。

1.3 可复现性与并行随机数

在科学计算中，可复现性至关重要。设置种子是最基本的方法，但在并行计算场景下，简单的种子设置会导致不同进程生成重叠的随机数序列。NumPy 提供了

机制来安全地派生独立的随机数流：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from numpy.random import SeedSequence



# 创建主种子序列

ss = SeedSequence(12345)



# 为4个并行任务派生独立的子种子

child_seeds = ss.spawn(4)

generators = [np.random.default_rng(s) for s in child_seeds]



# 每个生成器产生独立的随机数流

results = [rng.standard_normal(100) for rng in generators]



# 验证独立性：不同流之间无相关性

correlations = [np.corrcoef(results[0], results[i])[0,1] 

                for i in range(1, 4)]

print(f"Cross-stream correlations: {correlations}")

# 输出接近 [0.0, 0.0, 0.0]，确认独立性

这种

派生方式保证了即使并行任务数量变化，每个子流的随机数序列也不会重叠，是并行蒙特卡洛模拟的标准实践。

二、SciPy统计分布体系：参数估计与分布拟合

SciPy 的

模块提供了超过100种连续和离散统计分布的实现。每个分布都以面向对象的方式封装，提供统一的API接口进行PDF计算、采样、参数估计等操作。

2.1 连续分布的对象化操作

SciPy 中每个分布都是一个类，通过冻结参数可以创建具体的分布实例，这种设计非常优雅：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
from scipy import stats

import numpy as np



# 方式一：直接使用分布类（每次都要传参数）

pdf_val = stats.norm.pdf(1.96, loc=0, scale=1)

cdf_val = stats.norm.cdf(1.96, loc=0, scale=1)



# 方式二：冻结参数，创建分布对象（推荐）

dist = stats.norm(loc=100, scale=15)  # IQ分数分布



# 计算各统计量

print(f"PDF at 115: {dist.pdf(115):.6f}")

print(f"CDF at 115: {dist.cdf(115):.4f}")

print(f"SF at 115:  {dist.sf(115):.4f}")

print(f"PPF at 0.95: {dist.ppf(0.95):.2f}")

print(f"Mean: {dist.mean():.1f}, Std: {dist.std():.1f}")

冻结分布对象的写法不仅代码更简洁，还能避免反复传参导致的错误。在需要多次调用同一分布时，始终使用冻结对象。

2.2 分布拟合与参数估计

给定一组观测数据，如何确定它服从哪种分布以及分布参数？

方法使用最大似然估计（MLE）来完成这一任务：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
from scipy import stats

import numpy as np



rng = np.random.default_rng(42)



# 模拟真实数据：来自 Gamma(2, 0, 3) 分布

true_shape, true_loc, true_scale = 2.0, 0.0, 3.0

observed = rng.gamma(shape=true_shape, scale=true_scale, size=500)



# 用MLE拟合Gamma分布

fitted_shape, fitted_loc, fitted_scale = stats.gamma.fit(observed)

print(f"True:  shape={true_shape}, loc={true_loc}, scale={true_scale}")

print(f"Fitted: shape={fitted_shape:.3f}, loc={fitted_loc:.3f}, scale={fitted_scale:.3f}")



# 固定loc=0重新拟合（如果我们知道数据非负）

fitted_shape2, fitted_loc2, fitted_scale2 = stats.gamma.fit(observed, floc=0)

print(f"Fitted (loc=0): shape={fitted_shape2:.3f}, loc={fitted_loc2:.3f}, scale={fitted_scale2:.3f}")

当已知某些参数值时（如

固定位置参数），应使用

前缀的关键字参数传入固定值，这能显著提高其他参数的估计精度。

2.3 Kolmogorov-Smirnov拟合优度检验

拟合出分布后，需要检验数据是否确实来自该分布。K-S检验通过比较经验CDF与理论CDF的最大距离来判断：

1
2
3
4
5
6
7
8
9
10
# 对拟合的Gamma分布做K-S检验

ks_stat, ks_p = stats.kstest(observed, 'gamma', 

                              args=(fitted_shape2, fitted_loc2, fitted_scale2))

print(f"KS statistic: {ks_stat:.4f}")

print(f"P-value: {ks_p:.4f}")



if ks_p > 0.05:

    print("不能拒绝原假设：数据可能来自Gamma分布")

else:

    print("拒绝原假设：数据可能不来自Gamma分布")

K-S检验对位置和尺度参数敏感，如果用拟合参数做检验（而非预先指定的参数），p值会偏大。更严谨的做法是使用 Lilliefors 检验或 Anderson-Darling 检验。

三、假设检验实战：从t检验到非参数方法

假设检验是统计推断的核心工具。SciPy 提供了完整的参数检验和非参数检验工具集。选择正确的检验方法取决于数据特征和研究问题。

3.1 参数检验方法

参数检验假设数据来自特定分布（通常是正态分布）。以下是三种最常用的参数检验：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
from scipy import stats

import numpy as np



rng = np.random.default_rng(42)



# 模拟两组实验数据

group_a = rng.normal(loc=5.2, scale=1.5, size=50)

group_b = rng.normal(loc=5.8, scale=1.5, size=50)



# 1. 单样本t检验

t_stat, p_val = stats.ttest_1samp(group_a, popmean=5.0)

print(f"单样本t检验: t={t_stat:.3f}, p={p_val:.4f}")



# 2. 独立样本t检验

t_stat, p_val = stats.ttest_ind(group_a, group_b)

print(f"独立样本t检验: t={t_stat:.3f}, p={p_val:.4f}")



# 3. Welch t检验（不假设等方差，更稳健）

t_stat, p_val = stats.ttest_ind(group_a, group_b, equal_var=False)

print(f"Welch t检验: t={t_stat:.3f}, p={p_val:.4f}")



# 4. 配对t检验

before = rng.normal(loc=100, scale=15, size=30)

after = before + rng.normal(loc=-5, scale=8, size=30)

t_stat, p_val = stats.ttest_rel(before, after)

print(f"配对t检验: t={t_stat:.3f}, p={p_val:.4f}")

实际应用中，Welch t检验应作为默认选择，因为它不要求两组方差相等。经典的学生t检验在方差不等时可能给出错误的p值。

3.2 方差分析（ANOVA）

当需要比较三组或更多组的均值时，应使用方差分析而非多次t检验（避免多重比较问题）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 模拟三组实验数据

group1 = rng.normal(loc=10, scale=3, size=40)

group2 = rng.normal(loc=12, scale=3, size=40)

group3 = rng.normal(loc=11, scale=3, size=40)



# 单因素ANOVA

f_stat, p_val = stats.f_oneway(group1, group2, group3)

print(f"ANOVA: F={f_stat:.3f}, p={p_val:.4f}")



if p_val < 0.05:

    print("至少有一组均值与其他组显著不同")

    from scipy.stats import tukey_hsd

    result = tukey_hsd(group1, group2, group3)

    print(result)

ANOVA只能告诉你是否有差异，不能告诉你是哪组与哪组不同。事后检验（如Tukey HSD）用于定位具体差异来源，同时控制族错误率。

3.3 非参数检验：正态性不满足时的选择

当数据明显偏离正态分布或样本量过小时，参数检验的前提假设不再成立，需要转向非参数方法：

1
2
3
4
5
6
7
8
9
10
11
# Mann-Whitney U检验（独立样本的非参数替代）

u_stat, p_val = stats.mannwhitneyu(group_a, group_b, alternative='two-sided')

print(f"Mann-Whitney U: U={u_stat:.1f}, p={p_val:.4f}")



# Wilcoxon符号秩检验（配对样本的非参数替代）

w_stat, p_val = stats.wilcoxon(before, after)

print(f"Wilcoxon: W={w_stat:.1f}, p={p_val:.4f}")



# Kruskal-Wallis检验（ANOVA的非参数替代）

h_stat, p_val = stats.kruskal(group1, group2, group3)

print(f"Kruskal-Wallis: H={h_stat:.3f}, p={p_val:.4f}")

非参数检验不要求数据服从特定分布，而是基于秩（排名）进行比较。其统计功效通常略低于对应的参数检验，但在数据偏态严重或存在离群值时更为可靠。

四、蒙特卡洛方法：从数值积分到风险模拟

蒙特卡洛方法通过大量随机采样来近似求解确定性问题或估计随机系统的统计特征。它的理论基础是大数定律和中心极限定理，随着样本量增加，蒙特卡洛估计以 O(1/sqrt(N)) 的速率收敛。

4.1 蒙特卡洛数值积分

计算高维积分是蒙特卡洛方法最经典的应用。对于不规则区域或高维空间中的积分，蒙特卡洛方法往往比数值求积公式更高效：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
import numpy as np



rng = np.random.default_rng(42)



def estimate_pi(n_samples):

    # 在单位正方形内投点，落在单位圆内的比例约为 pi/4

    x = rng.random(n_samples)

    y = rng.random(n_samples)

    inside = (x**2 + y**2) <= 1.0

    pi_estimate = 4 * np.sum(inside) / n_samples

    pi_std = 4 * np.std(inside) / np.sqrt(n_samples)

    return pi_estimate, pi_std



for n in [1000, 10000, 100000, 1000000]:

    est, se = estimate_pi(n)

    print(f"N={n:>8d}: pi = {est:.6f} +/- {se:.6f} (error: {abs(est-np.pi):.6f})")

可以看到，样本量每增加100倍，估计误差大约减小10倍，这正是 O(1/sqrt(N)) 收敛速率的体现。

4.2 金融风险模拟：欧式期权定价

蒙特卡洛方法在金融工程中应用广泛。以欧式看涨期权定价为例，根据 Black-Scholes 模型，股票价格遵循几何布朗运动：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import numpy as np



def mc_european_call(S0, K, T, r, sigma, n_simulations=100000):

    # 蒙特卡洛模拟欧式看涨期权定价

    rng = np.random.default_rng(42)

    

    Z = rng.standard_normal(n_simulations)

    S_T = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)

    

    payoff = np.maximum(S_T - K, 0)

    discount_factor = np.exp(-r * T)

    option_price = discount_factor * np.mean(payoff)

    option_std = discount_factor * np.std(payoff) / np.sqrt(n_simulations)

    

    return option_price, option_std



# 定价示例

price, se = mc_european_call(S0=100, K=105, T=1.0, r=0.05, sigma=0.2)

print(f"期权价格: {price:.4f} +/- {se:.4f}")



# 与Black-Scholes解析解对比

from scipy.stats import norm

d1 = (np.log(100/105) + (0.05 + 0.5*0.2**2)*1.0) / (0.2*np.sqrt(1.0))

d2 = d1 - 0.2*np.sqrt(1.0)

bs_price = 100*norm.cdf(d1) - 105*np.exp(-0.05)*norm.cdf(d2)

print(f"BS解析解: {bs_price:.4f}")

print(f"误差: {abs(price - bs_price):.4f}")

10万次模拟的结果通常与解析解的误差在0.1元以内。对于没有解析解的复杂衍生品（如亚式期权、障碍期权），蒙特卡洛方法几乎是唯一可行的定价手段。

4.3 Bootstrap置信区间

Bootstrap是一种基于重采样的非参数统计方法，无需对数据分布做任何假设即可构建置信区间：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def bootstrap_ci(data, statistic=np.mean, n_bootstrap=10000, alpha=0.05):

    # 计算Bootstrap置信区间

    rng = np.random.default_rng(42)

    n = len(data)

    

    boot_stats = np.empty(n_bootstrap)

    for i in range(n_bootstrap):

        sample = rng.choice(data, size=n, replace=True)

        boot_stats[i] = statistic(sample)

    

    lower = np.percentile(boot_stats, 100 * alpha / 2)

    upper = np.percentile(boot_stats, 100 * (1 - alpha / 2))

    

    return lower, upper, boot_stats



# 示例：估计中位数的95%置信区间

rng = np.random.default_rng(42)

data = rng.exponential(scale=2.0, size=200)



lower, upper, boot_stats = bootstrap_ci(data, statistic=np.median)

print(f"中位数: {np.median(data):.3f}")

print(f"95% CI: [{lower:.3f}, {upper:.3f}]")

Bootstrap方法特别适合估计复杂统计量（如中位数、分位数、相关系数等）的置信区间，这些统计量的解析分布往往难以推导。

五、方差缩减技术：提升蒙特卡洛效率

标准蒙特卡洛的 O(1/sqrt(N)) 收敛速率意味着要将精度提高一倍，需要四倍的模拟次数。方差缩减技术通过改造采样策略，在不增加样本量的前提下显著降低估计方差。

5.1 对偶变量法

对偶变量法通过成对使用正负相关的随机数来抵消方差：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import numpy as np

from scipy.stats import norm



rng = np.random.default_rng(42)

n = 50000



# 标准蒙特卡洛

Z1 = rng.standard_normal(n)

payoff1 = np.maximum(100 * np.exp(0.05 - 0.5*0.04 + 0.2*Z1) - 105, 0)

mc_price = np.exp(-0.05) * np.mean(payoff1)

mc_var = np.exp(-0.1) * np.var(payoff1) / n



# 对偶变量法：每个Z配对使用-Z

Z = rng.standard_normal(n // 2)

S_T_pos = 100 * np.exp((0.05 - 0.5*0.04) + 0.2*Z)

S_T_neg = 100 * np.exp((0.05 - 0.5*0.04) - 0.2*Z)

payoff_pos = np.maximum(S_T_pos - 105, 0)

payoff_neg = np.maximum(S_T_neg - 105, 0)

av_price = np.exp(-0.05) * np.mean(0.5 * (payoff_pos + payoff_neg))

av_var = np.exp(-0.1) * np.var(0.5 * (payoff_pos + payoff_neg)) / (n // 2)



print(f"标准MC: 价格={mc_price:.4f}, 方差={mc_var:.6f}")

print(f"对偶变量: 价格={av_price:.4f}, 方差={av_var:.6f}")

print(f"方差缩减比: {mc_var/av_var:.1f}x")

对偶变量法在期权定价中通常能实现2-5倍的方差缩减，效果取决于收益函数的线性程度。

5.2 控制变量法

控制变量法利用与目标变量高度相关且期望已知的辅助变量来降低方差：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
# 控制变量法：用股票价格S_T作为控制变量

Z = rng.standard_normal(n)

S_T = 100 * np.exp((0.05 - 0.5*0.04) + 0.2*Z)

payoff = np.maximum(S_T - 105, 0)



# S_T的期望已知: E[S_T] = S0 * exp(r*T)

E_ST = 100 * np.exp(0.05)



# 计算最优控制系数c*

cov_matrix = np.cov(payoff, S_T)

c_star = cov_matrix[0, 1] / cov_matrix[1, 1]



# 调整后的估计

cv_estimate = np.exp(-0.05) * np.mean(payoff - c_star * (S_T - E_ST))

cv_var = np.exp(-0.1) * np.var(payoff - c_star * (S_T - E_ST)) / n



print(f"控制变量法: 价格={cv_estimate:.4f}, 方差={cv_var:.6f}")

print(f"相比标准MC方差缩减: {mc_var/cv_var:.1f}x")

控制变量法的关键在于选择与目标高度相关的控制变量。相关系数越高，方差缩减效果越显著。在期权定价中，标的资产价格是天然的控制变量，通常能实现10倍以上的方差缩减。

5.3 重要性采样

重要性采样通过改变采样分布，将更多样本集中在重要区域（如尾部事件），从而更高效地估计小概率事件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def importance_sampling_call(S0, K, T, r, sigma, n=50000):

    # 使用重要性采样的欧式看涨期权定价

    rng = np.random.default_rng(42)

    

    mu_shift = np.log(K / S0) / (sigma * np.sqrt(T))

    Z = rng.standard_normal(n) + mu_shift

    

    S_T = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)

    payoff = np.maximum(S_T - K, 0)

    

    likelihood_ratio = np.exp(-mu_shift*Z + 0.5*mu_shift**2)

    weighted_payoff = payoff * likelihood_ratio

    

    price = np.exp(-r*T) * np.mean(weighted_payoff)

    variance = np.exp(-2*r*T) * np.var(weighted_payoff) / n

    

    return price, variance



# 深虚值期权

price, var = importance_sampling_call(100, 150, 0.25, 0.05, 0.3)

print(f"IS定价: {price:.4f}, 方差: {var:.8f}")

对于深虚值期权，重要性采样可以将方差降低数十甚至数百倍。选择合适的偏移量是关键：一个经验法则是将采样分布的均值调整到使事件发生概率约为0.5的位置。

六、实战案例：工程质量可靠性评估

将上述技术综合应用于一个工程可靠性评估场景。假设我们需要评估某结构在随机荷载下的失效概率：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import numpy as np

from scipy import stats



rng = np.random.default_rng(42)



# 结构强度R ~ LogNormal(mu=5, sigma=0.2)

R = rng.lognormal(mean=5, sigma=0.2, size=100000)



# 荷载S ~ Gumbel(loc=80, scale=15)（极值分布，适合最大荷载）

S = rng.gumbel(loc=80, scale=15, size=100000)



# 安全裕度 Z = R - S

Z = R - S



# 方法1：直接蒙特卡洛估计失效概率

failure = Z < 0

pf_mc = np.mean(failure)

se_mc = np.std(failure) / np.sqrt(len(failure))

print(f"失效概率 (MC): {pf_mc:.6f} +/- {se_mc:.6f}")



# 方法2：用对数正态分布拟合R

R_shape, R_loc, R_scale = stats.lognorm.fit(R, floc=0)

print(f"R拟合参数: s={R_shape:.4f}, scale={R_scale:.4f}")



# 方法3：一阶可靠性方法（FORM近似）

mu_R, sigma_R = np.mean(R), np.std(R)

mu_S, sigma_S = np.mean(S), np.std(S)

beta = (mu_R - mu_S) / np.sqrt(sigma_R**2 + sigma_S**2)

pf_form = stats.norm.cdf(-beta)

print(f"失效概率 (FORM): {pf_form:.6f}")

print(f"可靠度指标 beta: {beta:.4f}")

这个案例展示了蒙特卡洛方法与解析近似方法的互补：蒙特卡洛适用于复杂非线性问题，FORM等方法则提供快速近似。在实际工程中，通常先用FORM快速筛选，再用蒙特卡洛精确验证。

七、性能优化与最佳实践

7.1 向量化 vs 循环

蒙特卡洛模拟中，向量化操作比逐样本循环快几个数量级：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
import numpy as np

import time



rng = np.random.default_rng(42)

n = 1_000_000



# 循环版本（极慢）

start = time.perf_counter()

result_loop = 0

Z = rng.standard_normal(n)

for z in Z:

    S = 100 * np.exp(0.03 + 0.2 * z)

    result_loop += max(S - 105, 0)

result_loop /= n

time_loop = time.perf_counter() - start



# 向量化版本（推荐）

start = time.perf_counter()

Z = rng.standard_normal(n)

S = 100 * np.exp(0.03 + 0.2 * Z)

payoff = np.maximum(S - 105, 0)

result_vec = np.mean(payoff)

time_vec = time.perf_counter() - start



print(f"循环: {result_loop:.4f}, 耗时{time_loop:.2f}s")

print(f"向量化: {result_vec:.4f}, 耗时{time_vec:.4f}s")

print(f"加速比: {time_loop/time_vec:.0f}x")

7.2 大规模模拟的内存管理

当模拟路径数极大时（如10亿条路径），不能将所有中间结果存入内存。分块处理是标准做法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def mc_chunked(n_total, chunk_size=1000000, **params):

    # 分块蒙特卡洛，支持任意大的总样本量

    rng = np.random.default_rng(42)

    S0, K, T, r, sigma = params['S0'], params['K'], params['T'], params['r'], params['sigma']

    

    total_payoff = 0.0

    total_payoff_sq = 0.0

    n_processed = 0

    

    while n_processed < n_total:

        n_chunk = min(chunk_size, n_total - n_processed)

        Z = rng.standard_normal(n_chunk)

        S_T = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)

        payoff = np.maximum(S_T - K, 0)

        

        total_payoff += np.sum(payoff)

        total_payoff_sq += np.sum(payoff**2)

        n_processed += n_chunk

    

    price = np.exp(-r*T) * total_payoff / n_total

    variance = np.exp(-2*r*T) * (total_payoff_sq/n_total - (total_payoff/n_total)**2) / n_total

    return price, np.sqrt(variance)



price, se = mc_chunked(100_000_000, S0=100, K=105, T=1.0, r=0.05, sigma=0.2)

print(f"1亿条路径: 价格={price:.4f} +/- {se:.6f}")

分块处理的关键是维护累积统计量（总和与平方和），而不是存储所有中间结果。这样即使总样本量为10亿，内存占用也仅取决于单个chunk的大小。

7.3 关键实践清单

• 始终使用

  1	default_rng(seed)

  而非

  1	np.random.xxx()

  ，确保可复现性

• 并行场景使用

  1	SeedSequence.spawn()

  派生独立种子流

• 拟合分布后用K-S检验验证拟合质量
• 假设检验优先选择Welch t检验而非经典t检验
• 小样本或偏态数据优先选择非参数检验
• 蒙特卡洛模拟始终报告标准误差，不只是点估计
• 向量化操作替代循环，性能提升可达100倍以上
• 大规模模拟使用分块处理，避免内存溢出
• 小概率事件估计使用重要性采样，避免样本浪费

总结

Python 科学计算中的概率统计工具链已经非常成熟。NumPy 的

提供了高性能、可并行化的随机数生成能力；SciPy 的

模块覆盖了从分布拟合到假设检验的完整统计推断流程；而蒙特卡洛方法配合方差缩减技术，则为复杂随机系统的建模提供了通用且强大的框架。掌握这些工具的组合使用——先用统计方法刻画数据特征，再用蒙特卡洛方法模拟系统行为，最后用假设检验验证结论——就能应对绝大多数科学计算中的概率统计问题。