如何在Python中利用最小值、最大值、均值和标准差创建适用于蒙特卡洛模拟的分布？

嘿，我完全懂你的困扰——你现在用的正态分布(np.random.normal)虽然上手简单，但它天生没有边界（理论上能取到正负无穷），还没法直接适配手里的**最小值、最大值、最可能值（众数）**这几个关键参数，这就是为啥10分位、90分位会和预期跑偏的核心原因。下面给你几个更贴合需求的方案：

1. PERT分布（优先推荐） #

PERT分布就是专门为这类有明确上下限和最可能值的场景设计的，它是Beta分布的变形，默认会根据min、max、mode拟合分布，还能通过调整参数控制集中度，进而匹配你的标准差。

代码示例 #

需要先安装scipy，然后用它自带的pert函数：

import numpy as np
from scipy.stats import pert

# 替换成你的实际参数
GCHmin = 10
GCHmax = 50
GCHmode = 30
GCHstdv = 8
num_reps = 10000

# 生成符合参数的样本
# c参数是mode在[min,max]区间内的相对位置
samples = pert.rvs(
    loc=GCHmin, 
    scale=GCHmax - GCHmin, 
    c=(GCHmode - GCHmin)/(GCHmax - GCHmin), 
    size=num_reps
)

# 验证分位数是否符合预期
print("10分位数:", np.percentile(samples, 10))
print("50分位数:", np.percentile(samples, 50))
print("90分位数:", np.percentile(samples, 90))

如果默认PERT的标准差和目标有差距，你可以调整形状参数，或者用数值方法反推合适的参数来精准匹配标准差。

2. 缩放平移的Beta分布 #

Beta分布本身定义在[0,1]区间，我们可以把它缩放平移到你的目标区间[GCHmin, GCHmax]，再通过调整形状参数a和b，同时匹配最可能值和标准差。

代码示例 #

import numpy as np
from scipy.stats import beta
from scipy.optimize import minimize

# 替换成你的实际参数
GCHmin = 10
GCHmax = 50
GCHmode = 30
GCHstdv = 8
num_reps = 10000

# 把最可能值转换到[0,1]的标准化区间
mode_standardized = (GCHmode - GCHmin) / (GCHmax - GCHmin)
# 目标标准差的标准化值
target_std_standardized = GCHstdv / (GCHmax - GCHmin)

# 定义误差函数：匹配众数和标准差的平方误差之和
def objective(params):
    a, b = params
    # 计算当前Beta分布的众数和标准差
    current_mode = (a - 1) / (a + b - 2) if a + b != 2 else 0.5
    current_std = np.sqrt(beta.stats(a, b, moments='v'))
    # 返回误差平方和
    return (current_mode - mode_standardized)**2 + (current_std - target_std_standardized)**2

# 初始猜测参数，用数值优化找到最优a和b
initial_guess = [2, 2]
result = minimize(objective, initial_guess, bounds=((1e-5, None), (1e-5, None)))
a_opt, b_opt = result.x

# 生成标准化样本，再缩放平移到目标区间
samples_standardized = beta.rvs(a_opt, b_opt, size=num_reps)
samples = GCHmin + (GCHmax - GCHmin) * samples_standardized

# 验证分位数
print("10分位数:", np.percentile(samples, 10))
print("50分位数:", np.percentile(samples, 50))
print("90分位数:", np.percentile(samples, 90))

这个方法能精准匹配你手里的四个参数，唯一的小门槛是需要用到数值优化来求解形状参数。

3. 截断正态分布 #

如果你的数据本身更接近正态分布，只是需要限制上下边界，可以用截断正态分布——把超出[GCHmin, GCHmax]的样本截断，重新采样补充。不过这个方法没法直接控制最可能值，只能调整均值和标准差来尽量贴近目标分位数。

代码示例 #

import numpy as np
from scipy.stats import truncnorm

# 替换成你的实际参数
GCHmin = 10
GCHmax = 50
GCHavg = 30
GCHstdv = 8
num_reps = 10000

# 计算截断正态分布的边界参数
lower = (GCHmin - GCHavg) / GCHstdv
upper = (GCHmax - GCHavg) / GCHstdv

# 生成截断后的样本
samples = truncnorm.rvs(lower, upper, loc=GCHavg, scale=GCHstdv, size=num_reps)

# 验证分位数
print("10分位数:", np.percentile(samples, 10))
print("50分位数:", np.percentile(samples, 50))
print("90分位数:", np.percentile(samples, 90))

总结建议 #

• 核心需求是匹配min、max、mode？优先选PERT分布，它完全贴合这类工程/风险分析场景；
• 必须精准匹配标准差？用缩放后的Beta分布，通过数值优化对齐所有参数；
• 数据本质接近正态只是需要加边界？截断正态分布最省心。

你可以根据自己的场景选合适的方法，再验证分位数是否符合预期~

内容的提问来源于stack exchange，提问作者James Black