测试矩阵
[A1,A2,...,Am] = gallery(matrixname,P1,P2,...,Pn) 生成由 matrixname 指定的一系列测试矩阵.
P1,P2,...,Pn 是单个矩阵系列要求的输入参数. 调用语法中使用的输入参数 P1,P2,...,Pn 的数目因矩阵而异.
A1,A2,...,Am — 输出系数、向量、矩阵.
A — 输出矩阵
binomial - 二项矩阵
hilbert - 对称矩阵
cauchy - 柯西矩阵
chebspec - 切比雪夫谱微分矩阵
chebvand - 切比雪夫多项式的类范德蒙类矩阵
chow - 奇异的托普利茨下黑森贝格矩阵
circul - 循环矩阵
clement - 具有零值对角线元的Clement三对角矩阵
compar - 比较矩阵
condex - 矩阵条件数估计量的反例
cycol - 其列周期性重复的矩阵
dorr - 对角占优,病态,三对角矩阵(稀疏矩阵)
dramadah - 由0和1组成的矩阵
forsythe - Forsythe矩阵或扰动Jordan块
frank - 具有病态特征值的弗兰克矩阵
gcdmat - 最大公约数矩阵
gearmat - Gear矩阵
grcar - 具有敏感特征值的托普利茨矩阵
hanowa - 其特征值位于复平面中的垂直线上的矩阵
house - Householder矩阵
invhess - 上黑森贝格矩阵的逆矩阵
invol - 对合矩阵(矩阵与其逆矩阵相同)
ipjfact - 含有阶乘元素的汉克尔矩阵
jordbloc - Jordan分块矩阵
kahan - 上梯形Kahan矩阵
kms - Kac-Murdock-Szegö 托普利茨矩阵
krylov - 克雷洛夫矩阵
lauchli - Lauchli 矩形矩阵
lehmer - Lehmer 对称正定矩阵
leslie - 由 Leslie 人口模型的出生率和成活率组成的矩阵
lesp - 含有实数敏感特征值的三对角矩阵
lotkin - Lotkin 矩阵
minij - 对称正定矩阵
moler - Moler 对称正定矩阵
neumann - 源自离散纽曼问题的奇异矩阵（稀疏矩阵）
normaldata - 从标准正态（高斯）分布中随机采样的数值组成的数组
orthog - 正交矩阵和准正交矩阵
parter - Parter矩阵
pei - Pei 矩阵
possion - 来自 Poisson 方程的块三对角矩阵（稀疏矩阵）
prolate - Prolate 矩阵
randcolu - 具有归一化列和指定的奇异值的随机矩阵
randcorr - 具有指定特征值的随机相关矩阵
randhess - 随机正交上黑森贝格矩阵
randjorth - 随机J正交矩阵
rando - 由元素 –1、0 或 1 组成的随机矩阵
randsvd - 具有预分配奇异值的随机矩阵
redheff - 由 1 和 0 组成的 Redheffer 矩阵
riemann - 与黎曼假设关联的矩阵
ris - Ris 矩阵
sampling - 具有病态整数特征值的非对称矩阵
smoke - 具有“烟环”伪谱的复矩阵
toeppd - 对称正定托普利茨矩阵
toeppen - 五对角托普利茨矩阵（稀疏矩阵）
tridiag - 三对角矩阵（稀疏矩阵）
triw - 威尔金森和其他人讨论的上三角矩阵
uniformdata - 标准均匀分布的随机采样数组成的数组
wathen - Wathen 矩阵（稀疏矩阵）
wilk - 威尔金森设计或讨论的各种矩阵
伽玛累积分布函数
p＝gamcdf（x，a）返回具有a中的形状参数的标准gamma分布的累积分布函数（cdf），在x中的值处进行评估。.
p=gamcdf（x，a，b）返回伽玛分布的cdf，其中a中的形状参数和b中的比例参数在x中的值处进行评估.
x是评估cdf的值,非负标量,非负标量数组
a是形状参数,正标量值,正标量值数组
b是比例参数,默认值是1,正标量值,正标量的数组.
如果输入参数x,a和b中的一个或多个是数组,则数组大小必须相同.在这种情况下,gamcdf将每个标量输入扩展为与数组输入大小相同的常量数组.p中的每个元素都是由a和b中的对应元素指定的分布的cdf值, 在x中的相应元素处进行评估.
示例1：等待时间分析
假设某服务台平均每分钟服务2个客户（k=2, θ=0.5分钟）
计算等待时间不超过1分钟的概率
result1 = gamcdf(1, 2, 0.5)
print(f"等待时间≤1分钟的概率: {float(result1):.4f} ({result1})")
#等待时间≤1分钟的概率: 0.5940 (0.593994150290162)
计算不同等待时间的概率
times = [0.5, 1.0, 1.5, 2.0]
for t in times:
prob = gamcdf(t, 2, 0.5)
print(f"等待时间≤{t}分钟的概率: {float(prob):.4f}")
#等待时间≤0.5分钟的概率: 0.2642
#等待时间≤1.0分钟的概率: 0.5940
#等待时间≤1.5分钟的概率: 0.8009
#等待时间≤2.0分钟的概率: 0.9084
#应用意义: 服务质量评估和资源规划
示例2：保险理赔分析
假设理赔金额服从伽马分布(k=3, θ=1000)
计算理赔金额不超过5000元的概率
result2 = gamcdf(5000, 3, 1000)
print(f"理赔金额≤5000元的概率: {float(result2):.4f}")
#理赔金额≤5000元的概率: 0.8753
计算不同金额的概率
amounts = [1000, 3000, 5000, 10000]
for amount in amounts:
prob = gamcdf(amount, 3, 1000)
print(f"理赔金额≤{amount}元的概率: {float(prob):.4f}")
#理赔金额≤1000元的概率: 0.0803
#理赔金额≤3000元的概率: 0.5768
#理赔金额≤5000元的概率: 0.8753
#理赔金额≤10000元的概率: 0.9972
#保险意义: 风险评估和保费定价
示例3：设备寿命预测
假设设备寿命服从伽马分布(k=2.5, θ=1000小时)
计算设备在2000小时内故障的概率
result3 = gamcdf(2000, 2.5, 1000)
print(f"设备在2000小时内故障的概率: {float(result3):.4f}")
#设备在2000小时内故障的概率: 0.4506
计算不同时间的可靠性
lifetimes = [500, 1000, 2000, 3000]
for time in lifetimes:
prob = gamcdf(time, 2.5, 1000)
reliability = 1 - float(prob)
print(f"设备运行{time}小时的可靠性: {reliability:.4f}")
#设备运行500小时的可靠性: 0.9626
#设备运行1000小时的可靠性: 0.8491
#设备运行2000小时的可靠性: 0.5494
#设备运行3000小时的可靠性: 0.3062
#工程意义: 预防性维护计划制定
示例4：软件缺陷发现模型
使用伽马分布模拟软件测试过程中缺陷的发现模式
k=3表示缺陷发现速率的变化，θ=50表示平均发现时间
testing_hours = [100, 200, 300, 400]
for hours in testing_hours:
prob = gamcdf(hours, 3, 50)
print(f"在{hours}小时测试内发现缺陷的概率: {float(prob):.4f}")
#在100小时测试内发现缺陷的概率: 0.3233
#在200小时测试内发现缺陷的概率: 0.7619
#在300小时测试内发现缺陷的概率: 0.9380
#在400小时测试内发现缺陷的概率: 0.9862
#软件工程意义: 测试进度评估和发布决策
示例5：极端损失概率计算
使用伽马分布模拟极端损失事件(k=2, θ=100000)
loss_thresholds = [100000, 200000, 500000, 1000000]
for loss in loss_thresholds:
prob = gamcdf(loss, 2, 100000)
print(f"损失≤{loss:,}元的概率: {float(prob):.6f}")
#损失≤100,000元的概率: 0.264241
#损失≤200,000元的概率: 0.593994
#损失≤500,000元的概率: 0.959572
#损失≤1,000,000元的概率: 0.999501
#金融意义: 风险资本计算和压力测试
示例6：信用风险建模
贷款违约时间可以用伽马分布建模
months = [6, 12, 24, 36]
for month in months:
# k=1.8, θ=18个月的平均违约时间
prob = gamcdf(month, 1.8, 18)
print(f"贷款在{month}个月内违约的概率: {float(prob):.4f}")
#贷款在6个月内违约的概率: 0.0669
#贷款在12个月内违约的概率: 0.1899
#贷款在24个月内违约的概率: 0.4494
#贷款在36个月内违约的概率: 0.6527
#银行意义: 贷款损失准备金计算
示例7：疾病潜伏期分析
疾病潜伏期常用伽马分布描述
days = [3, 7, 14, 21]
for day in days:
# k=5, θ=2表示平均潜伏期10天
prob = gamcdf(day, 5, 2)
print(f"潜伏期≤{day}天的概率: {float(prob):.4f}")
#潜伏期≤3天的概率: 0.0186
#潜伏期≤7天的概率: 0.2746
#潜伏期≤14天的概率: 0.8270
#潜伏期≤21天的概率: 0.9789
#流行病学意义: 隔离政策制定
示例8：药物治疗响应时间
药物生效时间分布
hours = [1, 2, 4, 8]
for hour in hours:
prob = gamcdf(hour, 2, 1.5)
print(f"药物在{hour}小时内生效的概率: {float(prob):.4f}")
#药物在1小时内生效的概率: 0.1443
#药物在2小时内生效的概率: 0.3849
#药物在4小时内生效的概率: 0.7452
#药物在8小时内生效的概率: 0.9694
#临床意义: 治疗方案优化
示例9：带符号参数的伽马CDF
symbolic_result = gamcdf(t, k, theta)
通用伽马CDF公式:
print(symbolic_result)
#Piecewise((lowergamma(k, t/theta)/gamma(k), t > 0), (0, True))
#数学意义: 理论分析和公式推导
示例10：特定参数下的符号结果
for k in [1, 2, 3, 4]:
result = gamcdf(t, k, 1)
print(f"k={k}, θ=1时的CDF: {result}")
#k=1, θ=1时的CDF: Piecewise((1.0 - exp(-t), t > 0), (0, True))
#k=2, θ=1时的CDF: Piecewise(((-t + 1.0*exp(t) - 1.0)*exp(-t), t > 0), (0, True))
#k=3, θ=1时的CDF: Piecewise(((-0.5*t**2 - t + 1.0*exp(t) - 1.0)*exp(-t), t > 0), (0, True))
#k=4, θ=1时的CDF: Piecewise(((-0.166666666666667*t**3 - 0.5*t**2 - t + 1.0*exp(t) - 1.0)*exp(-t), t > 0), (0, True))
#解析意义: 显示不同形状参数下的分布特性
示例11：批量概率计算
使用矩阵输入同时计算多个概率
times_matrix = [[1, 2, 3], [4, 5, 6]]
k_value = 2
theta_value = 1
时间矩阵:
for i, row in enumerate(times_matrix):
for j, time in enumerate(row):
prob = gamcdf(time, k_value, theta_value)
print(f"位置({i},{j}): 时间={time}, 概率={float(prob):.4f}")
#位置(0,0): 时间=1, 概率=0.2642
#位置(0,1): 时间=2, 概率=0.5940
#位置(0,2): 时间=3, 概率=0.8009
#位置(1,0): 时间=4, 概率=0.9084
#位置(1,1): 时间=5, 概率=0.9596
#位置(1,2): 时间=6, 概率=0.9826
#批量计算意义: 高效处理多个相似计算
示例12：贝叶斯分析中的共轭先验
伽马分布是指数族分布的共轭先验
先验分布: Gamma(k0, θ0)
k0, theta0 = 2, 1
print(f"先验分布: Gamma(k={k0}, θ={theta0})")
#先验分布: Gamma(k=2, θ=1)
后验分布计算示例
data = [1.2, 2.3, 0.8, 1.5, 2.1] # 观测数据
n = len(data)
sum_data = sum(data)
后验参数
k_posterior = k0 + n
theta_posterior = theta0 / (1 + theta0 * sum_data)
print(f"后验分布: Gamma(k={k_posterior}, θ={theta_posterior:.4f})")
#后验分布: Gamma(k=7, θ=0.1124)
#贝叶斯意义: 参数估计的不确定性量化
示例13：假设检验
检验假设: 平均等待时间是否小于2分钟
critical_value = 2.0
k_test, theta_test = 2, 0.8 # 假设的参数
p_value = float(gamcdf(critical_value, k_test, theta_test))
print(f"零假设下观测到等待时间≤{critical_value}分钟的概率: {p_value:.4f}")
#零假设下观测到等待时间≤2.0分钟的概率: 0.7127
alpha = 0.05 # 显著性水平
if p_value < alpha:
print("拒绝零假设: 平均等待时间显著小于2分钟")
else:
print("无法拒绝零假设: 没有足够证据表明平均等待时间小于2分钟")
#无法拒绝零假设: 没有足够证据表明平均等待时间小于2分钟
#统计意义: 基于分布的假设检验方法
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.stats import gamma as scipy_gamma
from sympy.stats import Gamma, cdf
def gamma_cumulative_distribution(input_str):
"""
伽马分布累积分布函数（CDF），支持符号变量和数值计算
特性：
1. 自动检测输入类型（符号/数值）
2. 符号计算使用sympy，数值计算使用scipy加速
3. 支持矩阵广播和元素级运算
输入格式：
"(x, k, theta)" 或 "(x, k)"（theta默认为1）
参数可以是标量、矩阵或包含符号的表达式
示例：
"(t, 2, 1)" -> 符号结果: 1 - (t + 1)*exp(-t)
"([[1,2],[3,4]], 2)" -> 数值矩阵结果
"""
try:
# 符号检测函数
expr = sp.sympify(input_str)
result = None
if isinstance(expr, tuple):
if len(expr) == 3:
z = expr[0]
k = expr[1] # k
theta = expr[2] # theta
elif len(expr) == 2:
z = expr[0]
k = expr[1] # k
theta = 1 # theta
else:
return f"输入错误:{input_str}"
if any(e.free_symbols for e in expr):
X = Gamma("gamma", k, theta)
result = cdf(X)(z)
else:
if k <= 0 or theta <= 0:
raise ValueError("k和theta必须>0")
result = scipy_gamma.cdf(float(z), a=float(k), scale=float(theta))
return result
except Exception as e:
return f"错误: {str(e)}"
# 测试案例
if __name__ == "__main__":
# 符号计算测试
print("符号测试:")
symbolic_input = "(t, 2, 1)"
symbolic_result = gamma_cumulative_distribution(symbolic_input)
print(f"输入: {symbolic_input}\n符号结果:", symbolic_result.simplify())
# 符号结果: Piecewise(((-t + exp(t) - 1)*exp(-t), t > 0), (0, True))
# 数值矩阵测试
print("\n数值矩阵测试:")
matrix_input = "(1, 2)"
matrix_result = gamma_cumulative_distribution(matrix_input)
print(f"输入: {matrix_input}\n数值结果:\n{np.array(matrix_result, dtype=float)}")
# 0.2642411176571153
# 混合输入测试
print("\n混合输入测试:")
mixed_input = "(x, 2, 1)"
mixed_result = gamma_cumulative_distribution(mixed_input)
print(f"输入: {mixed_input}\n混合结果:", mixed_result)
# 混合结果: Piecewise((-(x + 1)*exp(-x) + 1, x > 0), (0, True))
伽玛参数估计
[phat，pci]=gamfit（data）返回MLE和95%置信区间。pci的第一行是置信区间的下限；最后一行是上限.
示例1：设备故障时间数据分析
模拟设备故障时间数据（单位：小时）
np.random.seed(42) # 保证结果可重现
true_k, true_theta = 2.5, 800 # 真实参数：形状参数2.5，尺度参数800小时
failure_times = st.gamma.rvs(a=true_k, scale=true_theta, size=50)
print(f"故障时间数据样本（前10个）: {failure_times[:10].round(1)}")
#故障时间数据样本（前10个）: [2386.5 1575.6 1472. 1472. 4337.3 2803.2 1237.2 2453.9 2034.2 315.5]
estimates, ci = gamfit(failure_times)
print(f"真实参数: 形状k={true_k}, 尺度θ={true_theta}小时")
print(f"估计参数: 形状k={estimates[0]:.3f}, 尺度θ={estimates[1]:.3f}小时")
print(f"平均故障时间(MTTF): {estimates[0] * estimates[1]:.1f}小时")
print(f"形状参数95%置信区间: [{ci[0][0]:.3f}, {ci[0][1]:.3f}]")
print(f"尺度参数95%置信区间: [{ci[1][0]:.3f}, {ci[1][1]:.3f}]")
#真实参数: 形状k=2.5, 尺度θ=800小时
#估计参数: 形状k=4.044, 尺度θ=476.840小时
#平均故障时间(MTTF): 1928.1小时
#形状参数95%置信区间: [2.995, 6.127]
#尺度参数95%置信区间: [298.978, 660.475]
#工程意义: 设备可靠性评估和预防性维护计划
示例2：保险理赔金额分析
模拟保险理赔金额数据（单位：万元）
true_k_insurance, true_theta_insurance = 3.2, 15 # 真实参数
claim_amounts = st.gamma.rvs(a=true_k_insurance, scale=true_theta_insurance, size=200)
claim_amounts = np.maximum(claim_amounts, 0.1) # 确保最小理赔金额
estimates_insurance, ci_insurance = gamfit(claim_amounts)
print(f"理赔金额数据统计: 均值={claim_amounts.mean():.2f}万, 标准差={claim_amounts.std():.2f}万")
print(f"估计参数: 形状k={estimates_insurance[0]:.3f}, 尺度θ={estimates_insurance[1]:.3f}万")
print(f"平均理赔金额: {estimates_insurance[0] * estimates_insurance[1]:.2f}万")
print(f"形状参数95%置信区间: [{ci_insurance[0][0]:.3f}, {ci_insurance[0][1]:.3f}]")
print(f"尺度参数95%置信区间: [{ci_insurance[1][0]:.3f}, {ci_insurance[1][1]:.3f}]")
#理赔金额数据统计: 均值=47.13万, 标准差=27.04万
#估计参数: 形状k=2.945, 尺度θ=16.004万
#平均理赔金额: 47.13万
#形状参数95%置信区间: [2.499, 3.575]
#尺度参数95%置信区间: [13.070, 19.034]
#保险意义: 保费定价和准备金计算
示例3：电子元件寿命测试
模拟高温加速测试下的元件寿命数据（单位：小时）
np.random.seed(123)
true_k_component, true_theta_component = 1.8, 500
component_lifetimes = st.gamma.rvs(a=true_k_component, scale=true_theta_component, size=30)
estimates_component, ci_component = gamfit(component_lifetimes)
print(f"元件寿命数据: 最短={component_lifetimes.min():.1f}h, 最长={component_lifetimes.max():.1f}h")
print(f"估计参数: 形状k={estimates_component[0]:.3f}, 尺度θ={estimates_component[1]:.3f}小时")
#元件寿命数据: 最短=8.9h, 最长=3345.3h
#估计参数: 形状k=1.162, 尺度θ=872.768小时
计算可靠性指标
mttf = estimates_component[0] * estimates_component[1] # 平均故障时间
reliability_1000h = 1 - st.gamma.cdf(1000, estimates_component[0], scale=estimates_component[1])
print(f"平均寿命(MTTF): {mttf:.1f}小时")
print(f"1000小时可靠性: {reliability_1000h:.4f}")
print(f"形状参数95%置信区间: [{ci_component[0][0]:.3f}, {ci_component[0][1]:.3f}]")
#平均寿命(MTTF): 1014.4小时
#1000小时可靠性: 0.3829
#形状参数95%置信区间: [0.752, 2.381]
#质量工程意义: 产品寿命预测和质量改进
示例4：软件系统故障间隔时间
模拟软件系统故障间隔时间（单位：天）
true_k_software, true_theta_software = 2.0, 30
failure_intervals = st.gamma.rvs(a=true_k_software, scale=true_theta_software, size=25)
estimates_software, ci_software = gamfit(failure_intervals)
print(f"故障间隔数据: 平均间隔={failure_intervals.mean():.1f}天")
print(f"估计参数: 形状k={estimates_software[0]:.3f}, 尺度θ={estimates_software[1]:.3f}天")
print(f"平均故障间隔时间(MTBF): {estimates_software[0] * estimates_software[1]:.1f}天")
print(f"30天内无故障概率: {1 - st.gamma.cdf(30, estimates_software[0], scale=estimates_software[1]):.4f}")
#故障间隔数据: 平均间隔=50.1天
#估计参数: 形状k=2.463, 尺度θ=20.353天
#平均故障间隔时间(MTBF): 50.1天
#30天内无故障概率: 0.6986
#软件工程意义: 系统稳定性评估和发布决策
示例5：极端损失事件金额建模
模拟极端损失事件数据（单位：百万元）
np.random.seed(456)
true_k_loss, true_theta_loss = 2.5, 50 # 真实参数
extreme_losses = st.gamma.rvs(a=true_k_loss, scale=true_theta_loss, size=80)
estimates_loss, ci_loss = gamfit(extreme_losses)
print(f"极端损失数据: 次数={len(extreme_losses)}, 最大损失={extreme_losses.max():.1f}百万")
print(f"估计参数: 形状k={estimates_loss[0]:.3f}, 尺度θ={estimates_loss[1]:.3f}百万")
#极端损失数据: 次数=80, 最大损失=384.9百万
#估计参数: 形状k=3.034, 尺度θ=41.349百万
计算风险指标
var_95 = st.gamma.ppf(0.95, estimates_loss[0], scale=estimates_loss[1]) # 95% VaR
expected_shortfall = extreme_losses[extreme_losses > var_95].mean() if len(
extreme_losses[extreme_losses > var_95]) > 0 else var_95
print(f"95%风险价值(VaR): {var_95:.1f}百万")
print(f"期望短缺(Expected Shortfall): {expected_shortfall:.1f}百万")
print(f"形状参数95%置信区间: [{ci_loss[0][0]:.3f}, {ci_loss[0][1]:.3f}]")
#95%风险价值(VaR): 262.4百万
#期望短缺(Expected Shortfall): 325.8百万
#形状参数95%置信区间: [2.472, 4.018]
#风险管理意义: 资本充足率计算和压力测试
示例6：贷款违约损失分布
模拟贷款违约损失率数据（百分比）
true_k_loan, true_theta_loan = 4.0, 2.5
loss_rates = st.gamma.rvs(a=true_k_loan, scale=true_theta_loan, size=150) / 100 # 转换为小数
estimates_loan, ci_loan = gamfit(loss_rates * 100) # 保持百分比单位
print(f"违约损失率数据: 均值={loss_rates.mean() * 100:.2f}%, 标准差={loss_rates.std() * 100:.2f}%")
print(f"估计参数: 形状k={estimates_loan[0]:.3f}, 尺度θ={estimates_loan[1]:.3f}%")
print(f"平均损失率: {estimates_loan[0] * estimates_loan[1] / 100:.4f} ({estimates_loan[0] * estimates_loan[1]:.2f}%)")
#违约损失率数据: 均值=10.21%, 标准差=5.32%
#估计参数: 形状k=3.956, 尺度θ=2.582%
#平均损失率: 0.1021 (10.21%)
#信贷风险意义: 贷款损失准备金计算
示例7：疾病潜伏期分布估计
模拟疾病潜伏期数据（单位：天）
np.random.seed(789)
true_k_incubation, true_theta_incubation = 5.5, 1.8
incubation_periods = st.gamma.rvs(a=true_k_incubation, scale=true_theta_incubation, size=100)
estimates_incubation, ci_incubation = gamfit(incubation_periods)
print(f"潜伏期数据: 最短={incubation_periods.min():.1f}天, 最长={incubation_periods.max():.1f}天")
print(f"估计参数: 形状k={estimates_incubation[0]:.3f}, 尺度θ={estimates_incubation[1]:.3f}天")
#潜伏期数据: 最短=1.7天, 最长=21.6天
#估计参数: 形状k=5.401, 尺度θ=1.833天
计算重要百分位数
median_incubation = st.gamma.ppf(0.5, estimates_incubation[0], scale=estimates_incubation[1])
p95_incubation = st.gamma.ppf(0.95, estimates_incubation[0], scale=estimates_incubation[1])
print(f"潜伏期中位数: {median_incubation:.1f}天")
print(f"95%分位数: {p95_incubation:.1f}天")
print(f"形状参数95%置信区间: [{ci_incubation[0][0]:.3f}, {ci_incubation[0][1]:.3f}]")
#潜伏期中位数: 9.3天
#95%分位数: 17.8天
#形状参数95%置信区间: [4.286, 7.292]
#流行病学意义: 隔离政策制定和传播动力学建模
示例8：环境污染物浓度分析
模拟河流污染物浓度数据（单位：mg/L）
true_k_pollutant, true_theta_pollutant = 3.0, 0.8
pollutant_concentrations = st.gamma.rvs(a=true_k_pollutant, scale=true_theta_pollutant, size=60)
estimates_pollutant, ci_pollutant = gamfit(pollutant_concentrations)
regulatory_limit = 2.0 # 法规限值 2mg/L
exceedance_prob = 1 - st.gamma.cdf(regulatory_limit, estimates_pollutant[0], scale=estimates_pollutant[1])
print(f"污染物浓度数据: 均值={pollutant_concentrations.mean():.3f}mg/L")
print(f"估计参数: 形状k={estimates_pollutant[0]:.3f}, 尺度θ={estimates_pollutant[1]:.3f}mg/L")
print(f"超过法规限值({regulatory_limit}mg/L)的概率: {exceedance_prob:.4f}")
#污染物浓度数据: 均值=2.343mg/L
#估计参数: 形状k=2.381, 尺度θ=0.984mg/L
#超过法规限值(2.0mg/L)的概率: 0.5070
#环境科学意义: 环境风险评估和监管合规分析
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import scipy.stats as st
def gamma_parameter_estimates(data):
"""
估计Gamma分布的参数及置信区间，对标MATLAB的gamfit函数
参数:
data: list/array 输入数据样本
返回:
tuple: (参数估计值, 置信区间)
参数估计值格式: [形状参数, 尺度参数]
置信区间格式: [[形状参数CI], [尺度参数CI]]
如果出错返回错误信息字符串
"""
try:
# 将输入转换为numpy数组
data = np.array(data, dtype=float)
# 拟合Gamma分布参数（固定位置参数loc=0）
shape_est, _, scale_est = st.gamma.fit(data, floc=0)
# 自举法参数设置
n_bootstraps = 1000 # 通常建议1000次以上
shape_estimates = []
scale_estimates = []
# 执行自举采样
for _ in range(n_bootstraps):
sample = np.random.choice(data, size=len(data), replace=True)
shape, _, scale = st.gamma.fit(sample, floc=0)
shape_estimates.append(shape)
scale_estimates.append(scale)
# 计算95%置信区间
shape_ci = np.percentile(shape_estimates, [2.5, 97.5])
scale_ci = np.percentile(scale_estimates, [2.5, 97.5])
return [shape_est, scale_est], [shape_ci, scale_ci]
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 生成模拟数据（形状参数k=2，尺度参数θ=1）
true_k, true_theta = 2, 1
data = st.gamma.rvs(a=true_k, scale=true_theta, size=10)
# 进行参数估计
estimates, confidence_intervals = gamma_parameter_estimates(data)
# 结果展示
print(f"真实参数: 形状={true_k}, 尺度={true_theta}")
print(f"估计参数: 形状={estimates[0]:.3f}, 尺度={estimates[1]:.3f}")
print(f"形状参数95%置信区间: [{confidence_intervals[0][0]:.3f}, {confidence_intervals[0][1]:.3f}]")
print(f"尺度参数95%置信区间: [{confidence_intervals[1][0]:.3f}, {confidence_intervals[1][1]:.3f}]")
# 形状参数95%置信区间: [2.178, 8.166]
# 尺度参数95%置信区间: [0.219, 1.020]
伽玛分布的逆累积分布函数(gamma ppf = gamma cdf inverse)
x＝gaminv（p，a）返回具有形状参数a的标准伽玛分布的逆累积分布函数（icdf），该形状参数在p中的值处进行评估。.
x=gaminv（p，a，b）返回具有形状参数a和尺度参数b的伽玛分布的icdf，以p中的值进行评估。.
p是计算icdf的概率值,取值范围[0,1]
a是形状参数. 正标量值,正标量值数组
b是标度参数. 1是默认值,正标量值,正标量的数组.
示例1：金融风险价值(VaR)计算
假设损失分布服从Gamma(k=3, θ=10000)
confidence_levels = [0.90, 0.95, 0.99] # 置信水平
k_risk, theta_risk = 3.0, 10000
print(f"损失分布参数: 形状k={k_risk}, 尺度θ={theta_risk}")
#损失分布参数: 形状k=3.0, 尺度θ=10000
不同置信水平下的风险价值(VaR):
for cl in confidence_levels:
var = gaminv(cl, k_risk, theta_risk)
print(f"{cl * 100}% VaR: {var:,.0f} 元")
#90.0% VaR: 53,223 元
#95.0% VaR: 62,958 元
#99.0% VaR: 84,059 元
计算期望短缺(Expected Shortfall)
es_95 = st.gamma.expect(lambda x: x, args=(k_risk,), scale=theta_risk,
lb=gaminv(0.95, k_risk, theta_risk))
conditional_var = es_95 / (1 - 0.95) # 条件风险价值
print(f"95% 期望短缺: {conditional_var:,.0f} 元")
#95% 期望短缺: 76,017 元
#金融意义: 金融机构资本充足率计算和风险管理
示例2：保险理赔准备金计算
假设理赔分布为Gamma(k=4, θ=5000)
k_claim, theta_claim = 4.0, 5000
reserve_percentiles = [0.75, 0.90, 0.95, 0.99]
不同分位数下的理赔准备金估计:
for p in reserve_percentiles:
reserve = gaminv(p, k_claim, theta_claim)
print(f"{p * 100}% 分位数准备金: {reserve:,.0f} 元")
#75.0% 分位数准备金: 25,547 元
#90.0% 分位数准备金: 33,404 元
#95.0% 分位数准备金: 38,768 元
#99.0% 分位数准备金: 50,226 元
计算超额赔款再保险起赔点
attachment_point = gamma_inverse_ppf(f"0.95, {k_claim}, {theta_claim}")
print(f"超额赔款再保险起赔点: {attachment_point:,.0f} 元")
#超额赔款再保险起赔点: 38,768 元
#保险意义: 再保险定价和资本管理
示例3：生产过程容差限计算
假设产品尺寸变异服从Gamma(k=8, θ=0.1mm)
k_quality, theta_quality = 8.0, 0.1
tolerance_levels = [0.001, 0.01, 0.05, 0.10] # 允许的不合格品率
基于不同质量要求的过程容差限:
for alpha in tolerance_levels:
upper_limit = gaminv(1 - alpha, k_quality, theta_quality)
print(f"{alpha * 100}% 不合格率对应的上限: {upper_limit:.3f} mm")
#0.1% 不合格率对应的上限: 1.963 mm
#1.0% 不合格率对应的上限: 1.600 mm
#5.0% 不合格率对应的上限: 1.315 mm
#10.0% 不合格率对应的上限: 1.177 mm
六西格玛质量水平计算
six_sigma_limit = gaminv(0.9999966, k_quality, theta_quality)
print(f"六西格玛质量水平上限: {six_sigma_limit:.3f} mm")
#六西格玛质量水平上限: 2.756 mm
#质量工程意义: 统计过程控制和规格限设定
示例4：环境监测预警阈值
污染物浓度分布Gamma(k=2.5, θ=0.3mg/L)
k_pollutant, theta_pollutant = 2.5, 0.3
regulatory_limit = 2.0 # 法规限值
计算不同预警级别的阈值
warning_levels = [0.90, 0.95, 0.99]
不同预警级别的污染物浓度阈值:
for level in warning_levels:
threshold = gaminv(level, k_pollutant, theta_pollutant)
status = "超标" if threshold > regulatory_limit else "达标"
print(f"{level * 100}% 预警阈值: {threshold:.3f} mg/L ({status})")
#90.0% 预警阈值: 1.385 mg/L (达标)
#95.0% 预警阈值: 1.661 mg/L (达标)
#99.0% 预警阈值: 2.263 mg/L (超标)
#环境监测意义: 预警系统设计和合规管理
示例5：设备寿命保证期确定
设备寿命分布Gamma(k=2.5, θ=1000小时)
k_life, theta_life = 2.5, 1000
warranty_failure_rates = [0.01, 0.05, 0.10] # 允许的保修期内故障率
基于不同故障率要求的保修期:
for fr in warranty_failure_rates:
warranty_period = gaminv(fr, k_life, theta_life)
print(f"{fr * 100}% 故障率对应的保修期: {warranty_period:.0f} 小时")
#1.0% 故障率对应的保修期: 277 小时
#5.0% 故障率对应的保修期: 573 小时
#10.0% 故障率对应的保修期: 805 小时
计算可靠性寿命指标
b10_life = gaminv(0.10, k_life, theta_life) # B10寿命
median_life = gaminv(0.50, k_life, theta_life) # 中位寿命
print(f"B10寿命 (10%故障): {b10_life:.0f} 小时")
print(f"中位寿命: {median_life:.0f} 小时")
#B10寿命 (10%故障): 805 小时
#中位寿命: 2176 小时
#可靠性工程意义: 保修策略和产品寿命规划
示例6：预防性维护间隔优化
故障时间分布Gamma(k=3, θ=500小时)
k_maintenance, theta_maintenance = 3.0, 500
maintenance_strategies = {
"保守策略": 0.10, # 10%故障风险
"均衡策略": 0.20, # 20%故障风险
"经济策略": 0.30 # 30%故障风险
}
不同维护策略下的维护间隔:
for strategy, risk in maintenance_strategies.items():
interval = gaminv(risk, k_maintenance, theta_maintenance)
print(f"{strategy}: {interval:.0f} 小时 (故障风险: {risk * 100}%)")
#保守策略: 551 小时 (故障风险: 10.0%)
#均衡策略: 768 小时 (故障风险: 20.0%)
#经济策略: 957 小时 (故障风险: 30.0%)
计算成本最优维护间隔
preventive_cost = 1000 # 预防性维护成本
corrective_cost = 5000 # 故障修复成本
best_risk = preventive_cost / corrective_cost # 成本最优风险水平
optimal_interval = gaminv(best_risk, k_maintenance, theta_maintenance)
print(f"成本最优维护间隔: {optimal_interval:.0f} 小时 (风险: {best_risk * 100:.1f}%)")
#成本最优维护间隔: 768 小时 (风险: 20.0%)
#维护工程意义: 基于风险的维护优化
示例7：库存安全水平计算
需求分布Gamma(k=6, θ=100单位/周)
k_demand, theta_demand = 6.0, 100
service_levels = [0.90, 0.95, 0.98, 0.99] # 服务水平
lead_time = 2 # 提前期2周
average_demand = k_demand * theta_demand
print(f"平均周需求: {average_demand:.0f} 单位")
#平均周需求: 600 单位
不同服务水平下的安全库存:
for sl in service_levels:
demand_during_lead_time = gaminv(sl, k_demand * lead_time, theta_demand)
safety_stock = demand_during_lead_time - average_demand * lead_time
print(f"{sl * 100}% 服务水平: 安全库存 = {safety_stock:.0f} 单位")
#90.0% 服务水平: 安全库存 = 460 单位
#95.0% 服务水平: 安全库存 = 621 单位
#98.0% 服务水平: 安全库存 = 814 单位
#99.0% 服务水平: 安全库存 = 949 单位
#供应链意义: 库存优化和服务水平管理
示例8：供应商交货时间可靠性
交货时间分布Gamma(k=4, θ=2天)
k_delivery, theta_delivery = 4.0, 2
promised_delivery = 10 # 承诺交货时间10天
计算按时交货概率
on_time_prob = st.gamma.cdf(promised_delivery, k_delivery, scale=theta_delivery)
计算可靠交货时间
reliable_levels = [0.90, 0.95, 0.99]
print(f"承诺交货时间: {promised_delivery} 天")
print(f"按时交货概率: {on_time_prob:.1%}")
#承诺交货时间: 10 天
#按时交货概率: 73.5%
可靠交货时间估计:
for level in reliable_levels:
reliable_time = gaminv(level, k_delivery, theta_delivery)
print(f"{level * 100}% 可靠度: {reliable_time:.1f} 天")
#90.0% 可靠度: 13.4 天
#95.0% 可靠度: 15.5 天
#99.0% 可靠度: 20.1 天
#采购管理意义: 供应商绩效评估和采购决策
示例9：医疗资源需求规划
病人住院时间分布Gamma(k=3, θ=4天)
k_hospital, theta_hospital = 3.0, 4
occupancy_levels = [0.80, 0.90, 0.95] # 床位占用率目标
average_stay = k_hospital * theta_hospital
print(f"平均住院时间: {average_stay:.1f} 天")
#平均住院时间: 12.0 天
计算满足不同占用率目标的床位需求
daily_admissions = 20 # 每日入院人数
基于不同占用率目标的床位需求:
for occupancy in occupancy_levels:
# 计算住院时间分位数
stay_percentile = gaminv(occupancy, k_hospital, theta_hospital)
# 估算床位需求
bed_requirement = daily_admissions * stay_percentile
print(f"{occupancy * 100}% 占用率: {bed_requirement:.0f} 张床位")
#80.0% 占用率: 342 张床位
#90.0% 占用率: 426 张床位
#95.0% 占用率: 504 张床位
#医疗管理意义: 资源规划和效率优化
示例10：疾病治疗时间分析
治疗时间分布Gamma(k=5, θ=7天)
k_treatment, theta_treatment = 5.0, 7
recovery_probabilities = [0.50, 0.75, 0.90, 0.95]
患者康复时间预测:
for prob in recovery_probabilities:
recovery_time = gaminv(prob, k_treatment, theta_treatment)
print(f"{prob * 100}% 患者将在 {recovery_time:.0f} 天内康复")
#50.0% 患者将在 33 天内康复
#75.0% 患者将在 44 天内康复
#90.0% 患者将在 56 天内康复
#95.0% 患者将在 64 天内康复
计算治疗计划时间
treatment_plan = gamma_inverse_ppf(f"0.95, {k_treatment}, {theta_treatment}")
print(f"治疗计划时间(95%把握): {treatment_plan:.0f} 天")
#治疗计划时间(95%把握): 64 天
#临床意义: 治疗计划制定和患者沟通
示例11：置信区间计算
从样本数据估计参数后计算置信区间
np.random.seed(42)
true_k, true_theta = 2.5, 1.2
sample_data = st.gamma.rvs(a=true_k, scale=true_theta, size=100)
参数估计（使用矩估计法简化）
sample_mean = np.mean(sample_data)
sample_var = np.var(sample_data)
k_est = sample_mean ** 2 / sample_var
theta_est = sample_var / sample_mean
计算置信区间
alpha = 0.05
lower_ci = gaminv(alpha / 2, k_est, theta_est)
upper_ci = gaminv(1 - alpha / 2, k_est, theta_est)
print(f"样本均值: {sample_mean:.3f}, 样本方差: {sample_var:.3f}")
print(f"估计参数: k={k_est:.3f}, θ={theta_est:.3f}")
print(f"95% 置信区间: [{lower_ci:.3f}, {upper_ci:.3f}]")
print(f"真实均值: {true_k * true_theta:.3f} (包含在置信区间内: {lower_ci <= true_k * true_theta <= upper_ci})")
#样本均值: 2.905, 样本方差: 2.662
#估计参数: k=3.171, θ=0.916
#95% 置信区间: [0.635, 6.867]
#真实均值: 3.000 (包含在置信区间内: True)
#统计推断意义: 参数估计的不确定性量化
示例12：假设检验临界值
单侧检验的临界值计算
significance_level = 0.05
k_null, theta_null = 2.0, 1.0 # 零假设参数
critical_value = gaminv(1 - significance_level, k_null, theta_null)
print(f"显著性水平 {significance_level} 的单侧检验临界值: {critical_value:.3f}")
#显著性水平 0.05 的单侧检验临界值: 4.744
#假设检验意义: 检验决策规则的制定
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import scipy.stats as st
import sympy as sp
def gamma_inverse_ppf(input_str):
"""
Gamma分布逆累积分布函数（分位数函数），对标MATLAB的gaminv函数
参数:
p : float/array_like
累积概率值，范围需在[0,1]之间
a : float/array_like
形状参数（shape parameter），必须为正数
scale : float/array_like, 可选
尺度参数（scale parameter），默认为1，必须为正数
返回:
ndarray : 对应的分位数值，形状与输入广播后的形状一致
如果出错返回包含错误信息的字符串
示例:
>>> # 计算单个分位数
>>> gamma_inverse_ppf(0.5, 2, 1)
1.678346990016661
>>> # 计算多个分位数（向量化计算）
>>> gamma_inverse_ppf([0.1, 0.5, 0.9], 2, 1)
array([0.53181161, 1.67834699, 3.88972044])
>>> # 处理不同形状的输入
>>> gamma_inverse_ppf([[0.1], [0.9]], [[2, 3]], 1)
array([[0.53181161, 0.74008317],
[3.88972044, 5.41908035]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple):
if len(expr) == 3:
p = expr[0]
a = expr[1]
scale = expr[2]
elif len(expr) == 2:
p = expr[0]
a = expr[1]
scale = 1
else:
error = True
if isinstance(p, list):
p = sp.Matrix(p)
if isinstance(a, list):
a = sp.Matrix(a)
if isinstance(scale, list):
scale = sp.Matrix(scale)
# 将输入转换为numpy数组并进行广播
p_arr, a_arr, scale_arr = np.broadcast_arrays(
np.asarray(p, dtype=float),
np.asarray(a, dtype=float),
np.asarray(scale, dtype=float)
)
# 参数有效性检查
if np.any((p_arr < 0) | (p_arr > 1)):
raise ValueError("概率值p必须在[0, 1]范围内")
if np.any(a_arr <= 0):
raise ValueError("形状参数a必须为正数")
if np.any(scale_arr <= 0):
raise ValueError("尺度参数scale必须为正数")
# 初始化结果数组
result = np.empty_like(p_arr, dtype=float)
# 向量化计算分位数
for index in np.ndindex(p_arr.shape):
result[index] = st.gamma.ppf(
q=p_arr[index],
a=a_arr[index],
scale=scale_arr[index]
)
if result.size == 1:
result = float(result)
else:
result = sp.Matrix(result.tolist())
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例测试
if __name__ == "__main__":
# 测试用例1：标量输入
print("标量输入测试:")
print(gamma_inverse_ppf("0.5, 3, 5"))
# 13.370301568617794
# 测试用例2：向量输入
print("\n向量输入测试:")
print(gamma_inverse_ppf("[0.1, 0.5, 0.9], 2"))
# Matrix([[0.531811608389612],
# [1.67834699001666],
# [3.88972016986743]])
# 测试用例3：广播测试
print("\n广播测试:")
print(gamma_inverse_ppf("[[0.1], [0.9]], [[2, 3]]"))
# Matrix([[0.531811608389612, 1.10206532824932],
# [3.88972016986743, 5.32232033783421]])
伽玛分布的负对数似然函数
nlogL=gamlike（params，data）返回给定数据的参数params的gamma对数似然的负值. params（1）=A,形状参数,params（2）=B，比例参数。params中的参数必须全部为正.
示例1：保险索赔金额建模
模拟保险索赔数据（单位：万元）
np.random.seed(42)
claims_data = np.random.gamma(shape=2.5, scale=1.2, size=100)
假设我们已经通过某种方法得到了参数估计
estimated_params = [2.3, 1.1] # [形状参数, 尺度参数]
计算负对数似然值（衡量模型拟合优度）
nll = gamlike(estimated_params, claims_data)
print(f"索赔数据负对数似然值: {nll}")
#索赔数据负对数似然值: 183.46245646722977
比较不同参数模型的拟合效果
alt_params1 = [2.0, 1.0]
alt_params2 = [3.0, 1.5]
nll1 = gamlike(alt_params1, claims_data)
nll2 = gamlike(alt_params2, claims_data)
print(f"参数{alt_params1}的负对数似然: {nll1}")
print(f"参数{alt_params2}的负对数似然: {nll2}")
#参数[2.0, 1.0]的负对数似然: 199.76099175588814
#参数[3.0, 1.5]的负对数似然: 203.11551575829495
#负对数似然值越小，说明模型拟合越好
示例2：设备寿命分析
模拟设备故障时间数据（单位：小时）
failure_times = [125, 89, 156, 203, 78, 145, 167, 198, 112, 134]
不同维护策略下的参数估计
strategy_a = [2.8, 45] # 策略A：定期维护
strategy_b = [3.2, 38] # 策略B：预测性维护
nll_a = gamlike(strategy_a, failure_times)
nll_b = gamlike(strategy_b, failure_times)
print(f"策略A(定期维护)负对数似然: {nll_a}")
print(f"策略B(预测性维护)负对数似然: {nll_b}")
#策略A(定期维护)负对数似然: 54.75534678188527
#策略B(预测性维护)负对数似然: 54.40380198332201
if nll_a < nll_b:
print("策略A的模型拟合更好")
else:
print("策略B的模型拟合更好")
#策略B的模型拟合更好
示例3：降雨量建模
某地区月降雨量数据（单位：mm）
monthly_rainfall = [45.2, 78.9, 123.4, 156.7, 89.1, 67.8,
34.5, 56.7, 98.9, 112.3, 145.6, 87.9]
不同季节的参数（春、夏、秋、冬）
seasonal_params = [
[2.5, 35], # 春季
[3.2, 45], # 夏季
[2.8, 40], # 秋季
[2.0, 25] # 冬季
]
计算每个季节模型对全年数据的拟合程度
for i, params in enumerate(seasonal_params):
season_name = ["春季", "夏季", "秋季", "冬季"][i]
nll = gamlike(params, monthly_rainfall)
print(f"{season_name}模型负对数似然: {nll:.2f}")
#春季模型负对数似然: 61.76
#夏季模型负对数似然: 64.34
#秋季模型负对数似然: 61.98
#冬季模型负对数似然: 68.03
示例4：住院时间分析
不同科室病人住院天数
surgery_days = [5, 7, 3, 8, 6, 4, 9, 5, 7, 6]
internal_days = [10, 8, 12, 9, 11, 7, 13, 10, 9, 11]
两个科室的参数估计
surgery_params = [6.0, 1.2] # 外科
internal_params = [10.0, 1.1] # 内科
计算各自模型的拟合优度
nll_surgery = gamlike(surgery_params, surgery_days)
nll_internal = gamlike(internal_params, internal_days)
print(f"外科模型负对数似然: {nll_surgery:.2f}")
print(f"内科模型负对数似然: {nll_internal:.2f}")
#外科模型负对数似然: 21.54
#内科模型负对数似然: 22.62
交叉验证：用外科模型拟合内科数据
cross_nll = gamlike(surgery_params, internal_days)
print(f"外科模型拟合内科数据的负对数似然: {cross_nll:.2f}")
#外科模型拟合内科数据的负对数似然: 27.79
#交叉验证值明显更大，说明模型不适用于不同分布的数据
示例5：参数优化过程
模拟数据
data = np.random.gamma(shape=3.0, scale=2.0, size=50).tolist()
测试不同的形状参数
shape_values = [1.5, 2.0, 2.5, 3.0, 3.5, 4.0]
scale_fixed = 2.0
固定尺度参数=2.0，测试不同形状参数的拟合效果:
for shape in shape_values:
params = [shape, scale_fixed]
nll = gamlike(params, data)
print(f"形状参数={shape}, 负对数似然={nll:.2f}")
#形状参数=1.5, 负对数似然=148.73
#形状参数=2.0, 负对数似然=132.18
#形状参数=2.5, 负对数似然=123.83
#形状参数=3.0, 负对数似然=121.67
#形状参数=3.5, 负对数似然=124.48
#形状参数=4.0, 负对数似然=131.44
找到最佳参数（最小负对数似然）
best_shape = min(shape_values, key=lambda s: gamlike([s, 2.0], data))
print(f"\n最佳形状参数: {best_shape}")
#最佳形状参数: 3.0
示例6：异常检测
正常数据（设备正常运行时的能耗）
normal_energy = [25.1, 26.3, 24.8, 25.7, 26.1, 25.4, 25.9, 24.7, 26.2, 25.5]
拟合正常数据的Gamma分布参数
normal_params = [25.0, 1.0] # 假设已通过MLE估计得到
测试数据点（包含正常和异常）
test_points = [25.3, 26.1, 35.8, 24.9, 40.2, 25.6]
测试点的负对数似然（值越大越可能是异常）:
for point in test_points:
nll = gamlike(normal_params, [point]")
print(f"能耗{point}kW: 负对数似然={nll:.2f}")
#能耗25.3kW: 负对数似然=2.55
#能耗26.1kW: 负对数似然=2.60
#能耗35.8kW: 负对数似然=4.71
#能耗24.9kW: 负对数似然=2.53
#能耗40.2kW: 负对数似然=6.33
#能耗25.6kW: 负对数似然=2.56
设置阈值
threshold = 5.0
anomalies = [pt for pt in test_points if gamlike(normal_params, [pt]) > threshold]
print(f"\n检测到的异常值: {anomalies}")
#检测到的异常值: [40.2]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import scipy.stats as st
def gamma_neg_log_like(params, data):
"""
计算Gamma分布的负对数似然，对标MATLAB的gamlike函数
参数:
params : tuple/list/array_like
包含形状参数(shape)和尺度参数(scale)的数组，格式为[shape, scale]
data : list/array_like
观测数据样本，所有值必须为正数
返回:
float: 负对数似然值
如果出错返回包含错误信息的字符串
示例:
>>> # MATLAB等效调用：gamlike([2, 1], [1.2, 0.8, 2.3])
>>> params = [2, 1]
>>> data = [1.2, 0.8, 2.3]
>>> gamma_neg_log_like(params, data)
3.876492512616387
"""
try:
# 参数类型转换和验证
params = np.asarray(params, dtype=float)
data = np.asarray(data, dtype=float)
# 检查参数维度
if params.size != 2:
raise ValueError("参数必须包含两个元素：[形状参数, 尺度参数]")
# 解包参数
shape, scale = params[0], params[1]
# 参数有效性检查
if shape <= 0 or scale <= 0:
raise ValueError("形状参数和尺度参数必须大于0")
if np.any(data <= 0):
raise ValueError("所有数据值必须为正数")
# 计算对数似然
log_pdf = st.gamma.logpdf(data, a=shape, scale=scale)
return -np.sum(log_pdf)
except Exception as e:
return f"错误: {e}"
# 示例测试
if __name__ == "__main__":
# 测试用例1：正常情况
data = [1.2, 0.8, 2.3, 1.5]
params = [2.0, 1.0]
print(f"测试1结果: {gamma_neg_log_like(params, data)}")
# 4.602447763476986
# 测试用例2：多参数维度测试
matrix_params = [[2.0, 1.0], [3.0, 0.5]] # 两组不同参数
print(f"\n测试4结果: {[gamma_neg_log_like(p, data) for p in matrix_params]}")
# [4.602447763476986, 3.6597180824744107]
Y = gamma(X) 返回在X的元素处计算的gamma函数值.
X是标量,向量,矩阵,多维数组
示例1：概率论与组合数学
1.1 计算阶乘（Gamma(n+1) = n!）
for n in [5, 10, 15]:
factorial_gamma = gamma(n + 1)
print(f"Gamma({n + 1}) = {factorial_gamma} (等价于 {n}! = {np.math.factorial(n)})")
#Gamma(6) = (119.99999999999987+0j) (等价于 5! = 120)
#Gamma(11) = (3628800.0000000084+0j) (等价于 10! = 3628800)
#Gamma(16) = (1307674368000.003+0j) (等价于 15! = 1307674368000)
1.2 计算二项式系数 C(n, k) = n!/(k!(n-k)!) = Gamma(n+1)/(Gamma(k+1)*Gamma(n-k+1))
n, k = 10, 3
binomial_coeff = (gamma(n + 1) /
(gamma(k + 1) * gamma(n - k + 1)))
print(f"C({n}, {k}) = {binomial_coeff}")
#C(10, 3) = 120.00000000000061
1.3 计算半整数的Gamma值（在概率分布中常见）
half_integers = [0.5, 1.5, 2.5, 3.5]
for x in half_integers:
result = gamma(x)
print(f"Gamma({x}) = {result}")
#Gamma(0.5) = 1.7724538509055134
#Gamma(1.5) = 0.8862269254527566
#Gamma(2.5) = 1.3293403881791352
#Gamma(3.5) = 3.323350970447838
示例2：统计学中的分布函数
2.1 Gamma分布的归一化常数
def gamma_pdf(x, alpha, beta):
"""Gamma分布的概率密度函数"""
return (beta ** alpha / gamma(alpha)) * x ** (alpha - 1) * np.exp(-beta * x)
alpha, beta = 3, 2
x_values = np.linspace(0.1, 5, 100)
pdf_values = [gamma_pdf(x, alpha, beta) for x in x_values]
print(f"Gamma分布PDF在alpha={alpha}, beta={beta}时的归一化常数: 1/Gamma({alpha}) = {1 / gamma(str(alpha))}")
#Gamma分布PDF在alpha=3, beta=2时的归一化常数: 1/Gamma(3) = 0.5000000000000009
2.2 学生t分布的归一化常数
def student_t_constant(nu):
"""学生t分布的归一化常数"""
return gamma((nu+1)/2) / (np.sqrt(nu * np.pi) * gamma(nu / 2))
degrees_freedom = [1, 5, 10, 30]
学生t分布的归一化常数:
for nu in degrees_freedom:
const = student_t_constant(nu)
print(f"自由度 ν={nu}: 归一化常数 = {const}")
#自由度 ν=1: 归一化常数 = 0.3183098861837912
#自由度 ν=5: 归一化常数 = 0.37960668982249435
#自由度 ν=10: 归一化常数 = 0.3891083839660311
#自由度 ν=30: 归一化常数 = 0.39563218489409696
示例3：工程与物理学应用
3.1 计算球体的超体积（n维球体体积公式包含Gamma函数）
def n_sphere_volume(n, R):
"""计算n维球体的体积"""
return (np.pi ** (n / 2) * R ** n) / gamma(n / 2 + 1)
n维球体体积计算:
dimensions = [2, 3, 4, 5] # 2维圆，3维球，4维超球等
radius = 1.0
for n in dimensions:
volume = n_sphere_volume(n, radius)
print(f"{n}维球体体积 (R=1): {volume:.4f}")
#2维球体体积 (R=1): 3.1416
#3维球体体积 (R=1): 4.1888
#4维球体体积 (R=1): 4.9348
#5维球体体积 (R=1): 5.2638
3.2 量子力学中的波函数归一化
def hydrogen_wavefunction_constant(n, l):
"""氢原子波函数的归一化常数"""
# 简化公式，实际公式更复杂
return np.sqrt(gamma(n - l) / gamma(n + l + 1))
氢原子波函数归一化常数:
quantum_numbers = [(1, 0), (2, 0), (2, 1), (3, 1)]
for n, l in quantum_numbers:
const = hydrogen_wavefunction_constant(n, l)
print(f"主量子数n={n}, 角量子数l={l}: 归一化常数 ≈ {const:.4f}")
#主量子数n=1, 角量子数l=0: 归一化常数 ≈ 1.0000
#主量子数n=2, 角量子数l=0: 归一化常数 ≈ 0.7071
#主量子数n=2, 角量子数l=1: 归一化常数 ≈ 0.4082
#主量子数n=3, 角量子数l=1: 归一化常数 ≈ 0.2041
示例4：金融数学应用
4.1 Black-Scholes模型中的Gamma（二阶导数）
def black_scholes_gamma(S, K, T, r, sigma):
"""计算期权定价的Gamma值"""
from scipy.stats import norm
d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
return norm.pdf(d1) / (S * sigma * np.sqrt(T))
虽然这里用的是正态分布，但Gamma函数在更复杂的衍生品定价中出现
4.2 风险度量中的极端值理论
def extreme_value_constant(shape_param):
"""极值分布中的常数项（包含Gamma函数）"""
if shape_param != 0:
return gamma(1 - shape_param)
else:
return 1.0 # Gumbel分布情况
shape_params = [-0.5, 0, 0.3]
极值分布中的Gamma函数应用:
for xi in shape_params:
const = extreme_value_constant(xi)
print(f"形状参数 ξ={xi}: 相关常数 = {const}")
#形状参数 ξ=-0.5: 相关常数 = 0.8862269254527566+
#形状参数 ξ=0: 相关常数 = 1.0
#形状参数 ξ=0.3: 相关常数 = 1.2980553326475577
示例5：数值分析与特殊函数
5.1 计算Beta函数 B(x,y) = Gamma(x)Gamma(y)/Gamma(x+y)
def beta_function(x, y):
"""Beta函数计算"""
return (gamma(x) * gamma(y)) / gamma(x + y)
Beta函数计算:
pairs = [(2, 3), (0.5, 0.5), (3, 4)]
for x, y in pairs:
beta_val = beta_function(x, y)
print(f"B({x}, {y}) = {beta_val:.4f}")
#B(2, 3) = 0.0833
#B(0.5, 0.5) = 3.1416
#B(3, 4) = 0.0167
5.2 不完全Gamma函数的相关计算
def incomplete_gamma_ratio(a, x):
"""计算Gamma(a,x)/Gamma(a)的近似（用于不完全Gamma函数）"""
# 简化版本，实际需要更复杂的数值方法
return 1 - np.exp(-x) * x ** a / (gamma(a) * a)
不完全Gamma函数相关计算:
parameters = [(2, 1), (3, 2), (0.5, 0.1)]
for a, x in parameters:
ratio = incomplete_gamma_ratio(a, x)
print(f"P({a}, {x}) ≈ {ratio:.4f} (正则化不完全Gamma函数)")
#P(2, 1) ≈ 0.8161+0.0000j (正则化不完全Gamma函数)
#P(3, 2) ≈ 0.8196+0.0000j (正则化不完全Gamma函数)
#P(0.5, 0.1) ≈ 0.6771+0.0000j (正则化不完全Gamma函数)
示例6：符号计算与数学恒等式
6.1 反射公式: Gamma(z)Gamma(1-z) = π/sin(πz)
z_values = [0.25, 0.33, 0.5, 0.75]
Gamma函数反射公式验证:
for z in z_values:
left_side = gamma(z) * gamma(1 - z)
right_side = np.pi / np.sin(np.pi * z)
print(f"z={z}: Gamma({z})*Gamma({1 - z}) = {left_side:.6f}")
print(f" π/sin(π{z}) = {right_side:.6f}")
print(f" 差值: {abs(left_side - right_side):.2e}")
#z=0.25: Gamma(0.25)*Gamma(0.75) = 4.442883+0.000000j
π/sin(π0.25) = 4.442883
差值: 8.88e-16
#z=0.33: Gamma(0.33)*Gamma(0.6699999999999999) = 3.649866+0.000000j
π/sin(π0.33) = 3.649866
差值: 4.44e-15
#z=0.5: Gamma(0.5)*Gamma(0.5) = 3.141593+0.000000j
π/sin(π0.5) = 3.141593
差值: 8.88e-15
#z=0.75: Gamma(0.75)*Gamma(0.25) = 4.442883+0.000000j
π/sin(π0.75) = 4.442883
差值: 8.88e-16
6.2 递推关系: Gamma(z+1) = z * Gamma(z)
recurrence = gamma(z + 1) - symbolic_z * gamma(z)
print(f"Gamma(z+1) - z*Gamma(z) = {recurrence} (应该为0)")
#Gamma(z+1) - z*Gamma(z) = 0 (应该为0)
示例7：实际数据分析应用
7.1 拟合Gamma分布参数的最大似然估计
def gamma_mle_equations(data):
"""Gamma分布参数估计的方程（包含Gamma函数）"""
n = len(data)
mean_log = np.mean(np.log(data))
mean_data = np.mean(data)
# 最大似然方程涉及Gamma函数的导数（Digamma函数）
# 这里简化显示Gamma函数的作用
print("Gamma分布在参数估计中的作用:")
print(f"样本大小: {n}")
print(f"数据均值: {mean_data:.4f}")
print(f"对数均值: {mean_log:.4f}")
print("MLE方程涉及Gamma函数的对数导数")
# 模拟等待时间数据（通常服从Gamma分布）
waiting_times = np.random.gamma(shape=2.5, scale=1.0, size=100)
gamma_mle_equations(waiting_times)
#样本大小: 100
#数据均值: 2.4858
#对数均值: 0.7186
#MLE方程涉及Gamma函数的对数导数
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import gamma
def gamma_compute(input_str):
"""
对标 MATLAB 的 gamma 函数，支持标量、矩阵和符号计算。
参数:
input_str: 输入的数学表达式字符串。
返回:
计算结果或错误信息字符串。
示例:
>>> gamma_compute("5")
24.0
>>> gamma_compute("Matrix([[1, 2], [3, 4]])")
Matrix([[1, 1], [2, 6]])
>>> gamma_compute("x + 1")
x*gamma(x)
>>> gamma_compute("invalid")
'错误: ...'
"""
try:
expr = sp.sympify(input_str)
# 数值处理分支
if expr.is_number:
z = complex(expr)
return gamma(z)
# 符号表达式处理分支
elif expr.free_symbols:
return sp.expand_func(sp.gamma(expr))
# 未知类型处理
else:
return f"无法处理的输入类型: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例测试
if __name__ == "__main__":
# 标量测试
print(gamma_compute("2.5"))
# (1.3293403881791352+0j)
# 符号测试
print(gamma_compute("x + 1"))
# x*gamma(x)
正则化不完全gamma函数
Y = gammainc(X,A) 返回在 X 和 A 的元素处计算的正则化下不完全 gamma 函数。X 和 A 必须都为实数，A 必须为非负值.
Y = gammainc(X,A,type) 返回正则化下/上不完全 gamma 函数。type 的选项是 'lower'（默认值）和 'upper'.
X是标量,向量,矩阵,多维数组
A是标量,向量,矩阵. 多维数组, A的元素必须为非负实数.X和A必须大小相同,或者其中之一必须为标量.
type是正则化不完全gamma函数的类型,upper或lower(默认)
示例1. 通信系统中的误码率计算
Q函数与上不完全伽马函数的关系：Q(x) = 0.5 * Γ(0.5, x²/2)
def q_function(x):
"""计算Q函数（标准正态分布的右尾概率）"""
return 0.5 * gammainc(0.5, x**2/2, upper)
测试Q函数
print("Q(1) =", q_function(1))
print("Q(2) =", q_function(2))
#Q(1) = 0.282094791773878*uppergamma(0.5, x**2/2)
#Q(2) = 0.282094791773878*uppergamma(0.5, x**2/2)
BPSK调制的误码率
def bpsk_ber(snr_db):
"""计算BPSK调制的误码率"""
snr_linear = 10 ** (snr_db / 10)
return q_function(sp.sqrt(2 * snr_linear))
print("BPSK BER at SNR=10dB:", bpsk_ber(10))
#BPSK BER at SNR=10dB: 0.282094791773878*uppergamma(0.5, x**2/2)
示例2. 统计学中的卡方分布
卡方分布的累积分布函数
def chi_squared_cdf(x, k):
"""自由度为k的卡方分布CDF"""
# P(X ≤ x) = γ(k/2, x/2) / Γ(k/2)
return gammainc(k/2, x/2, lower)
计算卡方分布的p值
def chi_squared_p_value(x, k):
"""卡方检验的p值"""
return 1 - chi_squared_cdf(x, k)
假设检验示例
chi2_statistic = 15.0 # 卡方统计量
degrees_freedom = 10 # 自由度
p_value = chi_squared_p_value(chi2_statistic, degrees_freedom)
print(f"卡方统计量={chi2_statistic}, 自由度={degrees_freedom}")
print(f"p值 = {p_value}")
#卡方统计量=15.0, 自由度=10
#p值 = 0.1320618562877206
临界值计算（近似）
critical_value_95 = 18.3 # 95%置信水平的临界值
if chi2_statistic > critical_value_95:
print("拒绝原假设")
else:
print("不能拒绝原假设")
#不能拒绝原假设
示例3. 矩阵输入的批量计算
创建测试矩阵
a_matrix = [[1, 2], [3, 4]]
x_matrix = [[0.5, 1.5], [2.5, 3.5]]
矩阵形式的下不完全伽马函数
result_lower = gammainc(a_matrix, x_matrix, lower)
下不完全伽马函数（矩阵）:
print(result_lower)
#[[0.39346934 0.4421746 ]
[0.45618688 0.46336733]]
矩阵形式的上不完全伽马函数
result_upper = gammainc(a_matrix, x_matrix, upper)
上不完全伽马函数（矩阵）:
print(result_upper)
#[[0.60653066 0.5578254 ]
[0.54381312 0.53663267]]
示例4. 符号计算应用
符号计算示例
基本符号表达式
symbolic_result = gammainc(a, x, lower)
print(f"Γ(a, x)/Γ(a) = {symbolic_result}")
#Γ(a, x)/Γ(a) = lowergamma(a, x)/gamma(a)
特定参数情况
specific_case = gammainc(n+1, x, lower)
print(f"Γ(n+1, x)/Γ(n+1) = {specific_case}")
#Γ(n+1, x)/Γ(n+1) = lowergamma(n + 1, x)/gamma(n + 1)
半整数参数（常见于物理问题）
half_integer = gammainc(3/2, x, upper)
print(f"Γ(3/2, x)/Γ(3/2) = {half_integer}")
#Γ(3/2, x)/Γ(3/2) = 1.12837916709551*uppergamma(3/2, x)
示例5. 金融风险管理中的应用
伽马分布在金融风险建模中的应用
def gamma_var(alpha, beta, confidence_level=0.95):
"""基于伽马分布的风险价值计算"""
# VaR = F⁻¹(confidence_level)，其中F是伽马分布的CDF
# 使用不完全伽马函数求逆（这里简化处理）
return gammainc(alpha, beta*x, lower)
信用风险模型中的违约概率
def default_probability(pd, lgd, exposure):
"""简化版的信用风险计算"""
# 使用伽马分布建模违约损失
return gammainc(2, pd*exposure, upper) * lgd
pd = 0.02 # 违约概率
lgd = 0.6 # 违约损失率
exposure = 1000000 # 风险暴露
prob = default_probability(pd, lgd, exposure)
print(f"预期信用损失: {prob}")
#预期信用损失: 0.0
示例6. 物理学中的应用
黑体辐射的普朗克定律相关计算
def planck_integral(nu, T):
"""计算黑体辐射的谱辐射出射度（简化版）"""
# 使用不完全伽马函数计算特定频率范围的辐射
h = 6.626e-34 # 普朗克常数
k = 1.381e-23 # 玻尔兹曼常数
x = h * nu / (k * T)
return gammainc(f4, x, upper) * (8 * np.pi * k ** 4 * T ** 4) / (h ** 3 * 3.0e8 ** 3)
计算太阳表面温度下的辐射
T_sun = 5778 # 开尔文
nu_visible = 5e14 # 可见光频率
radiation = planck_integral(nu_visible, T_sun)
print(f"太阳在可见光频段的辐射出射度: {radiation:.2e} W/m²/Hz")
#太阳在可见光频段的辐射出射度: 5.25e-02 W/m²/Hz
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.special as sc
def gamma_incomplete_function(input_str):
"""
对标 MATLAB 的 gammainc 函数，计算正则化的下/上不完全伽马函数。
参数:
input_str: 输入的数学表达式字符串，格式为元组 (x, a, 'type')，其中 type 可选 'lower' 或 'upper'。
返回:
计算结果或错误信息字符串。
示例:
>>> gamma_incomplete_function("(2, 5, lower)")
0.052653017343711
>>> gamma_incomplete_function("(Matrix([[1,2]]), Matrix([[3,4]]), upper)")
Matrix([[0.849145, 0.238536]])
>>> gamma_incomplete_function("(x, a, lower)")
lowergamma(a, x)/gamma(a)
"""
try:
expr = sp.sympify(input_str, evaluate=False)
error = False
result = None
if isinstance(expr, tuple):
# 解析参数
if len(expr) == 3:
x, a, gammainc_type = expr
gammainc_type = str(gammainc_type)
elif len(expr) == 2:
x, a = expr
gammainc_type = 'lower'
else:
error = True
return f"输入错误: 参数数量错误，期望2或3个参数，实际{len(expr)}个"
# 检查类型参数有效性
if gammainc_type not in ['lower', 'upper']:
error = True
return f"输入错误: 类型参数必须为 'lower' 或 'upper'，实际为 {gammainc_type}"
if gammainc_type == "lower":
if any(e.free_symbols for e in [x, a]):
result = sp.lowergamma(a, x) / sp.gamma(a)
elif all(e.is_number for e in [x, a]):
result = sc.gammainc(float(a), float(x))
else:
if any(e.free_symbols for e in [x, a]):
result = sp.uppergamma(a, x) / sp.gamma(a)
elif all(e.is_number for e in [x, a]):
result = sc.gammaincc(float(a), float(x))
return result
else:
return f"输入错误: 输入必须是元组"
except Exception as e:
return f"错误: {e}"
# 示例测试
if __name__ == "__main__":
# 标量测试
print(gamma_incomplete_function("2, 5, lower"))
# 0.052653017343711125
print(gamma_incomplete_function("(5, 2, upper)"))
# 0.04042768199451279
# 符号测试
print(gamma_incomplete_function("(x, a, lower)"))
# lowergamma(a, x)/gamma(a)
正则化不完全gamma函数的逆函数
X = gammaincinv(Y,A)返回在Y和A元素处计算的正则化下不完全gamma函数的逆函数,满足Y = gammainc(X,A). Y和A都必须为实数. Y的元素必须在闭区间[0,1]内,A必须为非负值.
X = gammaincinv(Y,A,type) 返回正则化下或上不完全gamma函数的逆函数.type的选项是'lower'（默认值）和'upper'.
Y是标量,向量,矩阵,多维数组
A是标量,向量,矩阵,多维数组. A的元素必须为非负实数.Y和A必须大小相同,或者其中之一必须为标量.
type是逆不完全gamma函数的类型,upper或lower(默认)
示例1. 统计学中的置信区间计算
卡方分布的临界值计算
def chi_squared_critical_value(confidence_level, degrees_freedom):
"""
计算卡方分布的临界值（用于假设检验）
"""
# P(X ≤ x) = confidence_level => x = gammaincinv(confidence_level, df/2) * 2
return 2 * gammaincinv(confidence_level, degrees_freedom/2)
计算95%置信水平的卡方临界值
df = 10
critical_value = chi_squared_critical_value(0.95, df)
print(f"自由度为{df}的卡方分布在95%置信水平的临界值: {critical_value:.4f}")
#自由度为10的卡方分布在95%置信水平的临界值: 18.3070
不同自由度下的临界值比较
degrees = [5, 10, 15, 20]
conf_levels = [0.90, 0.95, 0.99]
不同自由度和置信水平的卡方临界值:
for df in degrees:
row = []
for cl in conf_levels:
cv = chi_squared_critical_value(cl, df)
row.append(f"{cv:.3f}")
print(f"自由度{df:2d}: {', '.join(row)}")
#自由度 5: 9.236, 11.070, 15.086
#自由度10: 15.987, 18.307, 23.209
#自由度15: 22.307, 24.996, 30.578
#自由度20: 28.412, 31.410, 37.566
#设备寿命的中位数: 267.41 小时
示例2. 可靠性工程中的寿命预测
伽马分布的分位数计算
def gamma_quantile(probability, shape_param, scale_param=1):
"""
计算伽马分布的分位数（逆CDF）
F(x) = P(X ≤ x) = γ(α, x/β)/Γ(α)
=> x = β * gammaincinv(P, α)
"""
return scale_param * gammaincinv(probability, shape_param)
设备寿命的中位数估计
shape = 3.0 # 形状参数
scale = 100 # 尺度参数（平均寿命）
median_life = gamma_quantile(0.5, shape, scale)
print(f"设备寿命的中位数: {median_life:.2f} 小时")
#设备寿命的中位数: 267.41 小时
计算不同分位点的寿命
quantiles = [0.1, 0.25, 0.5, 0.75, 0.9, 0.95]
伽马分布(α=3, β=100)的分位数:
for q in quantiles:
life = gamma_quantile(q, shape, scale)
print(f"P(X ≤ {life:6.2f}) = {q}")
#P(X ≤ 110.21) = 0.1
#P(X ≤ 172.73) = 0.25
#P(X ≤ 267.41) = 0.5
#P(X ≤ 392.04) = 0.75
#P(X ≤ 532.23) = 0.9
#P(X ≤ 629.58) = 0.95
保修期设计：确定90%设备能存活的时间
warranty_period = gamma_quantile(0.9, shape, scale)
print(f"\n建议保修期: {warranty_period:.1f} 小时 (90%设备在此时间内不会故障)")
#建议保修期: 532.2 小时 (90%设备在此时间内不会故障)
示例3. 金融风险管理中的VaR计算
基于伽马分布的风险价值计算
def gamma_var(shape_param, scale_param, confidence_level=0.95):
"""
计算伽马分布的风险价值（Value at Risk）
"""
return gamma_quantile(confidence_level, shape_param, scale_param)
def gamma_expected_shortfall(shape_param, scale_param, confidence_level=0.95):
"""
计算伽马分布的预期短缺（Expected Shortfall）
"""
var_level = gamma_var(shape_param, scale_param, confidence_level)
# 简化计算：使用条件期望公式
conditional_shape = shape_param + 1
conditional_expectation = scale_param * conditional_shape
# 更精确的计算需要积分，这里使用近似
return conditional_expectation * (1 - confidence_level) / (1 - confidence_level)
投资损失分布分析
shape_loss = 2.5 # 损失分布的形状参数
scale_loss = 10000 # 损失分布的尺度参数
var_95 = gamma_var(shape_loss, scale_loss, 0.95)
var_99 = gamma_var(shape_loss, scale_loss, 0.99)
投资损失风险分析:
print(f"95% VaR: ${var_95:,.2f}")
print(f"99% VaR: ${var_99:,.2f}")
#95% VaR: $55,352.49
#99% VaR: $75,431.36
风险资本分配
capital_required = var_99 * 1.1 # 附加10%缓冲
print(f"建议风险资本: ${capital_required:,.2f}")
#建议风险资本: $82,974.50
示例4. 医学统计中的生存分析
生存函数的分位数计算
def survival_quantile(survival_probability, shape_param, scale_param):
"""
计算生存时间的分位数（生存分析）
S(t) = 1 - F(t) = P(T > t)
=> t = gammaincinv(1 - survival_probability, α) * β
"""
return scale_param * gammaincinv(1 - survival_probability, shape_param)
癌症患者生存时间分析
shape_medical = 1.8 # 形状参数
scale_medical = 24 # 尺度参数（月）
计算中位生存时间
median_survival = survival_quantile(0.5, shape_medical, scale_medical)
print(f"中位生存时间: {median_survival:.1f} 个月")
#中位生存时间: 35.5 个月
不同生存概率对应的时间
survival_probs = [0.9, 0.75, 0.5, 0.25, 0.1]
患者生存时间分析:
for prob in survival_probs:
time = survival_quantile(prob, shape_medical, scale_medical)
print(f"{prob * 100:2.0f}%患者生存时间超过 {time:5.1f} 个月")
#90%患者生存时间超过 10.3 个月
#75%患者生存时间超过 19.6 个月
#50%患者生存时间超过 35.5 个月
#25%患者生存时间超过 58.6 个月
#10%患者生存时间超过 86.1 个月
临床试验设计：确定观察期
trial_duration = survival_quantile(0.8, shape_medical, scale_medical)
print(f"\n建议临床试验观察期: {trial_duration:.1f} 个月 (80%患者在此时间内生存)")
#建议临床试验观察期: 16.6 个月 (80%患者在此时间内生存)
示例5. 矩阵输入的批量计算
批量处理多个概率和参数组合
多个概率值
probabilities = [[0.1, 0.5], [0.9, 0.95]]
shape_params = [[2, 3], [4, 5]]
result = gammaincinv(probabilities, shape_params)
逆不完全伽马函数结果矩阵:
print(result)
#[[0.53181161 2.67406031]
[6.68078307 9.15351903]]
保险精算应用：不同风险等级的分位数
risk_levels = [0.5, 0.75, 0.9, 0.95, 0.99]
shape_parameters = [1.5, 2.0, 2.5, 3.0]
不同风险等级和形状参数的分位数:
for shape in shape_parameters:
row = []
for level in risk_levels:
quantile = gamma_quantile(level, shape, 1000)
row.append(f"{quantile:.0f}")
print(f"α={shape}: {', '.join(row)}")
#α=1.5: 1183, 2054, 3126, 3907, 5672
#α=2.0: 1678, 2693, 3890, 4744, 6638
#α=2.5: 2176, 3313, 4618, 5535, 7543
#α=3.0: 2674, 3920, 5322, 6296, 8406
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.special as sc
def gamma_incomplete_inverse_function(input_str):
"""
计算逆不完全伽马函数，类似MATLAB的gammaincinv函数。
参数:
input_str: 输入的字符串表达式，格式为元组(y, a, 'lower'/'upper')或(y, a)。
返回:
标量数值、SymPy矩阵或错误信息字符串。
示例:
>>> gamma_incomplete_inverse_function("(0.5, 1)")
0.6931471806
>>> gamma_incomplete_inverse_function("([[0.5], [0.6]], [[1], [2]])")
Matrix([[1.386294361], [1.42711696]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = sc.gammaincinv(*params[::-1])
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例测试
if __name__ == "__main__":
# 标量测试
print(gamma_incomplete_inverse_function("0.5, 0.6"))
# 0.31570201701610756
Y = gammaln(A)返回gamma函数的对数gammaln(A)=log(gamma(A)).
A 必须是非负数和实数. gammaln命令可避免直接使用log(gamma(A))计算时可能会出现的下溢和上溢.
示例1. 统计学中的最大似然估计
伽玛分布的对数似然函数
def gamma_log_likelihood(data, alpha, beta):
"""
计算伽玛分布的对数似然函数
"""
n = len(data)
sum_log_x = sum(log(data))
sum_x = sum(data)
# 使用gammaln避免数值溢出
log_likelihood = (alpha - 1) * sum_log_x - sum_x / beta - n * alpha * log(beta) - n * gammaln(alpha)
return log_likelihood
示例数据（模拟伽玛分布数据）
np.random.seed(42)
true_alpha, true_beta = 3.0, 2.0
sample_data = np.random.gamma(true_alpha, true_beta, 100)
计算不同参数的对数似然
alpha_values = np.linspace(2.0, 4.0, 50)
beta_values = np.linspace(1.5, 2.5, 50)
log_likelihoods = []
for alpha in alpha_values:
for beta in beta_values:
ll = gamma_log_likelihood(sample_data, alpha, beta)
log_likelihoods.append((alpha, beta, ll))
找到最大似然估计
best_fit = max(log_likelihoods, key=lambda x: x[2])
print(f"真实参数: α={true_alpha}, β={true_beta}")
print(f"最大似然估计: α={best_fit[0]:.3f}, β={best_fit[1]:.3f}")
print(f"最大对数似然值: {best_fit[2]:.3f}")
#真实参数: α=3.0, β=2.0
#最大似然估计: α=3.918, β=1.500
#最大对数似然值: -240.090
示例2. 贝叶斯统计中的模型比较
计算贝叶斯因子（用于模型比较）
def bayes_factor_gamma_models(data, alpha1, beta1, alpha2, beta2):
"""
计算两个伽玛分布模型的贝叶斯因子
"""
# 模型1的对数边际似然
log_marginal1 = gamma_log_likelihood(data, alpha1, beta1) + log(1.0) # 简化先验
# 模型2的对数边际似然
log_marginal2 = gamma_log_likelihood(data, alpha2, beta2) + log(1.0) # 简化先验
bayes_factor = exp(log_marginal1 - log_marginal2)
return bayes_factor
模型比较示例
model1_params = (2.5, 2.0) # (α, β)
model2_params = (3.0, 1.8) # (α, β)
bf = bayes_factor_gamma_models(sample_data, *model1_params, *model2_params)
print(f"贝叶斯因子（模型1/模型2）: {bf:.3f}")
#贝叶斯因子（模型1/模型2）: 0.004
if bf > 1:
print("证据支持模型1")
else:
print("证据支持模型2")
#证据支持模型2
示例3. 组合数学和大数计算
计算大数的组合数对数（避免阶乘溢出）
def log_combination(n, k):
"""
计算组合数C(n,k)的自然对数，避免大数阶乘溢出
"""
if k < 0 or k > n:
return -inf
# 使用Gamma函数性质: n! = Gamma(n+1)
# log(C(n,k)) = log(n!) - log(k!) - log((n-k)!)
return (gammaln(n + 1) -
gammaln(k + 1) -
gammaln(n - k + 1))
大数组合计算示例
n_large = 1000
k_large = 200
直接计算会溢出，但对数计算可行
log_comb = log_combination(n_large, k_large)
print(f"log(C({n_large}, {k_large})) = {log_comb}")
print(f"C({n_large}, {k_large}) ≈ {np.exp(log_comb):.2e}")
#log(C(1000, 200)) = 496.9454605977171
#C(1000, 200) ≈ 6.62e+215
二项分布的概率质量函数（对数形式）
def log_binomial_pmf(k, n, p):
"""
二项分布概率质量函数的对数
"""
if k < 0 or k > n or p < 0 or p > 1:
return -inf
log_comb = log_combination(n, k)
log_prob = log_comb + k * log(p) + (n - k) * log(1 - p)
return log_prob
计算罕见事件的概率
n_rare = 10000
p_rare = 0.001
k_rare = 15
log_p_rare = log_binomial_pmf(k_rare, n_rare, p_rare)
print(f"罕见事件概率: P(X={k_rare}) = {exp(log_p_rare):.2e}")
#罕见事件概率: P(X=15) = 3.47e-02
示例4. 物理化学中的统计力学应用
理想气体的熵计算（使用伽玛函数）
def ideal_gas_entropy(N, V, T, mass):
"""
计算单原子理想气体的萨库尔-特罗德方程（Sackur-Tetrode）熵
"""
import scipy.constants as const
# 热波长
lambda_thermal = const.h / sqrt(2 * pi * mass * const.k * T)
# 使用斯特林近似和伽玛函数
# S = N*k*[ln(V/Nλ³) + 5/2]
# 更精确的版本涉及伽玛函数
# 使用gammaln计算精确的构型数对数
log_omega = (N * log(V) - gammaln(N + 1) +
3 * N / 2 * log(2 * pi * mass * const.k * T / const.h ** 2) -
3 * N / 2)
entropy = const.k * log_omega
return entropy
# 计算氦气在标准条件下的熵
N_avogadro = 6.022e23 # 阿伏伽德罗常数
V_molar = 0.0224 # 标准状况下气体摩尔体积 (m³)
T_std = 273.15 # 标准温度 (K)
mass_he = 6.646e-27 # 氦原子质量 (kg)
S_molar = ideal_gas_entropy(N_avogadro, V_molar, T_std, mass_he)
print(f"氦气在标准状况下的摩尔熵: {S_molar:.2f} J/(mol·K)")
#氦气在标准状况下的摩尔熵: 99.27 J/(mol·K)
示例5. 金融工程中的极值理论
广义极值分布（GEV）的对数似然
def gev_log_likelihood(data, mu, sigma, xi):
"""
广义极值分布的对数似然函数
"""
n = len(data)
z = (data - mu) / sigma
if abs(xi) < 1e-10: # Gumbel case (xi = 0)
log_likelihood = -n * log(sigma) - sum(z) - sum(exp(-z))
else: # Fréchet or Weibull case (xi ≠ 0)
# 需要确保1 + xi*z > 0
valid = 1 + xi * z > 0
if not all(valid):
return -inf
log_likelihood = (-n * log(sigma) -
(1 + 1 / xi) * sum(log(1 + xi * z)) -
sum((1 + xi * z) ** (-1 / xi)))
return log_likelihood
极值风险分析示例（模拟极端损失数据）
np.random.seed(123)
extreme_losses = np.random.gumbel(loc=1000, scale=200, size=50) # 极端损失
使用gammaln辅助模型选择
def model_complexity_penalty(k, n):
"""
计算模型复杂度的惩罚项（使用伽玛函数）
"""
# BIC近似: 0.5 * k * log(n)
# 更精确的版本涉及伽玛函数
return 0.5 * (k * log(n) - 2 * gammaln(k / 2 + 1))
比较不同复杂度模型
models = [
{"name": "简单模型", "params": 2, "log_likelihood": -1500},
{"name": "复杂模型", "params": 5, "log_likelihood": -1490}
]
n_data = len(extreme_losses)
for model in models:
penalty = model_complexity_penalty(model["params"], n_data)
bic = -2 * model["log_likelihood"] + penalty
print(f"{model['name']}: BIC = {bic:.1f}")
#简单模型: BIC = 3003.9
#复杂模型: BIC = 2988.6
示例6. 生物信息学中的序列分析
DNA序列的熵计算（使用伽玛函数进行狄利克雷先验的贝叶斯估计）
def dna_sequence_entropy(sequence, pseudocount=1):
"""
计算DNA序列的熵，使用狄利克雷先验进行贝叶斯平滑
"""
from collections import Counter
# 计数各碱基出现次数
counts = Counter(sequence)
bases = ['A', 'C', 'G', 'T']
# 添加伪计数（狄利克雷先验）
total_count = len(sequence)
alpha = pseudocount # 对称狄利克雷参数
# 贝叶斯估计的概率
probs = {}
for base in bases:
count_base = counts.get(base, 0)
probs[base] = (count_base + alpha) / (total_count + 4 * alpha)
# 计算熵
entropy = -sum(p * log(p) for p in probs.values() if p > 0)
# 使用gammaln计算模型证据（用于模型比较）
log_evidence = (gammaln(4 * alpha) -
4 * gammaln(alpha) +
sum(gammaln(counts.get(base, 0) + alpha) for base in bases) -
gammaln(total_count + 4 * alpha))
return entropy, log_evidence
DNA序列分析示例
dna_sequence = "ATGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGCTAGC"
entropy, log_evidence = dna_sequence_entropy(dna_sequence)
print(f"DNA序列熵: {entropy:.3f}")
print(f"模型对数证据: {log_evidence:.3f}")
#DNA序列熵: 1.386
#模型对数证据: -59.324
比较不同序列的复杂性
sequences = [
"AAAAAATTTTTTGGGGGGCCCCCC", # 低复杂性
"ATCGATCGATCGATCGATCGATCG", # 中等复杂性
"ATGCTAGCTAGCTAGCTAGCTAGC" # 高复杂性
]
for i, seq in enumerate(sequences):
entropy, _ = dna_sequence_entropy(seq)
print(f"序列{i + 1}熵: {entropy:.3f}")
#序列1熵: 1.386
#序列2熵: 1.386
#序列3熵: 1.386
示例7. 机器学习中的变分推断
变分贝叶斯中的ELBO（证据下界）计算
def variational_elbo_gamma(data, alpha_q, beta_q, alpha_prior, beta_prior):
"""
计算伽玛分布变分近似的ELBO
"""
n = len(data)
# 期望对数似然
E_log_lambda = log(beta_q) + digamma(alpha_q) # digamma是Gamma函数的对数导数
E_lambda = alpha_q / beta_q
expected_log_likelihood = (n * E_log_lambda - E_lambda * sum(data) +
(E_log_lambda - 1) * sum(log(data)) -
n * gammaln(1)) # 简化
# KL散度（变分分布与先验之间）
kl_divergence = ((alpha_q - alpha_prior) * digamma(alpha_q) -
gammaln(alpha_q) + gammaln(alpha_prior) +
alpha_prior * (log(beta_q) - log(beta_prior)) +
alpha_q * (beta_prior - beta_q) / beta_q)
elbo = expected_log_likelihood - kl_divergence
return elbo
变分推断示例
variational_params = (2.5, 1.5) # (alpha_q, beta_q)
prior_params = (1.0, 1.0) # (alpha_prior, beta_prior)
elbo = variational_elbo_gamma(sample_data, *variational_params, *prior_params)
print(f"变分下界(ELBO): {elbo:.3f}")
#变分下界(ELBO): -845.937
示例8. 符号计算的高级应用
符号微分和级数展开
计算Gamma函数对数的导数（Digamma函数）
digamma_expr = diff(loggamma(z), z)
print(f"d(logΓ(z))/dz = {digamma_expr}")
#d(logΓ(z))/dz = polygamma(0, z)
泰勒级数展开
taylor_series = series(loggamma(z), z, 1, 5) # 在z=1处展开到5阶
print(f"logΓ(z)在z=1处的泰勒展开:\n{taylor_series}")
#logΓ(z)在z=1处的泰勒展开:
#-EulerGamma*(z - 1) + pi**2*(z - 1)**2/12 - (z - 1)**3*zeta(3)/3 + pi**4*(z - 1)**4/360 + O((z - 1)**5, (z, 1))
特殊值计算
special_values = [
(1, gammaln(1)),
(0.5, gammaln(0.5)),
(2, gammaln(2)),
]
for z_val, log_gamma_val in special_values:
print(f"logΓ({z_val}) = {log_gamma_val}")
#logΓ(1) = 0.0
#logΓ(0.5) = 0.5723649429247
#logΓ(2) = 0.0
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from scipy.special import gammaln
def gamma_logarithm(input_str):
"""
计算gamma函数的自然对数，对标MATLAB的gammaln函数。
参数:
input_str (str): 输入的数学表达式字符串，支持标量、矩阵或符号表达式。
对于数值输入，必须为非负实数；符号表达式需为合法SymPy表达式。
返回:
SymPy表达式/float/np.ndarray/str:
- 符号输入返回SymPy的loggamma表达式
- 标量数值输入返回float结果
- 矩阵输入返回NumPy数组
- 错误时返回字符串描述
示例:
>>> gamma_logarithm("5") # 标量数值输入
3.1780538303479458
>>> gamma_logarithm("z") # 符号输入
loggamma(z)
>>> gamma_logarithm("[[1, 2], [3, 4]]") # 矩阵输入
array([[0. , 0.69314718],
[1.79175947, 3.17805383]])
>>> gamma_logarithm("-1") # 错误输入
'错误：输入必须为非负实数。'
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if expr.free_symbols:
result = sp.loggamma(expr)
elif expr.is_number:
# 数值型输入验证
z = float(expr)
# 使用SciPy计算数值结果
result = gammaln(z)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except (sp.SympifyError, TypeError) as e:
return f"语法错误：无法解析输入 '{input_str}' -> {e}"
except ValueError as e:
return f"计算错误：{e}"
except Exception as e:
return f"未知错误：{type(e).__name__} - {str(e)}"
# 符号表达式
print(gamma_logarithm("x+y"))
# loggamma(x + y)
# 标量计算测试
print(gamma_logarithm("6"))
# 4.787491742782046
伽玛概率密度函数
p＝gampdf（x，a）返回具有a中的形状参数的标准gamma分布的概率密度函数（pdf），在x中的值处进行评估。.
y=gampdf（x，a，b）返回具有形状参数a和比例参数b的伽马分布的pdf,在x中的值处进行评估.
x是评估pdf的值,非负标量,非负标量数组
a是形状参数,正标量值,正标量值数组
b是比例参数,默认值是1,正标量值,正标量的数组.
如果输入参数x,a和b中的一个或多个是数组,则数组大小必须相同.在这种情况下,gampdf将每个标量输入扩展为与数组输入大小相同的常量数组.p中的每个元素都是由a和b中的对应元素指定的分布的cdf值, 在x中的相应元素处进行评估.
示例1. 可靠性工程中的设备寿命建模
import matplotlib.pyplot as plt
设备寿命的伽马分布建模
def plot_gamma_lifetimes():
"""绘制不同形状参数的伽马分布PDF，模拟设备寿命"""
x = np.linspace(0, 20, 1000)
# 不同形状参数（故障模式）
shape_params = [1, 2, 3, 5]
scale_param = 2 # 平均寿命参数
plt.figure(figsize=(10, 6))
for k in shape_params:
# 计算PDF值
pdf_values = [gampdf(xi, k, scale_param) for xi in x]
plt.plot(x, pdf_values, label=f'形状参数 k={k}')
plt.title('伽马分布模拟设备寿命概率密度')
plt.xlabel('时间 (千小时)')
plt.ylabel('概率密度')
plt.legend()
plt.grid(True)
plt.show()
计算特定设备的故障概率密度
def equipment_failure_density(operating_time, shape_param, scale_param):
"""计算设备在特定运行时间下的故障概率密度"""
return gampdf(operating_time, shape_param, scale_param)
示例：计算设备运行5000小时后的故障概率密度
operating_hours = 5 # 千小时
shape = 3.0 # 形状参数（与故障机制相关）
scale = 2.0 # 尺度参数（平均寿命）
density = equipment_failure_density(operating_hours, shape, scale)
print(f"设备运行{operating_hours}千小时后的故障概率密度: {density:.4f}")
#设备运行5千小时后的故障概率密度: 0.1283
plot_gamma_lifetimes()
示例2. 保险精算中的索赔金额分布
保险索赔金额的伽马分布建模
def insurance_claim_model():
"""模拟保险索赔金额的分布"""
claim_amounts = np.linspace(0, 50000, 1000) # 索赔金额范围
# 不同风险等级的索赔分布
risk_profiles = [
{"name": "低风险", "shape": 2.0, "scale": 5000},
{"name": "中风险", "shape": 1.5, "scale": 8000},
{"name": "高风险", "shape": 1.2, "scale": 12000}
]
plt.figure(figsize=(12, 8))
for profile in risk_profiles:
pdf_values = []
for amount in claim_amounts:
pdf_val = gampdf(amount, profile['shape'], profile['scale'])
pdf_values.append(pdf_val)
plt.plot(claim_amounts, pdf_values,
label=f"{profile['name']} (k={profile['shape']}, θ={profile['scale']})")
plt.title('不同风险等级的保险索赔金额概率密度')
plt.xlabel('索赔金额 ($)')
plt.ylabel('概率密度')
plt.legend()
plt.grid(True)
plt.xlim(0, 50000)
plt.show()
计算特定索赔金额的概率密度
def claim_probability_density(claim_amount, risk_level="medium"):
"""计算特定索赔金额的概率密度"""
risk_params = {
"low": (2.0, 5000),
"medium": (1.5, 8000),
"high": (1.2, 12000)
}
shape, scale = risk_params.get(risk_level, (1.5, 8000))
return gampdf(f"({claim_amount}, {shape}, {scale})")
示例计算
claim_10k = claim_probability_density(10000, "medium")
claim_30k = claim_probability_density(30000, "high")
print(f"中等风险索赔$10,000的概率密度: {claim_10k:.6f}")
print(f"高风险索赔$30,000的概率密度: {claim_30k:.6f}")
#中等风险索赔$10,000的概率密度: 0.000045
#高风险索赔$30,000的概率密度: 0.000009
insurance_claim_model()
示例3. 医学统计中的治疗时间分析
医疗治疗时间的伽马分布分析
def medical_treatment_analysis():
"""分析不同医疗程序的持续时间分布"""
treatment_times = np.linspace(0, 120, 1000) # 治疗时间（分钟）
# 不同医疗程序的时间分布
procedures = [
{"name": "常规检查", "shape": 3.0, "scale": 15},
{"name": "小型手术", "shape": 2.5, "scale": 25},
{"name": "复杂手术", "shape": 2.0, "scale": 40}
]
plt.figure(figsize=(12, 8))
for procedure in procedures:
pdf_values = []
for time in treatment_times:
pdf_val = gampdf(time, procedure['shape'], procedure['scale'])
pdf_values.append(pdf_val)
plt.plot(treatment_times, pdf_values,
label=f"{procedure['name']} (k={procedure['shape']}, θ={procedure['scale']})")
plt.title('不同医疗程序持续时间的概率密度分布')
plt.xlabel('治疗时间 (分钟)')
plt.ylabel('概率密度')
plt.legend()
plt.grid(True)
plt.show()
手术室调度优化
def operating_room_scheduling():
"""使用伽马分布优化手术室调度"""
# 计算不同手术时间的概率
surgery_types = [
("阑尾切除术", 2.2, 45),
("胆囊切除术", 2.5, 60),
("心脏搭桥", 1.8, 180)
]
print("手术时间概率分析:")
for surgery_name, shape, scale in surgery_types:
# 计算90分钟内的完成概率（使用CDF，这里演示PDF用途）
time_points = [30, 60, 90, 120]
densities = []
for t in time_points:
density_val = gampdf(f"({t}, {shape}, {scale})")
densities.append(density_val)
print(f"{surgery_name}在{t}分钟时的概率密度: {density_val:.6f}")
print("---")
return surgery_types
medical_treatment_analysis()
surgery_data = operating_room_scheduling()
#手术时间概率分析:
#阑尾切除术在30分钟时的概率密度: 0.006366
#阑尾切除术在60分钟时的概率密度: 0.007508
#阑尾切除术在90分钟时的概率密度: 0.006271
#阑尾切除术在120分钟时的概率密度: 0.004547
#---
#胆囊切除术在30分钟时的概率密度: 0.002689
#胆囊切除术在60分钟时的概率密度: 0.004612
#胆囊切除术在90分钟时的概率密度: 0.005139
#胆囊切除术在120分钟时的概率密度: 0.004799
#---
#心脏搭桥在30分钟时的概率密度: 0.001204
#心脏搭桥在60分钟时的概率密度: 0.001775
#心脏搭桥在90分钟时的概率密度: 0.002078
#心脏搭桥在120分钟时的概率密度: 0.002214
示例4. 金融中的股价波动建模
使用伽马分布建模金融资产回报率
def stock_volatility_model():
"""建模股价波动率的分布"""
volatilities = np.linspace(0, 0.5, 1000) # 波动率范围 (0-50%)
# 不同资产类别的波动率分布
asset_classes = [
{"name": "国债", "shape": 4.0, "scale": 0.02},
{"name": "蓝筹股", "shape": 2.5, "scale": 0.05},
{"name": "科技股", "shape": 1.8, "scale": 0.08},
{"name": "加密货币", "shape": 1.2, "scale": 0.15}
]
plt.figure(figsize=(12, 8))
for asset in asset_classes:
pdf_values = []
for vol in volatilities:
pdf_val = gampdf(vol, asset['shape'], asset['scale'])
pdf_values.append(pdf_val)
plt.plot(volatilities, pdf_values,
label=f"{asset['name']} (k={asset['shape']}, θ={asset['scale']:.3f})")
plt.title('不同资产类别的波动率概率密度分布')
plt.xlabel('年化波动率')
plt.ylabel('概率密度')
plt.legend()
plt.grid(True)
plt.show()
风险价值(VaR)计算辅助
def volatility_risk_assessment():
"""评估特定波动率水平的风险"""
volatility_levels = [0.1, 0.2, 0.3, 0.4] # 10%-40%波动率
print("科技股波动率风险评估:")
tech_shape, tech_scale = 1.8, 0.08
for vol in volatility_levels:
density_val = gampdf(vol, tech_shape, tech_scale)
risk_level = "低" if density_val > 2 else "中" if density_val > 0.5 else "高"
print(f"波动率{vol * 100:.0f}%: 概率密度={density_val:.4f} ({risk_level}风险)")
return tech_shape, tech_scale
stock_volatility_model()
tech_params = volatility_risk_assessment()
#科技股波动率风险评估:
#波动率10%: 概率密度=4.5966 (低风险)
#波动率20%: 概率密度=2.2930 (低风险)
#波动率30%: 概率密度=0.9087 (中风险)
#波动率40%: 概率密度=0.3277 (高风险)
示例5. 环境科学中的降雨量分析
降雨量的伽马分布建模
def rainfall_distribution_model():
"""建模月降雨量的概率分布"""
rainfall_mm = np.linspace(0, 300, 1000) # 月降雨量 (mm)
# 不同气候区域的降雨分布
climate_zones = [
{"name": "干旱区", "shape": 1.5, "scale": 20},
{"name": "温带区", "shape": 2.0, "scale": 50},
{"name": "湿润区", "shape": 2.5, "scale": 80}
]
plt.figure(figsize=(12, 8))
for zone in climate_zones:
pdf_values = []
for rain in rainfall_mm:
pdf_val = gampdf(rain, zone['shape'], zone['scale'])
pdf_values.append(pdf_val)
plt.plot(rainfall_mm, pdf_values,
label=f"{zone['name']} (k={zone['shape']}, θ={zone['scale']})")
plt.title('不同气候区域月降雨量的概率密度分布')
plt.xlabel('月降雨量 (mm)')
plt.ylabel('概率密度')
plt.legend()
plt.grid(True)
plt.show()
干旱风险评估
def drought_risk_assessment():
"""评估不同降雨量水平的干旱风险"""
rainfall_thresholds = [10, 30, 50, 100] # 降雨量阈值 (mm)
# 干旱区的参数
arid_shape, arid_scale = 1.5, 20
print("干旱区降雨量风险评估:")
for threshold in rainfall_thresholds:
density_val = gampdf(threshold, arid_shape, arid_scale)
if threshold <= 10:
risk = "极端干旱"
elif threshold <= 30:
risk = "严重干旱"
elif threshold <= 50:
risk = "中等干旱"
else:
risk = "正常"
print(f"降雨量{threshold}mm: 概率密度={density_val:.4f} ({risk})")
return arid_shape, arid_scale
rainfall_distribution_model()
arid_params = drought_risk_assessment()
#干旱区降雨量风险评估:
#降雨量10mm: 概率密度=0.0242 (极端干旱)
#降雨量30mm: 概率密度=0.0154 (严重干旱)
#降雨量50mm: 概率密度=0.0073 (中等干旱)
#降雨量100mm: 概率密度=0.0009 (正常)
示例6. 工业生产中的质量控制
产品尺寸偏差的伽马分布建模
def manufacturing_quality_control():
"""建模产品制造过程中的尺寸偏差分布"""
deviations = np.linspace(0, 0.5, 1000) # 尺寸偏差 (mm)
# 不同制造工艺的偏差分布
processes = [
{"name": "传统加工", "shape": 1.8, "scale": 0.1},
{"name": "数控加工", "shape": 2.5, "scale": 0.05},
{"name": "精密加工", "shape": 3.5, "scale": 0.02}
]
plt.figure(figsize=(12, 8))
for process in processes:
pdf_values = []
for dev in deviations:
pdf_val = gampdf(dev, process['shape'], process['scale'])
pdf_values.append(pdf_val)
plt.plot(deviations, pdf_values,
label=f"{process['name']} (k={process['shape']}, θ={process['scale']:.3f})")
plt.title('不同制造工艺的尺寸偏差概率密度分布')
plt.xlabel('尺寸偏差 (mm)')
plt.ylabel('概率密度')
plt.legend()
plt.grid(True)
plt.xlim(0, 0.3)
plt.show()
公差设计优化
def tolerance_design_analysis():
"""分析不同公差要求下的质量水平"""
tolerance_limits = [0.05, 0.1, 0.15, 0.2] # 公差限 (mm)
# 数控加工的参数
cnc_shape, cnc_scale = 2.5, 0.05
print("数控加工公差设计分析:")
for tolerance in tolerance_limits:
density_val = gampdf(tolerance, cnc_shape, cnc_scale)
# 估算合格率（简化计算）
if density_val > 1:
yield_estimate = ">99%"
elif density_val > 0.1:
yield_estimate = "95%-99%"
else:
yield_estimate = "<95%"
print(f"公差±{tolerance}mm: 概率密度={density_val:.4f} (预估合格率: {yield_estimate})")
return cnc_shape, cnc_scale
manufacturing_quality_control()
cnc_params = tolerance_design_analysis()
#数控加工公差设计分析:
#公差±0.05mm: 概率密度=5.5348 (预估合格率: >99%)
#公差±0.1mm: 概率密度=5.7590 (预估合格率: >99%)
#公差±0.15mm: 概率密度=3.8922 (预估合格率: >99%)
#公差±0.2mm: 概率密度=2.2045 (预估合格率: >99%)
示例7. 矩阵输入的批量计算示例
批量处理多个参数组合
def batch_gamma_pdf_calculation():
"""演示矩阵输入情况下的批量计算"""
# 创建测试数据矩阵
x_values = [[1, 2, 3], [4, 5, 6]] # 输入值矩阵
k_values = [[1, 1, 1], [2, 2, 2]] # 形状参数矩阵
theta_values = [[1, 2, 1], [1, 1, 2]] # 尺度参数矩阵
print("矩阵输入批量计算示例:")
print("x值矩阵:")
print(np.array(x_values))
print("\n形状参数k矩阵:")
print(np.array(k_values))
print("\n尺度参数θ矩阵:")
print(np.array(theta_values))
# 批量计算PDF值
result = gampdf(x_values, k_values, theta_values)
print("\n伽马分布PDF结果矩阵:")
print(result)
return result
实际应用：批量风险评估
def batch_risk_assessment():
"""批量评估多个风险场景"""
# 不同场景的参数
scenarios = {
"设备A_低负荷": (1000, 2.0, 500),
"设备A_高负荷": (1000, 1.5, 300),
"设备B_正常": (2000, 2.5, 800),
"设备B_老化": (2000, 1.8, 600)
}
print("\n设备故障风险批量评估:")
for scenario, params in scenarios.items():
operating_time, shape, scale = params
density_val = gampdf(operating_time, shape, scale)
risk_level = "低" if density_val < 0.001 else "中" if density_val < 0.005 else "高"
print(f"{scenario}: 概率密度={density_val:.6f} ({risk_level}风险)")
return scenarios
执行批量计算示例
matrix_results = batch_gamma_pdf_calculation()
#矩阵输入批量计算示例:
#x值矩阵:
#[[1 2 3]
[4 5 6]]
#形状参数k矩阵:
#[[1 1 1]
[2 2 2]]
#尺度参数θ矩阵:
#[[1 2 1]
[1 1 2]]
#伽马分布PDF结果矩阵:
#[[0.36787944 0.18393972 0.04978707]
[0.07326256 0.03368973 0.0746806 ]]
risk_scenarios = batch_risk_assessment()
#设备故障风险批量评估:
#设备A_低负荷: 概率密度=0.000541 (低风险)
#设备A_高负荷: 概率密度=0.000245 (低风险)
#设备B_正常: 概率密度=0.000305 (低风险)
#设备B_老化: 概率密度=0.000167 (低风险)
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Geometric, density, cdf, Gamma, E, variance
from scipy.stats import gamma as scipy_gamma
def gamma_pdf(input_str):
"""
计算伽马分布概率密度函数 (对标MATLAB gampdf)
参数:
x: 输入值/矩阵
k: 形状参数/矩阵
theta: 尺度参数/矩阵 (默认1)
返回:
SymPy矩阵/数值/错误信息
示例:
>>> gamma_pdf(2, 1) # 标量计算
0.1353352832366127
>>> gamma_pdf([[1,2],[3,4]], 2, 1) # 矩阵输入
Matrix([
[0.367879441171442, 0.270670566473225],
[0.149361205103592, 0.073262555554937]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) <= 3:
if len(expr) == 2:
theta = 1
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = scipy_gamma.pdf(*params, theta)
elif any(e.free_symbols for e in expr):
G = Gamma("gamma", expr[1], theta)
result = density(G)(expr[0])
elif len(expr) == 3:
if all(e.is_number for e in expr):
result = scipy_gamma.pdf(*expr)
elif any(e.free_symbols for e in expr):
G = Gamma("gamma", expr[1], expr[2])
result = density(G)(expr[0])
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 标量测试
print(gamma_pdf("2, 1"))
# 0.36787944117144233
# 符号变量测试
print(gamma_pdf("x, 1"))
# exp(-x)
伽玛随机数
r = gamrnd(a,b) 从具有形状参数 a 和尺度参数 b 的 gamma 分布中生成一个随机数.
r = gamrnd(a,b,sz) 从 gamma 分布中生成一个随机数数组，其中向量 sz 指定 size(r).
a是形状参数,正标量值,正标量值数组
b是比例参数,默认值是1,正标量值,正标量的数组.
如果输入参数a和b都是数组,则数组大小必须相同.如果a或b是标量,则gamrnd会将标量参量扩展为与另一个参量大小相同的常量数组.r中的每个元素均是从a和b中对应元素所指定的分布中生成的随机数.
示例1. 可靠性工程中的设备寿命模拟
import matplotlib.pyplot as plt
import seaborn as sns
np.random.seed(42) # 设置随机种子保证可重复性
不同设备的伽马分布参数
equipment_types = {
"服务器": {"shape": 2.5, "scale": 1000}, # 平均寿命2500小时
"网络设备": {"shape": 3.0, "scale": 800}, # 平均寿命2400小时
"存储设备": {"shape": 4.0, "scale": 600}, # 平均寿命2400小时
"工作站": {"shape": 2.0, "scale": 1200} # 平均寿命2400小时
}
为每种设备生成1000个寿命样本
samples = {}
for eq_type, params in equipment_types.items():
shape, scale = params["shape"], params["scale"]
samples[eq_type] = gamrnd(shape, scale, 1000)
绘制寿命分布直方图
plt.figure(figsize=(12, 8))
for i, (eq_type, data) in enumerate(samples.items()):
plt.subplot(2, 2, i + 1)
plt.hist(data, bins=50, alpha=0.7, density=True)
plt.title(f'{eq_type}寿命分布')
plt.xlabel('寿命 (小时)')
plt.ylabel('密度')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
设备寿命模拟结果:
for eq_type, data in samples.items():
mean_life = np.mean(data)
std_life = np.std(data)
failure_before_2000 = np.mean(data < 2000) * 100
print(f"{eq_type}: 平均寿命={mean_life:.0f}小时, "
f"标准差={std_life:.0f}小时, "
f"2000小时内故障概率={failure_before_2000:.1f}%")
#服务器: 平均寿命=2569小时, 标准差=1561小时, 2000小时内故障概率=42.6%
#网络设备: 平均寿命=2420小时, 标准差=1402小时, 2000小时内故障概率=44.1%
#存储设备: 平均寿命=2369小时, 标准差=1176小时, 2000小时内故障概率=43.2%
#工作站: 平均寿命=2406小时, 标准差=1727小时, 2000小时内故障概率=50.6%
示例2. 保险精算中的索赔金额模拟
np.random.seed(123)
不同风险等级的索赔分布参数
risk_categories = {
"低风险": {"shape": 3.0, "scale": 5000}, # 平均索赔15000
"中风险": {"shape": 2.0, "scale": 10000}, # 平均索赔20000
"高风险": {"shape": 1.5, "scale": 15000} # 平均索赔22500
}
生成索赔金额样本
n_claims = 500 # 每种风险500个索赔
claim_samples = {}
for category, params in risk_categories.items():
shape, scale = params["shape"], params["scale"]
claim_samples[category] = gamrnd(shape, scale, n_claims)
可视化索赔分布
plt.figure(figsize=(14, 6))
箱线图比较
plt.subplot(1, 2, 1)
data_to_plot = [claim_samples[cat] for cat in risk_categories.keys()]
plt.boxplot(data_to_plot, labels=risk_categories.keys())
plt.title('不同风险等级的索赔金额分布')
plt.ylabel('索赔金额 ($)')
plt.grid(True, alpha=0.3)
密度图
plt.subplot(1, 2, 2)
for category, data in claim_samples.items():
sns.kdeplot(data, label=category, alpha=0.7)
plt.title('索赔金额概率密度')
plt.xlabel('索赔金额 ($)')
plt.ylabel('密度')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
保险索赔风险评估:
for category, data in claim_samples.items():
avg_claim = np.mean(data)
p95_claim = np.percentile(data, 95) # 95%分位数
extreme_claims = np.sum(data > 100000) # 极端索赔数量
print(f"{category}: 平均索赔=${avg_claim:,.0f}, "
f"95%分位数=${p95_claim:,.0f}, "
f"极端索赔(>$100k)数量={extreme_claims}")
#低风险: 平均索赔=$15,222, 95%分位数=$31,256, 极端索赔(>$100k)数量=0
#中风险: 平均索赔=$20,324, 95%分位数=$46,446, 极端索赔(>$100k)数量=0
#高风险: 平均索赔=$22,899, 95%分位数=$58,527, 极端索赔(>$100k)数量=1
示例3. 医疗统计中的治疗时间模拟
np.random.seed(456)
medical_procedures = {
"常规检查": {"shape": 3.0, "scale": 10}, # 平均30分钟
"小型手术": {"shape": 2.5, "scale": 20}, # 平均50分钟
"专科诊疗": {"shape": 2.0, "scale": 25}, # 平均50分钟
"急诊处理": {"shape": 1.8, "scale": 15} # 平均27分钟
}
生成治疗时间样本
n_cases = 200
treatment_times = {}
for procedure, params in medical_procedures.items():
shape, scale = params["shape"], params["scale"]
treatment_times[procedure] = gamrnd(shape, scale, n_cases)
医疗治疗时间模拟分析:
plt.figure(figsize=(12, 8))
for i, (procedure, times) in enumerate(treatment_times.items()):
# 计算统计量
mean_time = np.mean(times)
median_time = np.median(times)
p90_time = np.percentile(times, 90)
print(f"{procedure}: 平均={mean_time:.1f}分钟, "
f"中位数={median_time:.1f}分钟, "
f"90%分位数={p90_time:.1f}分钟")
# 绘制累积分布函数
sorted_times = np.sort(times)
cdf = np.arange(1, len(sorted_times) + 1) / len(sorted_times)
plt.plot(sorted_times, cdf, label=procedure, linewidth=2)
#常规检查: 平均=29.9分钟, 中位数=25.5分钟, 90%分位数=52.9分钟
#小型手术: 平均=47.8分钟, 中位数=41.1分钟, 90%分位数=90.8分钟
#专科诊疗: 平均=50.4分钟, 中位数=42.6分钟, 90%分位数=94.1分钟
#急诊处理: 平均=26.2分钟, 中位数=20.2分钟, 90%分位数=51.8分钟
plt.title('不同医疗程序治疗时间的累积分布')
plt.xlabel('治疗时间 (分钟)')
plt.ylabel('累积概率')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
手术室利用率分析
operating_room_hours = 8 * 60 # 8小时工作制
procedure_sequence = ["常规检查", "小型手术", "专科诊疗", "急诊处理"]
手术室调度模拟:
current_time = 0
for procedure in procedure_sequence:
# 随机选择一个该程序的治疗时间
proc_time = np.random.choice(treatment_times[procedure])
current_time += proc_time
if current_time <= operating_room_hours:
status = "完成"
else:
status = "超时"
print(f"{procedure}: {proc_time:.1f}分钟, 累计{current_time:.1f}分钟 - {status}")
#常规检查: 44.7分钟, 累计44.7分钟 - 完成
#小型手术: 28.1分钟, 累计72.8分钟 - 完成
#专科诊疗: 34.6分钟, 累计107.4分钟 - 完成
#急诊处理: 6.9分钟, 累计114.3分钟 - 完成
示例4. 金融风险管理中的损失模拟
np.random.seed(789)
不同资产类别的损失分布参数
asset_classes = {
"债券投资": {"shape": 4.0, "scale": 10000}, # 低风险
"股票投资": {"shape": 2.5, "scale": 50000}, # 中等风险
"衍生品投资": {"shape": 1.5, "scale": 100000}, # 高风险
"风险投资": {"shape": 1.2, "scale": 200000} # 极高风险
}
生成损失样本
n_scenarios = 10000 # 10000个模拟场景
loss_samples = {}
for asset_class, params in asset_classes.items():
shape, scale = params["shape"], params["scale"]
loss_samples[asset_class] = gamrnd(shape, scale, n_scenarios)
风险价值(VaR)和预期短缺(ES)计算
confidence_levels = [0.95, 0.99]
金融风险度量分析:
plt.figure(figsize=(15, 10))
for i, (asset_class, losses) in enumerate(loss_samples.items()):
plt.subplot(2, 2, i + 1)
# 计算风险指标
var_values = {}
es_values = {}
for cl in confidence_levels:
var = np.percentile(losses, cl * 100)
# ES是超过VaR的条件期望
excess_losses = losses[losses > var]
es = np.mean(excess_losses) if len(excess_losses) > 0 else var
var_values[cl] = var
es_values[cl] = es
# 绘制损失分布和VaR线
plt.hist(losses, bins=100, alpha=0.7, density=True)
for cl in confidence_levels:
plt.axvline(var_values[cl], color='red', linestyle='--',
label=f'{cl * 100}% VaR = ${var_values[cl]:,.0f}')
plt.title(f'{asset_class}损失分布')
plt.xlabel('损失金额 ($)')
plt.ylabel('密度')
plt.legend()
plt.grid(True, alpha=0.3)
# 输出风险指标
print(f"\n{asset_class}:")
for cl in confidence_levels:
print(f" {cl * 100}% VaR: ${var_values[cl]:,.0f}")
print(f" {cl * 100}% ES: ${es_values[cl]:,.0f}")
#债券投资:
# 95.0% VaR: $78,271
# 95.0% ES: $93,316
# 99.0% VaR: $103,080
# 99.0% ES: $117,009
#股票投资:
# 95.0% VaR: $278,903
# 95.0% ES: $338,190
# 99.0% VaR: $368,382
# 99.0% ES: $429,035
#衍生品投资:
# 95.0% VaR: $393,148
# 95.0% ES: $504,418
# 99.0% VaR: $569,512
# 99.0% ES: $683,083
plt.tight_layout()
plt.show()
示例5. 环境科学中的降雨量模拟
np.random.seed(321)
climate_regions = {
"干旱地区": {"shape": 1.8, "scale": 15}, # 低降雨量
"半干旱地区": {"shape": 2.2, "scale": 30}, # 中等降雨量
"湿润地区": {"shape": 2.5, "scale": 60}, # 高降雨量
"热带雨林": {"shape": 3.0, "scale": 100} # 极高降雨量
}
生成120个月（10年）的降雨量数据
n_months = 120
rainfall_data = {}
for region, params in climate_regions.items():
shape, scale = params["shape"], params["scale"]
rainfall_data[region] = gamrnd(shape, scale, n_months)
降雨量模式分析:
plt.figure(figsize=(15, 10))
时间序列图
plt.subplot(2, 1, 1)
for region, rainfall in rainfall_data.items():
plt.plot(rainfall, label=region, alpha=0.7)
plt.axhline(y=30, color='red', linestyle='--', label='干旱阈值 (30mm)')
plt.title('月降雨量时间序列 (10年)')
plt.xlabel('月份')
plt.ylabel('降雨量 (mm)')
plt.legend()
plt.grid(True, alpha=0.3)
季节性分析（按月份分组）
plt.subplot(2, 1, 2)
months = np.arange(12)
for region, rainfall in rainfall_data.items():
seasonal_means = []
for month in months:
# 获取所有该月份的降雨量
month_data = rainfall[month::12]
seasonal_means.append(np.mean(month_data))
plt.plot(months, seasonal_means, 'o-', label=region, markersize=8)
plt.title('月平均降雨量季节性模式')
plt.xlabel('月份')
plt.ylabel('平均降雨量 (mm)')
plt.xticks(months, ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun',
'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
干旱分析
drought_threshold = 30 # 月降雨量低于30mm视为干旱
干旱频率分析:
for region, rainfall in rainfall_data.items():
drought_months = np.sum(rainfall < drought_threshold)
drought_frequency = drought_months / n_months * 100
# 最长连续干旱期
drought_sequence = rainfall < drought_threshold
max_consecutive = 0
current = 0
for is_drought in drought_sequence:
if is_drought:
current += 1
max_consecutive = max(max_consecutive, current)
else:
current = 0
print(f"{region}: 干旱月份={drought_months} ({drought_frequency:.1f}%), "
f"最长连续干旱={max_consecutive}个月")
#干旱地区: 干旱月份=82 (68.3%), 最长连续干旱=13个月
#半干旱地区: 干旱月份=27 (22.5%), 最长连续干旱=4个月
#湿润地区: 干旱月份=1 (0.8%), 最长连续干旱=1个月
#热带雨林: 干旱月份=1 (0.8%), 最长连续干旱=1个月
示例6. 工业生产中的质量控制模拟
np.random.seed(654)
manufacturing_processes = {
"传统加工": {"shape": 2.0, "scale": 0.05}, # 较高变异
"数控加工": {"shape": 3.0, "scale": 0.02}, # 中等变异
"精密加工": {"shape": 4.0, "scale": 0.01}, # 低变异
"超精密加工": {"shape": 5.0, "scale": 0.005} # 极低变异
}
每种工艺生成500个产品的尺寸偏差
n_products = 500
quality_data = {}
for process, params in manufacturing_processes.items():
shape, scale = params["shape"], params["scale"]
quality_data[process] = gamrnd(shape, scale, n_products)
制造过程质量分析:
规格限设定
usl = 0.1 # 上限规格限
lsl = -0.1 # 下限规格限（假设对称）
plt.figure(figsize=(14, 10))
偏差分布图
for i, (process, deviations) in enumerate(quality_data.items()):
plt.subplot(2, 2, i + 1)
# 添加正负偏差（实际制造中偏差可正可负）
symmetric_deviations = deviations * np.random.choice([-1, 1], len(deviations))
plt.hist(symmetric_deviations, bins=50, alpha=0.7)
plt.axvline(lsl, color='red', linestyle='--', label=f'规格限 ±{usl}')
plt.axvline(usl, color='red', linestyle='--')
plt.title(f'{process}尺寸偏差分布')
plt.xlabel('尺寸偏差 (mm)')
plt.ylabel('频数')
plt.legend()
plt.grid(True, alpha=0.3)
# 过程能力指标计算
cpk = min((usl - np.mean(symmetric_deviations)) / (3 * np.std(symmetric_deviations)),
(np.mean(symmetric_deviations) - lsl) / (3 * np.std(symmetric_deviations)))
yield_rate = np.mean((symmetric_deviations >= lsl) & (symmetric_deviations <= usl)) * 100
print(f"{process}: Cpk = {cpk:.2f}, 合格率 = {yield_rate:.1f}%")
#传统加工: Cpk = 0.27, 合格率 = 59.6%
#数控加工: Cpk = 0.48, 合格率 = 87.2%
#精密加工: Cpk = 0.73, 合格率 = 98.6%
#超精密加工: Cpk = 1.20, 合格率 = 100.0%
plt.tight_layout()
plt.show()
批量生产质量预测:
batch_sizes = [100, 1000, 10000]
for process, params in manufacturing_processes.items():
shape, scale = params["shape"], params["scale"]
print(f"\n{process}工艺:")
for batch_size in batch_sizes:
# 模拟批量生产
batch_deviations = gamrnd(shape, scale, batch_size)
symmetric_batch = batch_deviations * np.random.choice([-1, 1], batch_size)
defects = np.sum((symmetric_batch < lsl) | (symmetric_batch > usl))
defect_rate = defects / batch_size * 100
print(f" 批量大小 {batch_size}: 缺陷数={defects}, 缺陷率={defect_rate:.2f}%")
#传统加工工艺:
# 批量大小 100: 缺陷数=41, 缺陷率=41.00%
# 批量大小 1000: 缺陷数=423, 缺陷率=42.30%
# 批量大小 10000: 缺陷数=4082, 缺陷率=40.82%
#数控加工工艺:
# 批量大小 100: 缺陷数=11, 缺陷率=11.00%
# 批量大小 1000: 缺陷数=134, 缺陷率=13.40%
# 批量大小 10000: 缺陷数=1201, 缺陷率=12.01%
#精密加工工艺:
# 批量大小 100: 缺陷数=0, 缺陷率=0.00%
# 批量大小 1000: 缺陷数=11, 缺陷率=1.10%
# 批量大小 10000: 缺陷数=121, 缺陷率=1.21%
#超精密加工工艺:
# 批量大小 100: 缺陷数=0, 缺陷率=0.00%
# 批量大小 1000: 缺陷数=0, 缺陷率=0.00%
# 批量大小 10000: 缺陷数=0, 缺陷率=0.00%
示例7. 高级应用：蒙特卡洛风险模拟
np.random.seed(999)
模拟投资组合的多个风险因素
n_simulations = 10000
市场风险（股票回报）
market_risk = gamrnd(2.0, 0.02, n_simulations) - 0.02 # 中心化
信用风险（违约损失）
credit_risk = gamrnd(1.5, 0.01, n_simulations)
操作风险（意外损失）
operational_risk = gamrnd(2.5, 0.005, n_simulations)
流动性风险（交易成本）
liquidity_risk = gamrnd(3.0, 0.002, n_simulations)
组合风险计算（简化模型）
portfolio_returns = 0.08 + market_risk - credit_risk - operational_risk - liquidity_risk
风险分析
plt.figure(figsize=(15, 10))
回报分布
plt.subplot(2, 2, 1)
plt.hist(portfolio_returns, bins=100, alpha=0.7, density=True)
plt.axvline(np.mean(portfolio_returns), color='red', linestyle='--',
label=f'平均回报: {np.mean(portfolio_returns):.3f}')
plt.title('投资组合回报分布')
plt.xlabel('回报率')
plt.ylabel('密度')
plt.legend()
plt.grid(True, alpha=0.3)
风险贡献分析
risk_factors = {
'市场风险': market_risk,
'信用风险': credit_risk,
'操作风险': operational_risk,
'流动性风险': liquidity_risk
}
risk_contributions = {}
total_variance = np.var(portfolio_returns)
for factor_name, factor_data in risk_factors.items():
# 简化计算风险贡献
correlation = np.corrcoef(factor_data, portfolio_returns)[0, 1]
contribution = correlation ** 2 * np.var(factor_data) / total_variance * 100
risk_contributions[factor_name] = contribution
风险贡献饼图
plt.subplot(2, 2, 2)
plt.pie(risk_contributions.values(), labels=risk_contributions.keys(), autopct='%1.1f%%')
plt.title('风险贡献分析')
极端损失分析
plt.subplot(2, 2, 3)
negative_returns = portfolio_returns[portfolio_returns < 0]
if len(negative_returns) > 0:
plt.hist(-negative_returns, bins=50, alpha=0.7, density=True) # 损失为正数
plt.title('损失分布（负回报）')
plt.xlabel('损失率')
plt.ylabel('密度')
plt.grid(True, alpha=0.3)
风险指标表
plt.subplot(2, 2, 4)
plt.axis('off')
risk_metrics = {
'平均回报': f"{np.mean(portfolio_returns):.4f}",
'回报标准差': f"{np.std(portfolio_returns):.4f}",
'95% VaR': f"{-np.percentile(portfolio_returns, 5):.4f}",
'最大回撤': f"{-np.min(portfolio_returns):.4f}",
'夏普比率': f"{(np.mean(portfolio_returns) - 0.03) / np.std(portfolio_returns):.2f}"
}
table_data = [[metric, value] for metric, value in risk_metrics.items()]
table = plt.table(cellText=table_data,
colLabels=['风险指标', '数值'],
loc='center',
cellLoc='center')
table.auto_set_font_size(False)
table.set_fontsize(10)
table.scale(1, 2)
plt.tight_layout()
plt.show()
蒙特卡洛风险模拟结果:
for metric, value in risk_metrics.items():
print(f"{metric}: {value}")
#平均回报: 0.0668
#回报标准差: 0.0323
#95% VaR: -0.0211
#最大回撤: 0.0325
#夏普比率: 1.14
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def gamma_random_samples(input_str):
"""
根据输入参数生成伽马分布的随机样本，模拟MATLAB的gamrnd函数行为。
参数:
input_str (str): 输入字符串，应能解析为包含形状参数a、比例参数b及可选大小参数的元组。
例如："(2, 3)", "(1, 5, 10)", "([1,2], [3,4], 5, 2)"。
返回:
numpy.ndarray 或 float: 生成的随机样本数组或标量。若输入错误则返回错误信息字符串。
示例:
>>> gamma_random_samples("(2, 3)")
# 生成一个形状参数2、比例参数3的样本
>>> gamma_random_samples("(2, 5, 10)")
# 生成10个形状参数2、比例参数5的样本
>>> gamma_random_samples("([1,2], [3,4], 5, 2)")
# 生成5x2的数组，每个元素基于对应的形状和比例参数
"""
try:
expr = sp.sympify(input_str)
if isinstance(expr, tuple):
if len(expr) >= 2:
a, b = expr[0], expr[1]
size = tuple(expr[2:]) if len(expr) > 2 else None
# 生成伽马分布样本
return np.random.gamma(a, b, size=size)
else:
return f"输入错误: 元组长度不足，至少需要形状参数a和比例参数b。"
else:
return f"输入错误: 输入必须为包含a和b的元组。"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1: 生成单个样本
print("示例1 输出:", gamma_random_samples("(2, 5)"))
# 8.0401321476408
# 示例2: 生成5个样本
print("示例2 输出形状:", gamma_random_samples("(2, 5, 5)"))
# [22.3469496 7.22552538 13.73268475 1.51002786 9.18842614]
# 示例3: 生成5x2矩阵，参数为数组
samples = gamma_random_samples("([1, 2], [3, 4], 5, 2)")
print("示例3 输出形状:", samples)
# [[0.39950036 7.19800049]
# [5.23155375 5.40451038]
# [0.67944897 8.87343055]
# [2.5881992 3.66194587]
# [1.17191571 0.99974115]]
伽玛均值和方差
[M，V]=gamstat(A，B)返回形状参数在A中,尺度参数在B中的伽马分布的均值和方差.A和B可以是具有相同大小的向量.矩阵,这也是M和V的大小.
A或B的标量输入被扩展为与另一输入具有相同维度的常数数组.A和B中的参数必须为正.
参数A和B的伽马分布的平均值为AB.方差为A*B^2.
示例1. 可靠性工程中的设备寿命分析
import numpy as np
import matplotlib.pyplot as plt
不同设备的伽马分布参数
equipment_params = {
"服务器": {"shape": 2.5, "scale": 1000},
"网络设备": {"shape": 3.0, "scale": 800},
"存储设备": {"shape": 4.0, "scale": 600},
"工作站": {"shape": 2.0, "scale": 1200}
}
设备寿命分布统计分析:
results = {}
for equipment, params in equipment_params.items():
shape, scale = params["shape"], params["scale"]
mean, variance = gamstat(shape, scale)
std_dev = np.sqrt(variance)
cv = std_dev / mean # 变异系数
results[equipment] = {
"mean": mean,
"variance": variance,
"std_dev": std_dev,
"cv": cv
}
print(f"{equipment:>10}: 平均寿命={mean:>6.0f}小时, "
f"标准差={std_dev:>5.0f}小时, 变异系数={cv:.3f}")
#服务器: 平均寿命= 2500小时, 标准差= 1581小时, 变异系数=0.632
#网络设备: 平均寿命= 2400小时, 标准差= 1386小时, 变异系数=0.577
#存储设备: 平均寿命= 2400小时, 标准差= 1200小时, 变异系数=0.500
#工作站: 平均寿命= 2400小时, 标准差= 1697小时, 变异系数=0.707
可视化比较
equipment_names = list(equipment_params.keys())
means = [results[eq]["mean"] for eq in equipment_names]
std_devs = [results[eq]["std_dev"] for eq in equipment_names]
plt.figure(figsize=(12, 8))
平均寿命和标准差条形图
plt.subplot(2, 2, 1)
x_pos = np.arange(len(equipment_names))
plt.bar(x_pos, means, yerr=std_devs, capsize=5, alpha=0.7)
plt.title('设备平均寿命与标准差')
plt.ylabel('小时')
plt.xticks(x_pos, equipment_names, rotation=45)
plt.grid(True, alpha=0.3)
变异系数比较
plt.subplot(2, 2, 2)
cvs = [results[eq]["cv"] for eq in equipment_names]
plt.bar(x_pos, cvs, alpha=0.7)
plt.title('设备寿命变异系数')
plt.ylabel('变异系数')
plt.xticks(x_pos, equipment_names, rotation=45)
plt.grid(True, alpha=0.3)
可靠性曲线（生存函数）
plt.subplot(2, 2, 3)
time_points = np.linspace(0, 5000, 100)
for equipment, params in equipment_params.items():
shape, scale = params["shape"], params["scale"]
mean_life = results[equipment]["mean"]
# 简化生存概率计算（基于均值）
survival_probs = np.exp(-time_points / mean_life)
plt.plot(time_points, survival_probs, label=equipment, linewidth=2)
plt.title('设备可靠性曲线（简化模型）')
plt.xlabel('时间 (小时)')
plt.ylabel('生存概率')
plt.legend()
plt.grid(True, alpha=0.3)
维护策略分析
plt.subplot(2, 2, 4)
maintenance_intervals = [1000, 2000, 3000, 4000]
failure_risks = {}
for equipment, params in equipment_params.items():
shape, scale = params["shape"], params["scale"]
mean_life = results[equipment]["mean"]
risks = []
for interval in maintenance_intervals:
# 简化计算：在维护间隔内故障的概率
risk = 1 - np.exp(-interval / mean_life)
risks.append(risk * 100) # 转换为百分比
failure_risks[equipment] = risks
plt.plot(maintenance_intervals, risks, 'o-', label=equipment, markersize=6)
plt.title('不同维护间隔的故障风险')
plt.xlabel('维护间隔 (小时)')
plt.ylabel('故障风险 (%)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
示例2. 保险精算中的索赔分析
不同风险等级的伽马分布参数
risk_profiles = {
"低风险": {"shape": 3.0, "scale": 5000},
"中风险": {"shape": 2.0, "scale": 10000},
"高风险": {"shape": 1.5, "scale": 15000}
}
保险索赔分布统计分析:
claim_stats = {}
for risk_level, params in risk_profiles.items():
shape, scale = params["shape"], params["scale"]
mean, variance = gamstat(shape, scale)
std_dev = np.sqrt(variance)
claim_stats[risk_level] = {
"mean_claim": mean,
"variance": variance,
"std_dev": std_dev,
"cv": std_dev / mean
}
print(f"{risk_level:>8}: 平均索赔=${mean:>8,.0f}, "
f"标准差=${std_dev:>8,.0f}, 变异系数={std_dev / mean:.3f}")
#低风险: 平均索赔=$ 15,000, 标准差=$ 8,660, 变异系数=0.577
#中风险: 平均索赔=$ 20,000, 标准差=$ 14,142, 变异系数=0.707
#高风险: 平均索赔=$ 22,500, 标准差=$ 18,371, 变异系数=0.816
保费定价分析
for risk_level, stats in claim_stats.items():
base_premium = stats["mean_claim"] * 1.2 # 20%的安全边际
risk_load = stats["cv"] * 0.1 # 风险附加费率
total_premium = base_premium * (1 + risk_load)
print(f"{risk_level:>8}: 基础保费=${base_premium:>8,.0f}, "
f"风险附加={risk_load:.1%}, 总保费=${total_premium:>8,.0f}")
#低风险: 基础保费=$ 18,000, 风险附加=5.8%, 总保费=$ 19,039
#中风险: 基础保费=$ 24,000, 风险附加=7.1%, 总保费=$ 25,697
#高风险: 基础保费=$ 27,000, 风险附加=8.2%, 总保费=$ 29,205
风险资本计算
confidence_level = 0.99
for risk_level, stats in claim_stats.items():
# 简化风险资本计算（基于正态近似）
var_amount = stats["mean_claim"] + 2.33 * stats["std_dev"] # 99% VaR
capital_requirement = var_amount * 1.5 # 150%的VaR作为资本要求
print(f"{risk_level:>8}: 99% VaR=${var_amount:>8,.0f}, "
f"风险资本=${capital_requirement:>8,.0f}")
#低风险: 99% VaR=$ 35,178, 风险资本=$ 52,768
#中风险: 99% VaR=$ 52,951, 风险资本=$ 79,427
#高风险: 99% VaR=$ 65,305, 风险资本=$ 97,957
示例3. 医疗统计中的治疗时间分析
medical_procedures = {
"常规检查": {"shape": 3.0, "scale": 10},
"小型手术": {"shape": 2.5, "scale": 20},
"专科诊疗": {"shape": 2.0, "scale": 25},
"急诊处理": {"shape": 1.8, "scale": 15}
}
医疗程序时间分布分析:
treatment_stats = {}
for procedure, params in medical_procedures.items():
shape, scale = params["shape"], params["scale"]
mean, variance = gamstat(f"({shape}, {scale})")
std_dev = np.sqrt(variance)
treatment_stats[procedure] = {
"mean_time": mean,
"std_dev": std_dev,
"cv": std_dev / mean
}
print(f"{procedure:>8}: 平均时间={mean:>5.1f}分钟, "
f"标准差={std_dev:>4.1f}分钟, 变异系数={std_dev / mean:.3f}")
#常规检查: 平均时间= 30.0分钟, 标准差=17.3分钟, 变异系数=0.577
#小型手术: 平均时间= 50.0分钟, 标准差=31.6分钟, 变异系数=0.632
#专科诊疗: 平均时间= 50.0分钟, 标准差=35.4分钟, 变异系数=0.707
#急诊处理: 平均时间= 27.0分钟, 标准差=20.1分钟, 变异系数=0.745
手术室调度优化
working_hours = 8 * 60 # 8小时工作制（分钟）
buffer_factor = 0.2 # 20%的缓冲时间
for procedure, stats in treatment_stats.items():
scheduled_time = stats["mean_time"] * (1 + buffer_factor)
procedures_per_day = working_hours / scheduled_time
print(f"{procedure:>8}: 建议安排时间={scheduled_time:>4.0f}分钟, "
f"每日可安排{procedures_per_day:>3.1f}台")
#常规检查: 建议安排时间= 36分钟, 每日可安排13.3台
#小型手术: 建议安排时间= 60分钟, 每日可安排8.0台
#专科诊疗: 建议安排时间= 60分钟, 每日可安排8.0台
#急诊处理: 建议安排时间= 32分钟, 每日可安排14.8台
患者等待时间预测:
arrival_rate = 4 # 每小时到达4个患者
for procedure, stats in treatment_stats.items():
service_rate = 60 / stats["mean_time"] # 每小时服务患者数
utilization = arrival_rate / service_rate
if utilization < 1:
# 排队论中的平均等待时间（M/G/1队列）
avg_waiting_time = (utilization ** 2 + (stats["cv"]) ** 2) / (2 * service_rate * (1 - utilization)) * 60
print(f"{procedure:>8}: 利用率={utilization:.2f}, 平均等待时间={avg_waiting_time:.1f}分钟")
else:
print(f"{procedure:>8}: 利用率={utilization:.2f} (系统不稳定)")
#常规检查: 利用率=2.00 (系统不稳定)
#小型手术: 利用率=3.33 (系统不稳定)
#专科诊疗: 利用率=3.33 (系统不稳定)
#急诊处理: 利用率=1.80 (系统不稳定)
示例4. 金融风险管理中的波动率分析
asset_classes = {
"国债": {"shape": 4.0, "scale": 0.02},
"蓝筹股": {"shape": 2.5, "scale": 0.05},
"科技股": {"shape": 1.8, "scale": 0.08},
"加密货币": {"shape": 1.2, "scale": 0.15}
}
资产波动率分布分析:
volatility_stats = {}
for asset_class, params in asset_classes.items():
shape, scale = params["shape"], params["scale"]
mean, variance = gamstat(shape, scale)
std_dev = np.sqrt(variance)
volatility_stats[asset_class] = {
"mean_vol": mean,
"std_dev": std_dev,
"cv": std_dev / mean
}
print(f"{asset_class:>8}: 平均波动率={mean:>6.3f}, "
f"标准差={std_dev:>6.3f}, 变异系数={std_dev / mean:.3f}")
#国债: 平均波动率= 0.080, 标准差= 0.040, 变异系数=0.500
#蓝筹股: 平均波动率= 0.125, 标准差= 0.079, 变异系数=0.632
#科技股: 平均波动率= 0.144, 标准差= 0.107, 变异系数=0.745
#加密货币: 平均波动率= 0.180, 标准差= 0.164, 变异系数=0.913
风险调整回报分析
expected_returns = {
"国债": 0.03,
"蓝筹股": 0.08,
"科技股": 0.12,
"加密货币": 0.20
}
for asset_class, stats in volatility_stats.items():
sharpe_ratio = (expected_returns[asset_class] - 0.03) / stats["mean_vol"] # 无风险利率3%
print(f"{asset_class:>8}: 预期回报={expected_returns[asset_class]:.1%}, "
f"夏普比率={sharpe_ratio:.2f}")
#国债: 预期回报=3.0%, 夏普比率=0.00
#蓝筹股: 预期回报=8.0%, 夏普比率=0.40
#科技股: 预期回报=12.0%, 夏普比率=0.62
#加密货币: 预期回报=20.0%, 夏普比率=0.94
投资组合优化
total_investment = 1000000 # 100万美元
for asset_class, stats in volatility_stats.items():
# 基于波动率倒数分配权重（简化方法）
weight = 1 / stats["mean_vol"]
print(f"{asset_class:>8}: 相对权重={weight:.2f}")
#国债: 相对权重=12.50
#蓝筹股: 相对权重=8.00
#科技股: 相对权重=6.94
#加密货币: 相对权重=5.56
归一化权重
total_weight = sum(1 / stats["mean_vol"] for stats in volatility_stats.values())
实际投资金额:
for asset_class, stats in volatility_stats.items():
weight = (1 / stats["mean_vol"]) / total_weight
investment = total_investment * weight
print(f"{asset_class:>8}: 权重={weight:.1%}, 投资金额=${investment:,.0f}")
#国债: 权重=37.9%, 投资金额=$378,788
#蓝筹股: 权重=24.2%, 投资金额=$242,424
#科技股: 权重=21.0%, 投资金额=$210,438
#加密货币: 权重=16.8%, 投资金额=$168,350
示例5. 环境科学中的降雨量分析
climate_regions = {
"干旱地区": {"shape": 1.8, "scale": 15},
"半干旱地区": {"shape": 2.2, "scale": 30},
"湿润地区": {"shape": 2.5, "scale": 60},
"热带雨林": {"shape": 3.0, "scale": 100}
}
降雨量分布分析:
rainfall_stats = {}
for region, params in climate_regions.items():
shape, scale = params["shape"], params["scale"]
mean, variance = gamstat(shape, scale)
std_dev = np.sqrt(variance)
rainfall_stats[region] = {
"mean_rainfall": mean,
"std_dev": std_dev,
"cv": std_dev / mean
}
print(f"{region:>8}: 平均降雨量={mean:>5.1f}mm, "
f"标准差={std_dev:>4.1f}mm, 变异系数={std_dev / mean:.3f}")
#干旱地区: 平均降雨量= 27.0mm, 标准差=20.1mm, 变异系数=0.745
#半干旱地区: 平均降雨量= 66.0mm, 标准差=44.5mm, 变异系数=0.674
#湿润地区: 平均降雨量=150.0mm, 标准差=94.9mm, 变异系数=0.632
#热带雨林: 平均降雨量=300.0mm, 标准差=173.2mm, 变异系数=0.577
干旱风险评估
drought_thresholds = {
"轻度干旱": 0.5, # 平均降雨量的50%
"中度干旱": 0.3, # 平均降雨量的30%
"严重干旱": 0.1 # 平均降雨量的10%
}
for region, stats in rainfall_stats.items():
print(f"\n{region}:")
for drought_level, threshold_ratio in drought_thresholds.items():
threshold = stats["mean_rainfall"] * threshold_ratio
# 简化干旱概率计算（基于正态近似）
z_score = (threshold - stats["mean_rainfall"]) / stats["std_dev"]
drought_probability = max(0, 0.5 * (1 + sp.erf(z_score / np.sqrt(2))))
print(f" {drought_level:>10}: 阈值={threshold:>4.1f}mm, "
f"发生概率={drought_probability:.1%}")
#干旱地区:
# 轻度干旱: 阈值=13.5mm, 发生概率=25.1%
# 中度干旱: 阈值= 8.1mm, 发生概率=17.4%
# 严重干旱: 阈值= 2.7mm, 发生概率=11.4%
#半干旱地区:
# 轻度干旱: 阈值=33.0mm, 发生概率=22.9%
# 中度干旱: 阈值=19.8mm, 发生概率=15.0%
# 严重干旱: 阈值= 6.6mm, 发生概率=9.1%
#湿润地区:
# 轻度干旱: 阈值=75.0mm, 发生概率=21.5%
# 中度干旱: 阈值=45.0mm, 发生概率=13.4%
# 严重干旱: 阈值=15.0mm, 发生概率=7.7%
#热带雨林:
# 轻度干旱: 阈值=150.0mm, 发生概率=19.3%
# 中度干旱: 阈值=90.0mm, 发生概率=11.3%
# 严重干旱: 阈值=30.0mm, 发生概率=6.0%
水资源规划
population_density = {
"干旱地区": 50, # 人/平方公里
"半干旱地区": 100, # 人/平方公里
"湿润地区": 200, # 人/平方公里
"热带雨林": 150 # 人/平方公里
}
water_requirement = 100 # 每人每年100立方米
for region, stats in rainfall_stats.items():
# 简化水资源计算（假设10%的降雨可收集利用）
collectable_water = stats["mean_rainfall"] * 0.1 # mm -> 立方米/公顷
population_support = (collectable_water * 10000) / (water_requirement * population_density[region])
print(f"{region:>8}: 可收集水量={collectable_water:.0f}立方米/公顷, "
f"可支持人口密度={population_support:.0f}人/平方公里")
#干旱地区: 可收集水量=3立方米/公顷, 可支持人口密度=5人/平方公里
#半干旱地区: 可收集水量=7立方米/公顷, 可支持人口密度=7人/平方公里
#湿润地区: 可收集水量=15立方米/公顷, 可支持人口密度=8人/平方公里
#热带雨林: 可收集水量=30立方米/公顷, 可支持人口密度=20人/平方公里
示例6. 矩阵输入的批量计算示例
创建测试数据矩阵
shape_params = [[1.5, 2.0, 2.5], [3.0, 3.5, 4.0]] # 形状参数矩阵
scale_params = [[10, 15, 20], [25, 30, 35]] # 尺度参数矩阵
形状参数矩阵:
print(np.array(shape_params))
#[[1.5 2. 2.5]
[3. 3.5 4. ]]
尺度参数矩阵:
print(np.array(scale_params))
#[[10 15 20]
[25 30 35]]
批量计算均值和方差
means, variances = gamstat(shape_params, scale_params)
print("\n均值矩阵:")
print(means)
print("\n方差矩阵:")
print(variances)
#均值矩阵:
#[[ 15. 30. 50.]
[ 75. 105. 140.]]
#方差矩阵:
#[[ 150. 450. 1000.]
[1875. 3150. 4900.]]
计算标准差和变异系数
std_devs = np.sqrt(variances)
coefficients_of_variation = std_devs / means
print("\n标准差矩阵:")
print(std_devs)
print("\n变异系数矩阵:")
print(coefficients_of_variation)
#标准差矩阵:
#[[12.24744871 21.21320344 31.6227766 ]
[43.30127019 56.1248608 70. ]]
#变异系数矩阵:
#[[0.81649658 0.70710678 0.63245553]
[0.57735027 0.53452248 0.5 ]]
实际应用：批量风险评估
risk_categories = ["低风险", "中风险", "高风险", "极高风险"]
for i in range(means.shape[0]):
for j in range(means.shape[1]):
mean_val = means[i, j]
cv_val = coefficients_of_variation[i, j]
# 基于均值和变异系数的风险评估
if mean_val < 20 and cv_val < 0.5:
risk_level = "低风险"
elif mean_val < 50 and cv_val < 0.8:
risk_level = "中风险"
elif mean_val < 100 and cv_val < 1.2:
risk_level = "高风险"
else:
risk_level = "极高风险"
print(f"位置({i},{j}): 均值={mean_val:.1f}, 变异系数={cv_val:.3f} -> {risk_level}")
#位置(0,0): 均值=15.0, 变异系数=0.816 -> 高风险
#位置(0,1): 均值=30.0, 变异系数=0.707 -> 中风险
#位置(0,2): 均值=50.0, 变异系数=0.632 -> 高风险
#位置(1,0): 均值=75.0, 变异系数=0.577 -> 高风险
#位置(1,1): 均值=105.0, 变异系数=0.535 -> 极高风险
#位置(1,2): 均值=140.0, 变异系数=0.500 -> 极高风险
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Gamma, E, variance
import numpy as np
def gamstat(k, theta):
"""
计算Gamma分布的均值和方差，对标MATLAB的gamstat函数
参数:
k : 形状参数（shape parameter），可以是标量、向量或多维数组
theta : 尺度参数（scale parameter），可以是标量、向量或多维数组
返回:
(mean, var) : 包含均值和方差的元组，形状与广播后的输入一致
数学公式:
均值 = k * theta
方差 = k * theta^2
"""
# 转换为NumPy数组
k_arr = np.asarray(k, dtype=float)
theta_arr = np.asarray(theta, dtype=float)
# 检查参数有效性
if np.any(k_arr <= 0) or np.any(theta_arr <= 0):
raise ValueError("参数k和theta必须大于0")
# 使用广播机制计算
mean = k_arr * theta_arr
var = k_arr * (theta_arr ** 2)
return mean, var
def gamma_mean_variance(input_str):
"""
计算Gamma分布的均值和方差，对标MATLAB的gamstat函数。
参数:
input_str (str): 输入字符串，应能解析为包含形状参数a和比例参数b的元组。
支持标量、列表或矩阵形式，例如："(2, 3)", "([1,2], [3,4])"。
返回:
tuple: 包含均值矩阵和方差矩阵的元组，若输入错误则返回错误信息字符串。
示例:
>>> gamma_mean_variance("(2, 3)")
(6, 18)
>>> gamma_mean_variance("([1, 2], [3, 4])")
(Matrix([[3, 8]]), Matrix([[9, 32]]))
>>> gamma_mean_variance("([1, 2], 2)")
(Matrix([[2, 4]]), Matrix([[4, 16]]))
"""
try:
expr = sp.sympify(input_str)
M, V = None, None # 初始化均值和方差结果
# 定义符号参数
a_sym = sp.Symbol("a", positive=True)
b_sym = sp.Symbol("b", positive=True)
X = Gamma("X", a_sym, b_sym) # 定义Gamma随机变量
# 符号表达式计算均值和方差
mean_expr = E(X)
var_expr = variance(X)
def eval_mean_var(a_val, b_val):
"""代入数值计算均值和方差"""
return (
mean_expr.subs({a_sym: a_val, b_sym: b_val}).evalf(),
var_expr.subs({a_sym: a_val, b_sym: b_val}).evalf()
)
# 检查输入是否为二元组
if not isinstance(expr, tuple) or len(expr) != 2:
return f"输入错误: 需要形状参数a和比例参数b的二元组，当前输入: {input_str}"
# 情况4: a和b均为标量
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
M, V = gamstat(*params)
elif any(e.free_symbols for e in expr):
M, V = eval_mean_var(*expr)
else:
raise ValueError(f"输入错误:{input_str}")
return (M, V)
except Exception as e:
return f"错误: {str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1: 标量参数
print("示例1:", gamma_mean_variance("(2, 3)"))
# (6.0, 18.0)
广义自回归条件异方差模型拟合
是特别用来建立条件方差模型并对其进行预测的.
这些模型被运用于经济学领域，尤其是金融时间序列分析中.
示例1. 股票收益率波动性分析
模拟典型的股票日收益率数据（百分比形式）
stock_returns = [
0.52, -1.23, 0.78, 2.15, -0.89, -2.34, 1.67, 0.45, -1.56, 0.92,
-0.67, 1.34, -2.01, 0.23, 1.89, -1.12, 0.56, -0.78, 2.34, -1.45,
0.67, -2.12, 1.56, -0.89, 0.34, 1.23, -1.78, 0.45, -0.92, 2.67
]
stock_result = GarchFit(stock_returns)
print(f"参数估计: {stock_result['参数估计']}")
print(f"模型类型: {stock_result['模型类型']}")
print(f"AIC: {stock_result['AIC']:.4f}")
#参数估计: {'均值': 6.699949696678251, '常数项': 10037.738334961734, 'ARCH(1)系数': 0.0, 'GARCH(1)系数': 0.5092528370787296}
#模型类型: GARCH
#AIC: 390.3290
示例2. 汇率波动建模
模拟汇率日变化率（特征：持续性波动）
exchange_rate_changes = [
0.12, -0.08, 0.05, -0.15, 0.03, 0.07, -0.02, -0.11, 0.09, -0.06,
0.04, 0.01, -0.13, 0.08, -0.04, 0.11, -0.07, 0.02, -0.09, 0.06,
-0.03, 0.10, -0.05, 0.07, -0.01, -0.12, 0.04, 0.09, -0.08, 0.03
]
fx_result = GarchFit(exchange_rate_changes)
print(f"ARCH系数: {fx_result['参数估计'].get('ARCH(1)系数', 'N/A')}")
print(f"GARCH系数: {fx_result['参数估计'].get('GARCH(1)系数', 'N/A')}")
#ARCH系数: 6.219731413742339e-14
#GARCH系数: 0.9887817961831756
示例3. 商品价格波动分析
模拟商品价格收益率（特征：大幅波动）
commodity_returns = [
1.5, -2.3, 3.1, -1.8, 0.9, -3.2, 2.4, -0.7, 1.8, -2.1,
3.5, -1.2, 0.6, -2.8, 1.9, -0.5, 2.7, -1.6, 0.8, -3.1,
2.2, -0.9, 1.4, -2.5, 3.3, -1.1, 0.7, -2.9, 2.1, -0.8
]
commodity_result = GarchFit(commodity_returns)
print(f"条件波动率范围: {min(commodity_result['条件波动率']):.4f} - {max(commodity_result['条件波动率']):.4f}")
#条件波动率范围: 205.5241 - 211.8954
示例4. 市场恐慌期波动率建模
模拟市场恐慌期的收益率（特征：极端波动）
crisis_returns = [
-0.5, -1.2, -3.8, 2.1, -4.5, 1.8, -5.2, -2.7, 3.4, -6.1,
0.9, -3.9, -1.6, 4.2, -2.8, -5.7, 1.3, -4.8, -0.7, 2.9,
-3.5, 0.6, -2.4, -4.1, 1.7, -5.9, -1.8, 3.2, -4.7, 0.8
]
crisis_result = GarchFit(crisis_returns)
print(f"波动率持续性: {crisis_result['参数估计'].get('GARCH(1)系数', 0) + crisis_result['参数估计'].get('ARCH(1)系数', 0):.4f}")
#波动率持续性: 0.7127
示例5. 低波动环境建模
模拟低波动环境的收益率（特征：小幅波动）
low_vol_returns = [
0.1, -0.2, 0.3, -0.1, 0.2, -0.3, 0.1, 0.0, -0.2, 0.3,
-0.1, 0.2, 0.0, -0.3, 0.1, -0.2, 0.3, -0.1, 0.2, -0.3,
0.1, 0.0, -0.2, 0.3, -0.1, 0.2, 0.0, -0.3, 0.1, -0.2
]
low_vol_result = GarchFit(low_vol_returns)
print(f"常数项(omega): {low_vol_result['参数估计'].get('常数项', 'N/A')}")
#常数项(omega): 194.94450723700578
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from arch import arch_model
def garchfit_regression_expression(input_str,model_type='Garch', p=1, q=1, dist='normal', update_freq=0):
"""
对输入的时间序列数据进行 ARCH 模型拟合。
参数:
input_str: 输入的字符串表达式，表示时间序列数据（矩阵或列表）。
返回:
如果输入有效，则返回 ARCH 模型的拟合参数（以中文命名）；否则返回错误信息。
"""
try:
# 将输入字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 检查输入是否为矩阵
if isinstance(expr, list):
# 将 SymPy 矩阵转换为 NumPy 数组并展平
A = np.array(expr, dtype=float).ravel()
else:
return f"输入错误: {input_str}"
# 使用 ARCH 模型进行拟合
# 假设 p=1（ARCH(1) 模型），vol='ARCH'（ARCH 模型）
am = arch_model(A * 100, vol='Garch', p=p, dist=dist)
# 拟合模型，设置 update_freq=0 和 disp='off' 以关闭输出
res = am.fit(update_freq=update_freq, disp=update_freq > 0)
# 获取拟合参数
params_dict = dict(res.params)
# 创建一个新的字典，用于存储重命名后的参数
renamed_params_dict = {}
# 遍历原始字典，将键重命名为更具描述性的中文名称
# 根据模型类型重命名参数
for key, value in params_dict.items():
if key == 'mu':
renamed_params_dict['均值'] = value
elif key == 'omega':
renamed_params_dict['常数项'] = value
elif key.startswith('alpha'):
idx = key.split('[')[1].split(']')[0] if '[' in key else '1'
renamed_params_dict[f'ARCH({idx})系数'] = value
elif key.startswith('beta'):
idx = key.split('[')[1].split(']')[0] if '[' in key else '1'
renamed_params_dict[f'GARCH({idx})系数'] = value
elif key == 'nu' and dist == 't':
renamed_params_dict['自由度'] = value
elif key == 'lambda' and dist == 'skewt':
renamed_params_dict['偏度参数'] = value
else:
renamed_params_dict[key] = value
# 计算一些有用的统计量
volatility = res.conditional_volatility
residuals = res.resid
# 返回结果字典
result = {
'参数估计': renamed_params_dict,
'模型类型': model_type,
'阶数': f"ARCH({p})" if model_type == 'ARCH' else f"GARCH({p},{q})",
'误差分布': dist,
'对数似然值': res.loglikelihood,
'AIC': res.aic,
'BIC': res.bic,
'条件波动率': volatility.tolist(),
'残差': residuals.tolist(),
'模型摘要': str(res.summary())
}
return result
except Exception as e:
return f"错误: {e}"
def main():
"""
主函数，用于演示 archfit_regression_expression 函数的使用。
"""
# 示范代码
input_examples = [
"[[1, 2, 3, 4, 5]]", # 时间序列数据
# 结果: {'均值': 2.999870510602842, '方差偏移': 1.9871738584685994, '方差系数': 0.00810550114746725}
"[[0.1, 0.2, 0.3, 0.4, 0.5]]", # 时间序列数据
# 结果: {'均值': 0.2999917328352241, '方差偏移': 0.019852350048059936, '方差系数': 0.00921009134817112}
]
for input_str in input_examples:
print(f"输入: {input_str}")
result = garchfit_regression_expression(input_str)
print(f"结果: {result}")
print("-" * 40)
if __name__ == "__main__":
main()
高斯拟合
gauss1: Y = a1*exp(-((x-b1)/c1)^2)
gauss2: Y = a1*exp(-((x-b1)/c1)^2)+a2*exp(-((x-b2)/c2)^2)
gauss3: Y = a1*exp(-((x-b1)/c1)^2)+...+a3*exp(-((x-b3)/c3)^2)
gauss8: Y = a1*exp(-((x-b1)/c1)^2)+...+a8*exp(-((x-b8)/c8)^2)
以此类推，最大到 gauss8
示例1. 光谱分析 - 原子发射光谱
模拟原子发射光谱数据（如钠双线）
模拟钠的D线：589.0nm和589.6nm
wavelength = np.linspace(580, 600, 200)
# 两个高斯峰叠加，模拟钠双线
intensity = (2.5 * np.exp(-((wavelength - 589.0) / 0.3) ** 2) +
2.0 * np.exp(-((wavelength - 589.6) / 0.25) ** 2) +
np.random.normal(0, 0.05, len(wavelength)))
result = gaussfit(wavelength, intensity, 2)
原子发射光谱拟合结果:
print(f"拟合函数: {result}")
#拟合函数: 2.0296*exp(-24176.6183185059*(0.00222861984692947*x - 1)**2) + 1.9706*exp(-64682.0349195574*(0.00136322559897747*x - 1)**2)
绘制结果
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.plot(wavelength, intensity, 'b.', label='实验数据', alpha=0.6)
计算拟合曲线
x_fit = np.linspace(580, 600, 400)
y_fit = [float(result.subs('x', xi)) for xi in x_fit]
plt.plot(x_fit, y_fit, 'r-', label='高斯拟合', linewidth=2)
plt.xlabel('波长 (nm)')
plt.ylabel('强度')
plt.title('钠原子发射光谱的高斯拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
示例2. 物理实验 - X射线衍射峰
模拟X射线衍射数据，晶体结构分析
two_theta = np.linspace(10, 80, 500)
模拟多晶材料的衍射峰
intensity = (8.0 * np.exp(-((two_theta - 28.4) / 0.4) ** 2) +
6.5 * np.exp(-((two_theta - 32.8) / 0.35) ** 2) +
4.2 * np.exp(-((two_theta - 47.5) / 0.5) ** 2) +
3.8 * np.exp(-((two_theta - 56.2) / 0.45) ** 2) +
2.1 * np.exp(-((two_theta - 68.7) / 0.6) ** 2) +
np.random.normal(0, 0.2, len(two_theta)))
result = gaussfit(two_theta, intensity, 5)
X射线衍射分析拟合结果:
print(f"拟合函数: {result}")
#拟合函数: 7.8525*exp(-4687.60959974145*(0.0352200726237897*x - 1)**2) +
6.522*exp(-8653.06115141254*(0.0304796274170345*x - 1)**2) +
4.232*exp(-9901.89693628553*(0.0210547148875783*x - 1)**2) +
3.5839*exp(-13661.0783125254*(0.0177947975129991*x - 1)**2) +
1.9857*exp(-11578.3865102754*(0.0145459834903087*x - 1)**2)
示例3. 生物医学 - 蛋白质电泳条带分析
模拟蛋白质电泳条带的强度分布
migration_distance = np.linspace(0, 100, 200)
模拟不同分子量蛋白质的条带
band_intensity = (15 * np.exp(-((migration_distance - 25) / 4) ** 2) +
12 * np.exp(-((migration_distance - 45) / 5) ** 2) +
8 * np.exp(-((migration_distance - 65) / 6) ** 2) +
20 * np.exp(-((migration_distance - 80) / 3) ** 2) +
np.random.normal(0, 0.3, len(migration_distance)))
result = gaussfit(migration_distance, band_intensity, 4)
蛋白质电泳条带分析结果:
print(f"拟合函数: {result}")
#拟合函数: 6.1472*exp(-1.7122351336583*(0.0182527716833801*x - 1)**2) +
1.4598*exp(-24.2254686864556*(0.00824973600844773*x + 1)**2) +
14.9839*exp(-621.758688213008*(0.00541764635771635*x - 1)**2) +
17.188*exp(-5807.31813286032*(0.00224091795170239*x - 1)**2)
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.optimize import curve_fit
def gauss_fit_nonlinear(input_str):
"""
使用非线性最小二乘法进行高斯曲线拟合，支持1到8阶高斯函数。
参数:
input_str: 字符串形式的输入数据，格式为"(x_data, y_data, n)"，其中
x_data 和 y_data 是数据点，n 是要拟合的高斯函数阶数。
返回:
拟合后的高斯函数表达式（SymPy 表达式）或错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
maxfev = 10000 # 最大函数评估次数
# 定义1到8阶高斯函数
def gauss1(x, a1, b1, c1):
return a1 * np.exp(-((x - b1) / c1) ** 2)
def gauss2(x, a1, b1, c1, a2, b2, c2):
return gauss1(x, a1, b1, c1) + gauss1(x, a2, b2, c2)
def gauss3(x, a1, b1, c1, a2, b2, c2, a3, b3, c3):
return gauss2(x, a1, b1, c1, a2, b2, c2) + gauss1(x, a3, b3, c3)
def gauss4(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4):
return gauss3(x, a1, b1, c1, a2, b2, c2, a3, b3, c3) + gauss1(x, a4, b4, c4)
def gauss5(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4, a5, b5, c5):
return gauss4(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4) + gauss1(x, a5, b5, c5)
def gauss6(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4, a5, b5, c5, a6, b6, c6):
return gauss5(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4, a5, b5, c5) + gauss1(x, a6, b6, c6)
def gauss7(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4, a5, b5, c5, a6, b6, c6, a7, b7, c7):
return gauss6(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4, a5, b5, c5, a6, b6, c6) + gauss1(x, a7, b7,
c7)
def gauss8(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4, a5, b5, c5, a6, b6, c6, a7, b7, c7, a8, b8, c8):
return gauss7(x, a1, b1, c1, a2, b2, c2, a3, b3, c3, a4, b4, c4, a5, b5, c5, a6, b6, c6, a7, b7,
c7) + gauss1(x, a8, b8, c8)
# 解析输入数据
if isinstance(expr, tuple):
if len(expr) > 2:
x, y, n = expr[0], expr[1], int(expr[2])
else:
x, y = expr[0], expr[1]
n = 1
# 转换为NumPy数组
if x is None or y is None:
raise ValueError("Invalid input format")
x_data = np.array(x, dtype=float).ravel()
y_data = np.array(y, dtype=float).ravel()
else:
raise ValueError("Input must be a tuple")
# 参数边界设置 (a≥0, b∈ℝ, c≥0)
lower_bounds = [0, -np.inf, 0] * n
upper_bounds = [np.inf, np.inf, np.inf] * n
param_bounds = (lower_bounds, upper_bounds)
# 生成初始猜测参数
initial_guess = []
sorted_indices = np.argsort(x_data)
sorted_x = x_data[sorted_indices]
sorted_y = y_data[sorted_indices]
# 通过寻找峰值改进初始猜测
peak_indices = np.argsort(sorted_y)[-n:][::-1] # 取最大的n个y值对应的索引
peak_x = sorted_x[peak_indices]
peak_x.sort()
for i in range(n):
a_guess = sorted_y[peak_indices[i]] - np.min(sorted_y)
b_guess = peak_x[i]
c_guess = (np.max(sorted_x) - np.min(sorted_x)) / (2 * n)
initial_guess.extend([a_guess, b_guess, c_guess])
# 选择对应阶数的高斯函数进行拟合
gauss_funcs = [None, gauss1, gauss2, gauss3, gauss4,
gauss5, gauss6, gauss7, gauss8]
if n < 1 or n > 8:
raise ValueError("n must be between 1 and 8")
params, _ = curve_fit(gauss_funcs[n], x_data, y_data,
p0=initial_guess, maxfev=maxfev,
bounds=param_bounds)
# 生成拟合表达式
terms = []
for i in range(n):
a, b, c = params[3 * i: 3 * i + 3]
terms.append(f"{a:.4f}*exp(-((x - {b:.4f})/{c:.4f})**2)")
return sp.sympify("+".join(terms))
except Exception as e:
return f"Error: {str(e)}"
# 示例使用
if __name__ == "__main__":
# 生成测试数据（两个高斯峰叠加）
x = np.linspace(0, 10, 100)
y = 3 * np.exp(-(x - 2) ** 2 / 1) + 2 * np.exp(-(x - 7) ** 2 / 0.5) + np.random.normal(0, 0.1, 100)
# 构造输入字符串（注意需要转换为Python列表）
input_str = f"({x.tolist()}, {y.tolist()}, 2)"
print(input_str)
# 进行二阶高斯拟合
result = gauss_fit_nonlinear(input_str)
print("拟合结果:", result)
# 3.0246*exp(-3.9764864633667*(0.500375281461096*x - 1)**2) + 2.0451*exp(-93.9532479488075*(0.143209027897119*x - 1)**2)
高斯窗
w = gausswin(L,std,sym=1) 返回一个长度为L个点的修正的高斯窗 默认sym=1
当sym=1, 生成一个对称窗口，用于滤波器设计.
当sym=0, 生成一个周期性窗口，用于光谱分析.
示例1. 信号处理 - 频谱分析中的加窗
在频谱分析中使用高斯窗减少频谱泄漏
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
生成测试信号：两个频率分量
fs = 1000 # 采样频率
t = np.linspace(0, 1, fs, endpoint=False)
signal_data = (1.0 * np.sin(2 * np.pi * 50 * t) +
0.5 * np.sin(2 * np.pi * 120 * t) +
0.2 * np.random.normal(0, 0.1, len(t)))
应用不同参数的高斯窗
windows = {
'矩形窗(无窗)': np.ones(len(signal_data)),
'高斯窗(α=2.5)': gausswin(len(t),2.5),
'高斯窗(α=1.0)': gausswin(len(t),1.0),
'高斯窗(α=0.4)': gausswin(len(t),0.4)
}
plt.figure(figsize=(12, 8))
for i, (name, window) in enumerate(windows.items()):
# 加窗
windowed_signal = signal_data * window
# 计算FFT
fft_vals = fft(windowed_signal)
freqs = fftfreq(len(t), 1 / fs)
# 绘制频谱
plt.subplot(2, 2, i + 1)
plt.plot(freqs[:len(freqs) // 2], np.abs(fft_vals[:len(fft_vals) // 2]))
plt.title(f'{name} - 频谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅度')
plt.grid(True, alpha=0.3)
# 标记主要频率分量
plt.axvline(x=50, color='r', linestyle='--', alpha=0.7, label='50Hz')
plt.axvline(x=120, color='g', linestyle='--', alpha=0.7, label='120Hz')
plt.legend()
plt.tight_layout()
plt.show()
示例2. 图像处理 - 高斯模糊滤波器
import matplotlib.pyplot as plt
from scipy import ndimage
创建测试图像（简单的几何形状）
image = np.zeros((100, 100))
在图像中心画一个矩形
image[30:70, 30:70] = 1.0
添加一些噪声
image += 0.1 * np.random.normal(0, 1, image.shape)
生成2D高斯窗（可分离滤波器）
window_size = 15
gauss_1d = gausswin(window_size,2.5)
创建2D高斯滤波器
gauss_2d = np.outer(gauss_1d, gauss_1d)
gauss_2d = gauss_2d / np.sum(gauss_2d) # 归一化
应用高斯模糊
blurred_image = ndimage.convolve(image, gauss_2d, mode='constant', cval=0.0)
显示结果
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
axes[0].imshow(image, cmap='gray')
axes[0].set_title('原始图像')
axes[0].axis('off')
axes[1].imshow(gauss_2d, cmap='hot')
axes[1].set_title('高斯滤波器')
axes[1].axis('off')
axes[2].imshow(blurred_image, cmap='gray')
axes[2].set_title('高斯模糊后')
axes[2].axis('off')
plt.tight_layout()
plt.show()
示例3. 音频处理 - 语音信号的短时傅里叶变换
import matplotlib.pyplot as plt
生成模拟语音信号（啁啾信号）
fs = 8000 # 采样率
duration = 2 # 秒
t = np.linspace(0, duration, int(fs * duration))
创建频率变化的信号
freq_sweep = 100 + 200 * t # 频率从100Hz线性增加到500Hz
audio_signal = np.sin(2 * np.pi * freq_sweep * t)
添加一些突发信号模拟语音特性
audio_signal[5000:5200] += 0.5 * np.sin(2 * np.pi * 800 * t[:200])
audio_signal[12000:12200] += 0.3 * np.sin(2 * np.pi * 1200 * t[:200])
STFT参数
window_size = 256
hop_size = 64
生成高斯窗
gauss_window = gausswin(window_size,3.0)
计算STFT
def compute_stft(signal, window, hop):
num_frames = (len(signal) - window_size) // hop + 1
stft_matrix = np.zeros((window_size // 2 + 1, num_frames), dtype=complex)
for i in range(num_frames):
start = i * hop
frame = signal[start:start + window_size] * window
spectrum = np.fft.rfft(frame)
stft_matrix[:, i] = spectrum
return stft_matrix
stft_result = compute_stft(audio_signal, gauss_window, hop_size)
绘制时频谱
plt.figure(figsize=(12, 8))
plt.subplot(2, 1, 1)
plt.plot(t, audio_signal)
plt.title('原始音频信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅度')
plt.grid(True, alpha=0.3)
plt.subplot(2, 1, 2)
显示频谱图
spectrogram = 20 * np.log10(np.abs(stft_result) + 1e-10)
plt.imshow(spectrogram, aspect='auto', origin='lower',
extent=[0, duration, 0, fs / 2], cmap='viridis')
plt.title('STFT时频谱（高斯窗）')
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')
plt.colorbar(label='幅度 (dB)')
plt.tight_layout()
plt.show()
示例4. 数字滤波器设计 - FIR滤波器
import matplotlib.pyplot as plt
滤波器规格
filter_order = 64 # 滤波器阶数
cutoff_freq = 0.2 # 归一化截止频率 (0.5对应奈奎斯特频率)
设计理想低通滤波器的冲激响应
n = np.arange(-filter_order // 2, filter_order // 2 + 1)
ideal_response = 2 * cutoff_freq * np.sinc(2 * cutoff_freq * n)
应用不同参数的高斯窗
windows = {
'矩形窗': np.ones(len(n)),
'高斯窗(α=2.5)': gausswin(len(n),2.5),
'高斯窗(α=1.5)': gausswin(len(n),1.5),
'高斯窗(α=0.8)': gausswin(len(n),0.8)
}
plt.figure(figsize=(12, 10))
for i, (name, window) in enumerate(windows.items()):
# 加窗得到实际滤波器系数
fir_coeff = ideal_response * window
# 计算频率响应
w, h = signal.freqz(fir_coeff)
freq = w / np.pi # 归一化频率
# 绘制频率响应
plt.subplot(2, 2, i + 1)
plt.plot(freq, 20 * np.log10(np.abs(h) + 1e-10))
plt.axvline(x=cutoff_freq, color='r', linestyle='--', alpha=0.7, label='截止频率')
plt.title(f'{name} - 频率响应')
plt.xlabel('归一化频率 (×π rad/sample)')
plt.ylabel('幅度 (dB)')
plt.grid(True, alpha=0.3)
plt.legend()
plt.ylim([-100, 5])
plt.tight_layout()
plt.show()
示例5. 数据平滑 - 时间序列滤波
import matplotlib.pyplot as plt
生成带噪声的时间序列数据（模拟股票价格）
np.random.seed(42)
t = np.linspace(0, 100, 500)
基础趋势 + 季节性 + 噪声
trend = 0.02 * t
seasonal = 2 * np.sin(2 * np.pi * t / 20)
noise = 0.5 * np.random.normal(0, 1, len(t))
time_series = 50 + trend + seasonal + noise
创建不同宽度的高斯平滑滤波器
window_sizes = [10, 25, 50]
plt.figure(figsize=(12, 8))
plt.plot(t, time_series, 'b-', alpha=0.3, label='原始数据', linewidth=1)
for window_size in window_sizes:
# 生成高斯窗
gauss_window = gausswin(window_size,2.0)
gauss_window = gauss_window / np.sum(gauss_window) # 归一化
# 应用卷积进行平滑
smoothed = np.convolve(time_series, gauss_window, mode='same')
plt.plot(t, smoothed, linewidth=2,
label=f'高斯平滑 (窗口大小={window_size})')
plt.title('时间序列数据的高斯平滑')
plt.xlabel('时间')
plt.ylabel('数值')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
示例6. 峰值检测 - 信号预处理
import matplotlib.pyplot as plt
生成带多个峰值的信号
t = np.linspace(0, 10, 1000)
创建多个高斯峰叠加的信号
peaks = [
(1.0, 2.0, 0.3), # (幅度, 位置, 宽度)
(0.8, 4.0, 0.2),
(1.2, 6.0, 0.4),
(0.9, 8.0, 0.25)
]
clean_signal = np.zeros_like(t)
for amp, pos, width in peaks:
clean_signal += amp * np.exp(-((t - pos) / width) ** 2)
# 添加噪声
noisy_signal = clean_signal + 0.1 * np.random.normal(0, 1, len(t))
使用高斯窗进行平滑以改善峰值检测
gauss_window = gausswin(51,1.5)
gauss_window = gauss_window / np.sum(gauss_window)
smoothed_signal = np.convolve(noisy_signal, gauss_window, mode='same')
峰值检测
peaks_clean, _ = signal.find_peaks(clean_signal, height=0.5)
peaks_noisy, _ = signal.find_peaks(noisy_signal, height=0.5)
peaks_smoothed, _ = signal.find_peaks(smoothed_signal, height=0.5)
绘制结果
plt.figure(figsize=(12, 8))
plt.subplot(3, 1, 1)
plt.plot(t, clean_signal, 'b-', label='原始信号')
plt.plot(t[peaks_clean], clean_signal[peaks_clean], 'ro', label='检测到的峰值')
plt.title('无噪声信号峰值检测')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(3, 1, 2)
plt.plot(t, noisy_signal, 'g-', label='带噪声信号', alpha=0.7)
plt.plot(t[peaks_noisy], noisy_signal[peaks_noisy], 'ro', label='检测到的峰值')
plt.title('带噪声信号峰值检测')
plt.legend()
plt.grid(True, alpha=0.3)
plt.subplot(3, 1, 3)
plt.plot(t, smoothed_signal, 'm-', label='高斯平滑后信号')
plt.plot(t[peaks_smoothed], smoothed_signal[peaks_smoothed], 'ro', label='检测到的峰值')
plt.title('高斯平滑后峰值检测')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.signal as signal
def gausswin_window(input_str):
"""
对标MATLAB的gausswin函数，生成高斯窗口。
参数：
input_str (str): 输入的参数字符串，格式应为 "M", "M, alpha" 或 "M, alpha, sym_flag"
- M: 窗口长度（正整数）
- alpha: 高斯窗的alpha参数（正数，默认2.5），对应MATLAB中的alpha，与标准差σ的关系为σ = (M-1)/(2*alpha)
- sym_flag: 可选参数，0表示非对称窗口（sym=False），其他值表示对称窗口（默认，sym=True）
返回：
list: 高斯窗口数组，或错误信息字符串。
"""
try:
expr = sp.sympify(input_str)
m = None
alpha = 2.5 # MATLAB默认alpha为2.5
sym = True # 默认对称窗口
# 解析输入参数
if isinstance(expr, (int, sp.Integer, float, sp.Float)):
# 单个参数，例如 "64"
m = int(expr)
elif isinstance(expr, tuple):
# 元组参数，例如 "64, 0.4" 或 "64, 0.4, 0"
if len(expr) >= 1:
m = int(expr[0])
if len(expr) >= 2:
alpha = float(expr[1])
if len(expr) >= 3:
# 处理第三个参数（sym_flag）
sym_flag = str(expr[2])
sym = False if sym_flag == '0' else True
if len(expr) > 3:
raise ValueError("参数过多，最多支持3个参数")
else:
return f"输入错误：无法解析的参数格式 '{input_str}'"
# 参数有效性检查
if m is None or m <= 0:
return "错误：窗口长度M必须为正整数"
if alpha <= 0:
return "错误：alpha必须为正数"
# 计算标准差σ，根据MATLAB公式 σ = (M-1)/(2*alpha)
std = (m - 1) / (2 * alpha)
# 生成高斯窗口
window = signal.windows.gaussian(m, std, sym)
return list(window)
except Exception as e:
return f"错误：{e}"
# 示例代码
if __name__ == "__main__":
# 示例1：默认参数（M=64, alpha=2.5, sym=True）
print("示例1 默认参数:", gausswin_window("64")[:5])
# [0.04393693362340742, 0.053411098929097026, 0.06452050267116646, 0.0774512497553297, 0.09238970420181072]
# 示例2：自定义alpha（M=64, alpha=0.4, sym=True）
print("示例2 自定义alpha:", gausswin_window("64,0.4")[:5])
# [0.9231163463866358, 0.9277423175922188, 0.9322411350082428, 0.936610727772961, 0.9408490778060681]
# 示例3：非对称窗口（M=64, alpha=0.4, sym=False）
print("示例3 非对称窗口:", gausswin_window("64,0.4,0")[:5])
# [0.9207563392996183, 0.9254450947656505, 0.9300077511714193, 0.9344422118574613, 0.938746432156527]
最大公约数
G = gcd(A,B,v=0) 返回 A 和 B 的元素的最大公约数. G 中的元素始终是非负值, gcd(0,0) 返回 0. 此语法支持任何数值类型的输入
[G,U,V] = gcd(A,B,v=1) 还返回 Bézout 系数 U 和 V，它们满足：A.*U + B.*V = G。Bézout 系数用于对 Diophantine 方程求解。
A,B — 输入值,标量,向量或实整数值数组
示例1. 分数简化 - 数学计算
fractions = [
(24, 36),
(48, 60),
(125, 75),
(81, 27),
(17, 51)
]
for numerator, denominator in fractions:
# 计算GCD
gcd_result = gcd(numerator, denominator)
# 简化分数
simplified_num = numerator // gcd_result
simplified_den = denominator // gcd_result
print(f"{numerator}/{denominator} = {simplified_num}/{simplified_den} (GCD: {gcd_result})")
#24/36 = 2/3 (GCD: 12)
#48/60 = 4/5 (GCD: 12)
#125/75 = 5/3 (GCD: 25)
#81/27 = 3/1 (GCD: 27)
#17/51 = 1/3 (GCD: 17)
示例2. 时间周期同步 - 物理系统
三个行星的公转周期（天）
planet_periods = [88, 225, 365] # 水星、金星、地球
计算它们同时回到起点的最小时间（最小公倍数）
def lcm(a, b):
return abs(a * b) // gcd(a, b)
def lcm_multiple(numbers):
result = 1
for num in numbers:
result = lcm(result, num)
return result
sync_time = lcm_multiple(planet_periods)
print(f"行星公转周期: {planet_periods} 天")
print(f"它们同时回到起点需要: {sync_time} 天")
print(f"即大约 {sync_time / 365:.1f} 年")
#行星公转周期: [88, 225, 365] 天
#它们同时回到起点需要: 1445400 天
#即大约 3960.0 年
验证GCD计算
for i in range(len(planet_periods)):
for j in range(i + 1, len(planet_periods)):
gcd_val = gcd(planet_periods[i], planet_periods[j])
print(f"GCD({planet_periods[i]}, {planet_periods[j]}) = {gcd_val}")
#GCD(88, 225) = 1
#GCD(88, 365) = 1
#GCD(225, 365) = 5
示例3. 密码学 - RSA算法中的密钥生成
选择两个质数
p, q = 61, 53
计算n和φ(n)
n = p * q
phi_n = (p - 1) * (q - 1)
选择的质数: p = {p}, q = {q}
print(f"n = p × q = {n}")
print(f"φ(n) = (p-1)(q-1) = {phi_n}")
#n = p × q = 3233
#φ(n) = (p-1)(q-1) = 3120
选择加密指数e，需要满足 gcd(e, φ(n)) = 1
possible_e_values = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31]
寻找合适的加密指数e:
for e in possible_e_values:
gcd_val = gcd(e, phi_n)
status = "✓ 合适" if gcd_val == 1 else "✗ 不合适"
print(f"e = {e}: gcd({e}, {phi_n}) = {gcd_val} {status}")
#e = 3: gcd(3, 3120) = 3 ✗ 不合适
#e = 5: gcd(5, 3120) = 5 ✗ 不合适
#e = 7: gcd(7, 3120) = 1 ✓ 合适
#e = 11: gcd(11, 3120) = 1 ✓ 合适
#e = 13: gcd(13, 3120) = 13 ✗ 不合适
#e = 17: gcd(17, 3120) = 1 ✓ 合适
#e = 19: gcd(19, 3120) = 1 ✓ 合适
#e = 23: gcd(23, 3120) = 1 ✓ 合适
#e = 29: gcd(29, 3120) = 1 ✓ 合适
#e = 31: gcd(31, 3120) = 1 ✓ 合适
选择第一个合适的e
e = next(e for e in possible_e_values if gcd(e, phi_n) == 1)
print(f"\n选择的加密指数: e = {e}")
#选择的加密指数: e = 7
使用扩展欧几里得算法计算解密指数d
gcd_result = gcd(e, phi_n, 1)
g, d, _ = gcd_result
# 确保d是正数
d = d[0] if d[0] > 0 else d[0] + phi_n
print(f"计算的解密指数: d = {d}")
print(f"验证: e × d mod φ(n) = {e} × {d} mod {phi_n} = {(e * d) % phi_n}")
#计算的解密指数: d = 1783
#验证: e × d mod φ(n) = 7 × 1783 mod 3120 = 1
示例4. 工程测量 - 齿轮传动比
齿轮系统的设计需求
required_ratio_numerator = 120
required_ratio_denominator = 75
计算GCD来简化传动比
gcd_val = gcd(required_ratio_numerator, required_ratio_denominator)
simplified_numerator = required_ratio_numerator // gcd_val
simplified_denominator = required_ratio_denominator // gcd_val
print(f"原始传动比: {required_ratio_numerator}:{required_ratio_denominator}")
print(f"最大公约数: {gcd_val}")
print(f"简化传动比: {simplified_numerator}:{simplified_denominator}")
#原始传动比: 120:75
#最大公约数: 15
#简化传动比: 8:5
实际齿轮齿数选择（考虑制造限制）
min_teeth = 15
max_teeth = 120
寻找满足比例且在限制范围内的齿数
multiplier = 1
while True:
gear1 = simplified_numerator * multiplier
gear2 = simplified_denominator * multiplier
if gear1 <= max_teeth and gear2 <= max_teeth and gear1 >= min_teeth and gear2 >= min_teeth:
break
multiplier += 1
推荐的齿轮齿数:
print(f"主动轮: {gear1} 齿")
print(f"从动轮: {gear2} 齿")
print(f"实际传动比: {gear1}:{gear2} = {gear1 / gear2:.3f}")
#主动轮: 24 齿
#从动轮: 15 齿
#实际传动比: 24:15 = 1.600
验证GCD
verification_gcd = gcd_number(f"{gear1}, {gear2}")
print(f"验证GCD({gear1}, {gear2}) = {verification_gcd}")
#验证GCD(24, 15) = 3
示例5. 计算机图形学 - 像素比例简化
常见屏幕分辨率
resolutions = [
(1920, 1080),
(2560, 1440),
(3840, 2160),
(2560, 1080),
(3440, 1440),
(1920, 1200),
(1680, 1050),
]
屏幕分辨率比例分析:
for width, height in resolutions:
gcd_val = gcd(width, height)
aspect_width = width // gcd_val
aspect_height = height // gcd_val
print(f"{width}×{height} → 比例 {aspect_width}:{aspect_height} (GCD: {gcd_val})")
#1920×1080 → 比例 16:9 (GCD: 120)
#2560×1440 → 比例 16:9 (GCD: 160)
#3840×2160 → 比例 16:9 (GCD: 240)
#2560×1080 → 比例 64:27 (GCD: 40)
#3440×1440 → 比例 43:18 (GCD: 80)
#1920×1200 → 比例 8:5 (GCD: 240)
#1680×1050 → 比例 8:5 (GCD: 210)
图像缩放时的比例保持
original_width, original_height = 1200, 800
target_width = 600
gcd_val = gcd(original_width, original_height)
aspect_ratio = original_width / original_height
计算保持比例的目标高度
target_height = int(target_width / aspect_ratio)
print(f"原始尺寸: {original_width}×{original_height} (比例 {original_width // gcd_val}:{original_height // gcd_val})")
print(f"缩放宽度到: {target_width}")
print(f"保持比例的高度: {target_height}")
#原始尺寸: 1200×800 (比例 3:2)
#缩放宽度到: 600
#保持比例的高度: 400
验证新尺寸的GCD
new_gcd = gcd(target_width, target_height)
print(f"新尺寸比例: {target_width // new_gcd}:{target_height // new_gcd}")
#新尺寸比例: 3:2
示例6. 音乐理论 - 节拍和节奏
不同音符的时值（以最小时间单位为1）
note_durations = {
"全音符": 16,
"二分音符": 8,
"四分音符": 4,
"八分音符": 2,
"十六分音符": 1
}
计算各种音符组合的GCD来确定最小时间单位
durations = list(note_durations.values())
计算所有时值的GCD
current_gcd = durations[0]
for duration in durations[1:]:
current_gcd = gcd(current_gcd, duration)
音符时值系统:
for note, duration in note_durations.items():
relative_duration = duration // current_gcd
print(f"{note}: {duration} → 相对时值 {relative_duration}")
#全音符: 16 → 相对时值 16
#二分音符: 8 → 相对时值 8
#四分音符: 4 → 相对时值 4
#八分音符: 2 → 相对时值 2
#十六分音符: 1 → 相对时值 1
print(f"\n最小时间单位: {current_gcd}")
#最小时间单位: 1
节奏模式分析
rhythm_patterns = [
[4, 4, 4, 4], # 4/4拍
[4, 4, 2], # 3/4拍
[8, 8, 4, 8, 8], # 复合节奏
[2, 2, 2, 2, 2, 2], # 6/8拍
]
节奏模式分析:
for i, pattern in enumerate(rhythm_patterns):
pattern_gcd = pattern[0]
for duration in pattern[1:]:
pattern_gcd = gcd(pattern_gcd, duration)
simplified_pattern = [d // pattern_gcd for d in pattern]
print(f"模式 {i + 1}: {pattern} → 简化 {simplified_pattern} (GCD: {pattern_gcd})")
#模式 1: [4, 4, 4, 4] → 简化 [1, 1, 1, 1] (GCD: 4)
#模式 2: [4, 4, 2] → 简化 [2, 2, 1] (GCD: 2)
#模式 3: [8, 8, 4, 8, 8] → 简化 [2, 2, 1, 2, 2] (GCD: 4)
#模式 4: [2, 2, 2, 2, 2, 2] → 简化 [1, 1, 1, 1, 1, 1] (GCD: 2)
示例7. 资源分配 - 最优包装方案
仓库包装问题：有不同尺寸的箱子，要装满固定容量的货柜
box_sizes = [12, 18, 24] # 三种箱子的容量
container_capacity = 360 # 货柜总容量
计算所有箱子尺寸的GCD
gcd_all = box_sizes[0]
for size in box_sizes[1:]:
gcd_all = gcd(gcd_all, size)
print(f"箱子尺寸: {box_sizes}")
print(f"所有箱子尺寸的GCD: {gcd_all}")
print(f"货柜容量: {container_capacity}")
#箱子尺寸: [12, 18, 24]
#所有箱子尺寸的GCD: 6
#货柜容量: 360
检查货柜容量是否能被GCD整除
if container_capacity % gcd_all == 0:
print("✓ 货柜可以被完全装满")
else:
print("✗ 货柜无法被完全装满")
print(f"最大可利用容量: {container_capacity - (container_capacity % gcd_all)}")
#✓ 货柜可以被完全装满
寻找最优包装方案
solutions = []
for count1 in range(container_capacity // box_sizes[0] + 1):
for count2 in range(container_capacity // box_sizes[1] + 1):
for count3 in range(container_capacity // box_sizes[2] + 1):
total = count1 * box_sizes[0] + count2 * box_sizes[1] + count3 * box_sizes[2]
if total == container_capacity:
solutions.append((count1, count2, count3))
if solutions:
print(f"找到 {len(solutions)} 种包装方案")
for i, sol in enumerate(solutions[:3]): # 只显示前3种方案
print(f"方案 {i + 1}: {sol[0]}个{box_sizes[0]}箱 + {sol[1]}个{box_sizes[1]}箱 + {sol[2]}个{box_sizes[2]}箱")
else:
print("没有找到完全装满的方案")
#找到 91 种包装方案
#方案 1: 0个12箱 + 0个18箱 + 15个24箱
#方案 2: 0个12箱 + 4个18箱 + 12个24箱
#方案 3: 0个12箱 + 8个18箱 + 9个24箱
示例8. 数论问题 - 扩展欧几里得算法应用
线性同余方程求解：ax ≡ b (mod m)
equations = [
(7, 3, 11), # 7x ≡ 3 (mod 11)
(15, 6, 21), # 15x ≡ 6 (mod 21)
(8, 5, 12), # 8x ≡ 5 (mod 12)
]
for a, b, m in equations:
print(f"\n求解方程: {a}x ≡ {b} (mod {m})")
# 计算 gcd(a, m)
gcd_am = gcd(a, m)
print(f"gcd({a}, {m}) = {gcd_am}")
# 检查解的存在性
if b % gcd_am != 0:
print("无解，因为 gcd(a,m) 不能整除 b")
continue
# 使用扩展欧几里得算法求特解
result = gcd(f"{a}, {m}, 1")
if isinstance(result, tuple):
g, u, v = result
u_val, v_val = u[0], v[0]
print(f"扩展欧几里得结果: {a}×({u_val}) + {m}×({v_val}) = {g}")
# 计算特解
x0 = (u_val * (b // g)) % m
if x0 < 0:
x0 += m
print(f"特解: x₀ = {x0}")
# 所有解
if gcd_am > 1:
print(f"所有解: x ≡ {x0} + k×{m // gcd_am} (mod {m}), k=0,1,...,{gcd_am - 1}")
solutions = [(x0 + k * (m // gcd_am)) % m for k in range(gcd_am)]
print(f"具体解: {solutions}")
else:
print(f"唯一解: x ≡ {x0} (mod {m})")
# 验证
verification = (a * x0) % m
print(f"验证: {a}×{x0} mod {m} = {verification} ≡ {b} (mod {m})")
#求解方程: 7x ≡ 3 (mod 11)
#gcd(7, 11) = 1
#扩展欧几里得结果: 7×(-3) + 11×(2) = [1]
#特解: x₀ = [2]
#唯一解: x ≡ [2] (mod 11)
#验证: 7×[2] mod 11 = [3] ≡ 3 (mod 11)
#求解方程: 15x ≡ 6 (mod 21)
#gcd(15, 21) = 3
#扩展欧几里得结果: 15×(3) + 21×(-2) = [3]
#特解: x₀ = [6]
#所有解: x ≡ [6] + k×7 (mod 21), k=0,1,...,2
#具体解: [array([6]), array([13]), array([20])]
#验证: 15×[6] mod 21 = [6] ≡ 6 (mod 21)
#求解方程: 8x ≡ 5 (mod 12)
#gcd(8, 12) = 4
#无解，因为 gcd(a,m) 不能整除 b
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def gcd_number(input_str):
"""
对标MATLAB的gcd函数，计算最大公约数及Bézout系数
参数：
input_str (str): 输入的参数字符串，格式应为 "a, b"（支持标量或矩阵）
返回：
tuple: (G, U, V) 包含最大公约数和Bézout系数的矩阵/标量
str: 错误信息（格式错误/非整数输入/矩阵形状不匹配等）
功能说明：
1. 支持标量计算：gcd(6,15) → (3, -2, 1)
2. 支持矩阵计算：逐元素计算gcd
3. 自动处理标量与矩阵的广播
4. 返回Bézout系数满足 a*u + b*v = gcd(a,b)
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def extended_gcd(a, b):
# 转换为 NumPy 数组并确保是整数类型
a0 = np.array(a, dtype=np.int64, ndmin=1, copy=False)
b0 = np.array(b, dtype=np.int64, ndmin=1, copy=False)
# 应用广播规则
a_broadcast, b_broadcast = np.broadcast_arrays(a0, b0)
# 初始化结果数组
shape = a_broadcast.shape
g_arr = np.zeros(shape, dtype=np.int64)
u_arr = np.zeros(shape, dtype=np.int64)
v_arr = np.zeros(shape, dtype=np.int64)
# 对每个元素单独计算扩展欧几里得算法
for idx in np.ndindex(shape):
a_val = a_broadcast[idx]
b_val = b_broadcast[idx]
sign_a = np.sign(a_val) if a_val != 0 else 1
sign_b = np.sign(b_val) if b_val != 0 else 1
a1 = np.abs(a_val)
b1 = np.abs(b_val)
u, v = 1, 0
u1, v1 = 0, 1
while b1 != 0:
q = a1 // b1
r = a1 % b1
new_u = u - q * u1
new_v = v - q * v1
a1, b1 = b1, r
u, v = u1, v1
u1, v1 = new_u, new_v
g_arr[idx] = a1
u_arr[idx] = u * sign_a
v_arr[idx] = v * sign_b
return g_arr, u_arr, v_arr
# 检查传入的 expr 是否为元组类型
if isinstance(expr, tuple):
# 如果 expr 是一个长度为 3 的元组
if len(expr) == 3:
# 从元组中提取第一个元素作为底数 n（这里注释提到 base b 可能有误，推测应该是 base n）
n = expr[0]
# 从元组中提取第二个元素
x = expr[1]
# 将元组中的第三个元素转换为整数类型
v = int(expr[2])
# 检查 v 是否既不等于 0 也不等于 1
if v != 0 and v != 1:
# 如果不满足条件，抛出 ValueError 异常，提示参数必须是 0 或 1
raise ValueError("参数必须是0或1")
# 如果 expr 是一个长度为 2 的元组
elif len(expr) == 2:
# 从元组中提取第一个元素作为底数 n
n = expr[0]
# 从元组中提取第二个元素
x = expr[1]
# 将 v 的值设为 0
v = 0
if v == 0:
result = np.gcd(n, x)
else:
result = extended_gcd(n, x)
# 如果 expr 不是元组类型
else:
# 将错误标志设为 True
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"输入错误:{input_str}"
# 示例代码 ======================================================================
if __name__ == "__main__":
# 示例1：标量计算
print("示例1：gcd(6, 15, 1)")
result = gcd_number("6, 15, 1")
print(result)
print("\n")
# (array([3]), array([-2]), array([1]))
# 示例2：
test_input = "12,15"
print("输入：", test_input)
result = gcd_number(test_input)
print(result)
print("\n")
# 3
test_input = "5, 15, 1"
print("输入：", test_input)
result = gcd_number(test_input)
print(result)
print("\n")
# (array([5]), array([1]), array([0]))
最大公约数矩阵
A = gcdmat(n) 返回 n×n 矩阵，其 A(i,j) 等于 gcd(i,j)
对于所有非负r, 矩阵 A 是对称正定矩阵, A^r 是对称半正定矩阵
n — 输入,标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def generate_gcd_matrix(n):
"""
生成n×n的GCD矩阵，其中元素A[i][j] = gcd(i+1, j+1)
参数:
n (int): 矩阵维度（正整数）
返回:
np.ndarray: 生成的GCD矩阵
"""
# 初始化n×n的零矩阵
A = np.zeros((n, n), dtype=int)
# 遍历每个元素计算gcd（注意i和j从0开始，对应1-based编号）
for i in range(n):
for j in range(n):
A[i][j] = sp.gcd(i + 1, j + 1)
return A
def gcdmat_matrix(input_str):
"""
根据输入字符串生成GCD矩阵
参数:
input_str (str): 输入的矩阵维度字符串表达式（需能解析为正整数）
返回:
sp.Matrix | str: 生成的Sympy矩阵对象或错误信息
"""
try:
# 将输入字符串解析为Sympy表达式
expr = sp.sympify(input_str)
# 检查是否为有效正整数
if not expr.is_integer: # 检查数学上是否为整数
return f"错误：输入必须为整数，但收到 {input_str}"
n = int(expr) # 转换为Python整数
if n <= 0:
return f"错误：维度必须为正整数，但收到 {n}"
# 生成矩阵并转换为Sympy矩阵对象
return sp.Matrix(generate_gcd_matrix(n))
except sp.SympifyError:
return f"解析错误：'{input_str}' 不是有效的数学表达式"
except Exception as e:
return f"运行时错误：{str(e)}"
if __name__ == "__main__":
# 示例测试代码
test_cases = [
"3",
# Matrix([[1, 1, 1],
# [1, 2, 1],
# [1, 1, 3]])
"4"
# Matrix([[1, 1, 1, 1],
# [1, 2, 1, 2],
# [1, 1, 3, 1],
# [1, 2, 1, 4]])
]
for case in test_cases:
print(f"输入 '{case}':")
result = gcdmat_matrix(case)
if isinstance(result, sp.Matrix):
print(f"生成的 {result.rows}x{result.cols} 矩阵:")
print(result)
else:
print(result)
print("\n" + "-" * 50 + "\n")
Gear矩阵
A = gearmat(n,i,j) 返回 n×n 矩阵, 其中下对角线和上对角线上为 1, (1,abs(i)) 位置为 sign(i), (n,n+1-abs(j)) 位置为 sign(j), 其余所有位置为 0.
参量 i 和 j 的默认值分别为 n 和 -n. 它们必须为 -n 到 n 范围内的整数
矩阵 A 是奇异矩阵，可具有双重和三重特征值，并且可以是亏损矩阵
所有特征值的形式均为 2*cos(a)，特征向量的形式为 [sin(w+a), sin(w+2a), …, sin(w+na)]
n,i,j — 输入,标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def generate_gear_matrix(n, i, j):
"""
生成n×n的齿轮矩阵，具有以下特征：
1. 次对角线（上下对角线）填充1
2. (1,|i|)位置填充sign(i)
3. (n, n+1-|j|)位置填充sign(j)
参数:
n (int): 矩阵维度（>=2的正整数）
i (int): 第一行特征位置（非零整数，-n <= i <=n）
j (int): 最后一行特征位置（非零整数，-n <= j <=n）
返回:
np.ndarray: 生成的齿轮矩阵
异常:
ValueError: 参数不符合要求时抛出
"""
# 参数验证
if n < 1:
raise ValueError("矩阵维度n必须是正整数")
if i == 0 or j == 0:
raise ValueError("参数i和j不能为零")
if not (-n <= i <= n) or not (-n <= j <= n):
raise ValueError(f"参数i和j必须在[-{n}, {n}]范围内")
# 创建基础矩阵
A = np.zeros((n, n), dtype=int)
# 填充次对角线（Python使用0-based索引）
# 下对角线（主对角线下方）
np.fill_diagonal(A[1:], 1)
# 上对角线（主对角线上方）
np.fill_diagonal(A[:, 1:], 1)
# 设置特征元素（转换为0-based索引）
# 第一行：位置(0, abs(i)-1)
A[0, abs(i) - 1] = np.sign(i)
# 最后一行：位置(n-1, n - abs(j) -1)
# 计算公式调整为：n - abs(j) 对应MATLAB的n+1-|j|
A[n - 1, n - abs(j)] = np.sign(j)
return A
def gearmat_matrix(input_str):
"""
解析输入字符串并生成齿轮矩阵
参数:
input_str (str): 支持两种格式：
- 单整数："n" 生成默认矩阵（i=n, j=-n）
- 元组："(n,i,j)" 或 "n,i,j" 生成自定义矩阵
返回:
sp.Matrix | str: Sympy矩阵对象或错误信息
"""
try:
expr = sp.sympify(input_str)
# 处理元组输入 (n, i, j)
if isinstance(expr, tuple):
if len(expr) != 3:
return f"需要3个参数，但收到{len(expr)}个"
# 验证所有元素为整数
if not all(e.is_integer for e in expr):
return "所有参数必须为整数"
n, i, j = map(int, expr)
matrix = generate_gear_matrix(n, i, j)
# 处理单整数输入
elif expr.is_integer:
n = int(expr)
if n < 1:
return f"维度必须≥1，但收到{n}"
matrix = generate_gear_matrix(n, i=n, j=-n) # 默认参数
else:
return "无效的输入格式"
return sp.Matrix(matrix)
except (sp.SympifyError, ValueError) as e:
return f"输入错误：{str(e)}"
except Exception as e:
return f"运行时错误：{str(e)}"
if __name__ == "__main__":
# 示例测试
test_cases = [
"3",
# Matrix([[0, 1, 1],
# [1, 0, 1],
# [-1, 1, 0]])
"(4, 2, -3)",
# Matrix([[0, 1, 0, 0],
# [1, 0, 1, 0],
# [0, 1, 0, 1],
# [0, -1, 1, 0]])
]
for case in test_cases:
print(f"输入: {case}")
result = gearmat_matrix(case)
if isinstance(result, sp.Matrix):
print(f"生成 {result.rows}x{result.cols} 齿轮矩阵:")
print(result)
else:
print(result)
print("\n" + "=" * 50 + "\n")
几何分布的累积分布函数
y=geocdf（x，p）返回几何分布的累积分布函数（cdf），使用p中的相应概率在x中的每个值处进行评估。
x -- 整数标量, 向量，矩阵，多维数组，是第一次成功之前的失败次数
p -- 标量，向量，矩阵，多维数组，是成功的概率
y -- cdf值，作为标量或[0,1]范围内的标量数组返回。在任何必要的标量展开后，y与x和p的大小相同。
示例1. 质量控制 - 检测第一个次品
假设生产线的次品率为1%
defect_rate = 0.01
print(f"生产线次品率: {defect_rate * 100}%")
#生产线次品率: 1.0%
在前N个产品中发现第一个次品的概率:
for n in [10, 50, 100, 200, 500]:
prob = geocdf(n, defect_rate)
print(f"前{n}个产品中发现第一个次品的概率: {prob:.4f} ({prob * 100:.2f}%)")
#前10个产品中发现第一个次品的概率: 0.0956 (9.56%)
#前50个产品中发现第一个次品的概率: 0.3950 (39.50%)
#前100个产品中发现第一个次品的概率: 0.6340 (63.40%)
#前200个产品中发现第一个次品的概率: 0.8660 (86.60%)
#前500个产品中发现第一个次品的概率: 0.9934 (99.34%)
计算达到特定置信度需要的检查数量
confidence_levels = [0.5, 0.9, 0.95, 0.99]
达到特定置信度需要检查的产品数量:
for confidence in confidence_levels:
# 解方程: 1 - (1-p)^n = confidence
n_needed = np.ceil(np.log(1 - confidence) / np.log(1 - defect_rate))
actual_prob = geocdf(n_needed, defect_rate)
print(f"{confidence * 100}% 置信度: 需要检查 {int(n_needed)} 个产品 (实际概率: {actual_prob:.4f})")
#50.0% 置信度: 需要检查 69 个产品 (实际概率: 0.5002)
#90.0% 置信度: 需要检查 230 个产品 (实际概率: 0.9009)
#95.0% 置信度: 需要检查 299 个产品 (实际概率: 0.9505)
#99.0% 置信度: 需要检查 459 个产品 (实际概率: 0.9901)
示例2. 网络传输 - 数据包重传分析
假设网络丢包率为5%
packet_loss_rate = 0.05
success_rate = 1 - packet_loss_rate
print(f"网络丢包率: {packet_loss_rate * 100}%")
print(f"单次传输成功率: {success_rate * 100}%")
#网络丢包率: 5.0%
#单次传输成功率: 95.0%
分析不同重传次数下的成功概率
max_attempts = 10
for attempts in range(1, max_attempts + 1):
prob = geocdf(attempts, success_rate)
print(f"前{attempts}次尝试内传输成功的概率: {prob:.4f} ({prob * 100:.2f}%)")
#前1次尝试内传输成功的概率: 0.9500 (95.00%)
#前2次尝试内传输成功的概率: 0.9975 (99.75%)
#前3次尝试内传输成功的概率: 0.9999 (99.99%)
#前4次尝试内传输成功的概率: 1.0000 (100.00%)
#前5次尝试内传输成功的概率: 1.0000 (100.00%)
#前6次尝试内传输成功的概率: 1.0000 (100.00%)
#前7次尝试内传输成功的概率: 1.0000 (100.00%)
#前8次尝试内传输成功的概率: 1.0000 (100.00%)
#前9次尝试内传输成功的概率: 1.0000 (100.00%)
#前10次尝试内传输成功的概率: 1.0000 (100.00%)
计算达到99%成功概率所需的最大重传次数
target_prob = 0.99
max_retries = np.ceil(np.log(1 - target_prob) / np.log(packet_loss_rate))
actual_prob = geocdf(max_retries, success_rate)
print(f"\n达到{target_prob * 100}%成功概率需要最多 {int(max_retries)} 次重传")
print(f"实际概率: {actual_prob:.4f}")
#达到99.0%成功概率需要最多 2 次重传
#实际概率: 0.9975
不同网络质量下的比较
network_conditions = {
"优秀网络 (丢包率1%)": 0.01,
"良好网络 (丢包率5%)": 0.05,
"一般网络 (丢包率10%)": 0.10,
"较差网络 (丢包率20%)": 0.20
}
for desc, loss_rate in network_conditions.items():
success_prob = geocdf(5, 1 - loss_rate) # 5次尝试内的成功概率
print(f"{desc}: 5次尝试内成功概率 = {success_prob:.4f}")
#优秀网络 (丢包率1%): 5次尝试内成功概率 = 1.0000
#良好网络 (丢包率5%): 5次尝试内成功概率 = 1.0000
#一般网络 (丢包率10%): 5次尝试内成功概率 = 1.0000
#较差网络 (丢包率20%): 5次尝试内成功概率 = 0.9997
示例3. 市场营销 - 客户转化分析
假设每次营销接触的转化率为2%
conversion_rate = 0.02
print(f"单次营销接触转化率: {conversion_rate * 100}%")
#单次营销接触转化率: 2.0%
分析不同接触次数下的累积转化概率
contact_points = [1, 3, 5, 7, 10, 15, 20]
不同营销接触次数下的累积转化概率:
for contacts in contact_points:
prob = geocdf(contacts, conversion_rate)
print(f"前{contacts:2d}次接触内转化的概率: {prob:.4f} ({prob * 100:6.2f}%)")
#前 1次接触内转化的概率: 0.0200 ( 2.00%)
#前 3次接触内转化的概率: 0.0588 ( 5.88%)
#前 5次接触内转化的概率: 0.0961 ( 9.61%)
#前 7次接触内转化的概率: 0.1319 ( 13.19%)
#前10次接触内转化的概率: 0.1829 ( 18.29%)
#前15次接触内转化的概率: 0.2614 ( 26.14%)
#前20次接触内转化的概率: 0.3324 ( 33.24%)
计算最优营销策略
成本假设：每次接触成本1单位，转化价值100单位
contact_cost = 1
conversion_value = 100
optimal_contacts = 0
max_profit = 0
for contacts in range(1, 31):
conversion_prob = geocdf(contacts, conversion_rate)
expected_value = conversion_prob * conversion_value
total_cost = contacts * contact_cost
expected_profit = expected_value - total_cost
if expected_profit > max_profit:
max_profit = expected_profit
optimal_contacts = contacts
if contacts in [5, 10, 15, 20, 25, 30]:
print(f"{contacts:2d}次接触: 期望利润 = {expected_profit:6.2f}")
# 5次接触: 期望利润 = 4.61
#10次接触: 期望利润 = 8.29
#15次接触: 期望利润 = 11.14
#20次接触: 期望利润 = 13.24
#25次接触: 期望利润 = 14.65
#30次接触: 期望利润 = 15.45
print(f"\n最优策略: 进行 {optimal_contacts} 次营销接触")
print(f"最大期望利润: {max_profit:.2f}")
#最优策略: 进行 30 次营销接触
#最大期望利润: 15.45
示例4. 医学研究 - 疾病检测灵敏度
假设某种检测方法的灵敏度为90%
test_sensitivity = 0.90
print(f"单次检测灵敏度: {test_sensitivity * 100}%")
#单次检测灵敏度: 90.0% (即对真实患者检测出阳性的概率)
分析重复检测对灵敏度的提升
for num_tests in [1, 2, 3, 4, 5]:
# 至少一次检测为阳性的概率
overall_sensitivity = geocdf(num_tests, test_sensitivity)
print(f"{num_tests}次检测的总体灵敏度: {overall_sensitivity:.6f} ({overall_sensitivity * 100:.4f}%)")
#1次检测的总体灵敏度: 0.900000 (90.0000%)
#2次检测的总体灵敏度: 0.990000 (99.0000%)
#3次检测的总体灵敏度: 0.999000 (99.9000%)
#4次检测的总体灵敏度: 0.999900 (99.9900%)
#5次检测的总体灵敏度: 0.999990 (99.9990%)
不同疾病检测场景
scenarios = {
"高灵敏度检测 (95%)": 0.95,
"标准检测 (90%)": 0.90,
"快速检测 (85%)": 0.85,
"筛查检测 (80%)": 0.80
}
for desc, sensitivity in scenarios.items():
# 3次检测内的累积灵敏度
cumul_sensitivity = geocdf(3, sensitivity)
print(f"{desc}: 3次检测累积灵敏度 = {cumul_sensitivity:.6f}")
#高灵敏度检测 (95%): 3次检测累积灵敏度 = 0.999875
#标准检测 (90%): 3次检测累积灵敏度 = 0.999000
#快速检测 (85%): 3次检测累积灵敏度 = 0.996625
#筛查检测 (80%): 3次检测累积灵敏度 = 0.992000
计算达到99.9%灵敏度所需的检测次数
target_sensitivity = 0.999
tests_needed = np.ceil(np.log(1 - target_sensitivity) / np.log(1 - test_sensitivity))
actual_sensitivity = geocdf(f"{tests_needed}, {test_sensitivity}")
print(f"\n达到{target_sensitivity * 100}%灵敏度需要 {int(tests_needed)} 次检测")
print(f"实际灵敏度: {actual_sensitivity:.6f}")
#达到99.9%灵敏度需要 3 次检测
#实际灵敏度: 0.999000
示例5. 游戏设计 - 稀有物品掉落机制
假设稀有物品的掉落率为1%
drop_rate = 0.01
print(f"稀有物品掉落率: {drop_rate * 100}%")
#稀有物品掉落率: 1.0%
分析玩家在不同尝试次数下获得物品的概率
attempts_range = [10, 50, 100, 200, 500, 1000]
在不同尝试次数下获得至少一个稀有物品的概率:
for attempts in attempts_range:
prob = geocdf(attempts, drop_rate)
print(f"{attempts:4d} 次尝试: {prob:.4f} ({prob * 100:6.2f}%)")
# 10 次尝试: 0.0956 ( 9.56%)
# 50 次尝试: 0.3950 ( 39.50%)
# 100 次尝试: 0.6340 ( 63.40%)
# 200 次尝试: 0.8660 ( 86.60%)
# 500 次尝试: 0.9934 ( 99.34%)
#1000 次尝试: 1.0000 (100.00%)
计算"保底"机制的设计
pity_system_thresholds = [100, 200, 300, 500]
for threshold in pity_system_thresholds:
natural_prob = geocdf(threshold, drop_rate)
print(f"在{threshold}次尝试内自然获得的概率: {natural_prob:.4f}")
print(f" 需要保底机制覆盖的玩家比例: {(1 - natural_prob) * 100:.2f}%")
#在100次尝试内自然获得的概率: 0.6340
# 需要保底机制覆盖的玩家比例: 36.60%
#在200次尝试内自然获得的概率: 0.8660
# 需要保底机制覆盖的玩家比例: 13.40%
#在300次尝试内自然获得的概率: 0.9510
# 需要保底机制覆盖的玩家比例: 4.90%
#在500次尝试内自然获得的概率: 0.9934
# 需要保底机制覆盖的玩家比例: 0.66%
玩家体验分析
milestone_attempts = [50, 100, 150, 200]
for attempts in milestone_attempts:
prob = geocdf(attempts, drop_rate)
unlucky_players = (1 - prob) * 100
print(f"尝试{attempts}次仍未获得的玩家: {unlucky_players:.2f}%")
#尝试50次仍未获得的玩家: 60.50%
#尝试100次仍未获得的玩家: 36.60%
#尝试150次仍未获得的玩家: 22.15%
#尝试200次仍未获得的玩家: 13.40%
示例6. 可靠性工程 - 设备故障分析
假设设备每天发生故障的概率为0.5%
daily_failure_rate = 0.005
print(f"设备每日故障率: {daily_failure_rate * 100}%")
print(f"平均无故障时间: {1 / daily_failure_rate:.1f} 天")
#设备每日故障率: 0.5%
#平均无故障时间: 200.0 天
分析不同时间范围内的故障概率
time_periods = [30, 90, 180, 365, 730] # 天数
在不同时间范围内发生首次故障的概率:
for days in time_periods:
prob = geocdf(days, daily_failure_rate)
print(f"{days:3d} 天内发生故障的概率: {prob:.4f} ({prob * 100:.2f}%)")
# 30 天内发生故障的概率: 0.1396 (13.96%)
# 90 天内发生故障的概率: 0.3631 (36.31%)
#180 天内发生故障的概率: 0.5943 (59.43%)
#365 天内发生故障的概率: 0.8395 (83.95%)
#730 天内发生故障的概率: 0.9742 (97.42%)
不同可靠性要求的维护策略
reliability_targets = {
"关键设备 (99.9%)": 0.999,
"重要设备 (99%)": 0.99,
"一般设备 (95%)": 0.95
}
for desc, target in reliability_targets.items():
# 计算达到可靠性目标的最大运行天数
max_days = np.floor(np.log(1 - target) / np.log(1 - daily_failure_rate))
actual_reliability = 1 - geocdf(max_days, daily_failure_rate)
print(f"{desc}: 最多运行 {int(max_days)} 天 (实际可靠性: {actual_reliability:.6f})")
#关键设备 (99.9%): 最多运行 1378 天 (实际可靠性: 0.001000)
#重要设备 (99%): 最多运行 918 天 (实际可靠性: 0.010037)
#一般设备 (95%): 最多运行 597 天 (实际可靠性: 0.050163)
不同质量设备的比较
failure_rates = {
"高质量设备": 0.001,
"标准设备": 0.005,
"经济型设备": 0.01
}
for desc, rate in failure_rates.items():
yearly_failure_prob = geocdf(365, rate)
print(f"{desc}: 一年内故障概率 = {yearly_failure_prob:.4f}")
#高质量设备: 一年内故障概率 = 0.3059
#标准设备: 一年内故障概率 = 0.8395
#经济型设备: 一年内故障概率 = 0.9745
示例7. 金融服务 - 贷款违约预测
假设月度违约率为0.2%
monthly_default_rate = 0.002
print(f"贷款月度违约率: {monthly_default_rate * 100}%")
#贷款月度违约率: 0.2%
分析不同贷款期限内的累计违约概率
loan_terms = [12, 24, 36, 48, 60] # 月数
不同贷款期限内的累计违约概率:
for months in loan_terms:
prob = geocdf(months, monthly_default_rate)
print(f"{months:2d} 个月内的违约概率: {prob:.4f} ({prob * 100:.3f}%)")
#12 个月内的违约概率: 0.0237 (2.374%)
#24 个月内的违约概率: 0.0469 (4.691%)
#36 个月内的违约概率: 0.0695 (6.954%)
#48 个月内的违约概率: 0.0916 (9.162%)
#60 个月内的违约概率: 0.1132 (11.319%)
基于风险的贷款定价分析
risk_categories = {
"优质客户 (违约率0.1%)": 0.001,
"标准客户 (违约率0.2%)": 0.002,
"次级客户 (违约率0.5%)": 0.005,
"高风险客户 (违约率1%)": 0.01
}
loan_term = 36 # 3年贷款
base_interest_rate = 0.05 # 基准利率5%
for desc, default_rate in risk_categories.items():
default_prob = geocdf(loan_term, default_rate)
# 简化风险溢价计算
risk_premium = default_prob * 2 # 假设风险溢价为违约概率的2倍
total_interest = base_interest_rate + risk_premium
print(f"{desc}: {loan_term}个月违约概率 = {default_prob:.4f}")
print(f" 建议利率: {(total_interest * 100):.2f}%")
#优质客户 (违约率0.1%): 36个月违约概率 = 0.0354
# 建议利率: 12.08%
#标准客户 (违约率0.2%): 36个月违约概率 = 0.0695
# 建议利率: 18.91%
#次级客户 (违约率0.5%): 36个月违约概率 = 0.1651
# 建议利率: 38.02%
#高风险客户 (违约率1%): 36个月违约概率 = 0.3036
# 建议利率: 65.72%
示例8. 符号计算应用 - 理论分析
使用符号变量进行理论推导
计算几何分布的CDF符号表达式
symbolic_result = geocdf(k, p)
print(f"几何分布CDF的符号表达式: {symbolic_result}")
#几何分布CDF的符号表达式:
#Piecewise((p*Piecewise((floor(k)/(1.0 - p), Eq(p, 0)),
((-p - (1.0 - p)**(floor(k) + 1.0) + 1.0)/(p*(1.0 - p)), True)), k >= 1), (0, True))
分析不同参数下的行为
specific_cases = [
("p = 0.5", geocdf(k, 0.5)),
("p = 0.1", geocdf(k, 0.1)),
("p = 0.01", geocdf(k, 0.01))
]
for desc, expr in specific_cases:
print(f"{desc}: {expr}")
#p = 0.5: Piecewise((1.0 - 2.0*0.5**(floor(k) + 1.0), k >= 1), (0, True))
#p = 0.1: Piecewise((1.0 - 1.11111111111111*0.9**(floor(k) + 1.0), k >= 1), (0, True))
#p = 0.01: Piecewise((0.999999999999999 - 1.01010101010101*0.99**(floor(k) + 1.0), k >= 1), (0, True))
期望值和方差的理论计算
期望值 E[X] = 1/p
expected_value = 1 / p
方差 Var[X] = (1-p)/p²
variance = (1 - p) / (p ** 2)
print(f"期望值 E[X] = {expected_value}")
print(f"方差 Var[X] = {variance}")
#期望值 E[X] = 1/p
#方差 Var[X] = (1 - p)/p**2
可视化不同参数下的CDF形状
try:
import matplotlib.pyplot as plt
# 数值计算示例
k_values = np.arange(0, 50)
p_values = [0.1, 0.3, 0.5, 0.7, 0.9]
plt.figure(figsize=(10, 6))
for p_val in p_values:
cdf_values = [float(geometric_cdf(f"{k}, {p_val}")) for k in k_values]
plt.plot(k_values, cdf_values, marker='o', label=f'p = {p_val}')
plt.xlabel('试验次数 k')
plt.ylabel('累积概率 F(k)')
plt.title('几何分布累积分布函数')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
except ImportError:
print("(需要matplotlib来显示图形)")
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import scipy.stats as st
import sympy as sp
from sympy.stats import Geometric, cdf
# 几何分布的累积分布函数
def geometric_cdf(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
# 计算卡方分布的概率密度函数 (PDF)
# 转换为 SymPy 数组（仅在符号分支使用）
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
k = int(expr[0]) # 实验的次数，标量或数组
p = float(expr[1]) # 成功的概率
# 计算概率值，处理k为负值的情况（设置概率为0）
result = np.where(k >= 0, st.geom.cdf(k, p), 0)
elif any(e.free_symbols for e in expr):
k, p = expr
X = Geometric('geom', p)
cdf_expr = cdf(X)(k)
result = cdf_expr.evalf()
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例使用:
# 标量输入
print(geometric_cdf("3, 0.5"))
# 0.875
# 符号变量输入
print(geometric_cdf("x, 0.5"))
# Piecewise((1.0 - 2.0*0.5**(floor(x) + 1.0), x >= 1), (0, True))
几何分布的概率质量函数
y=geopdf（x，p）使用p中的相应概率返回x中每个值处的几何分布的概率密度函数（pdf）。x和p可以是向量、矩阵或多维数组。标量输入被扩展为与其他输入具有相同维度的常数数组。p中的参数必须位于[0,1]区间内。
x -- 整数标量, 向量，矩阵，多维数组，是第一次成功之前的失败次数
p -- 标量，向量，矩阵，多维数组，是成功的概率
示例1. 质量控制 - 第一次成功检测分析
分析在生产线中第一次检测到次品的概率分布
假设次品率为2%
defect_rate = 0.02
print(f"生产线次品率: {defect_rate * 100}%")
#生产线次品率: 2.0%
在前N个产品中第一次发现次品的概率分布:
计算前20个产品中第一次发现次品的概率
k_values = range(1, 21)
probabilities = [geopdf(k, defect_rate) for k in k_values]
for k, prob in zip(k_values, probabilities):
print(f"在第{k:2d}个产品第一次发现次品的概率: {prob:.4f} ({prob * 100:.2f}%)")
#在第 1个产品第一次发现次品的概率: 0.0200 (2.00%)
#在第 2个产品第一次发现次品的概率: 0.0196 (1.96%)
#在第 3个产品第一次发现次品的概率: 0.0192 (1.92%)
#在第 4个产品第一次发现次品的概率: 0.0188 (1.88%)
#在第 5个产品第一次发现次品的概率: 0.0184 (1.84%)
#在第 6个产品第一次发现次品的概率: 0.0181 (1.81%)
#在第 7个产品第一次发现次品的概率: 0.0177 (1.77%)
#在第 8个产品第一次发现次品的概率: 0.0174 (1.74%)
#在第 9个产品第一次发现次品的概率: 0.0170 (1.70%)
#在第10个产品第一次发现次品的概率: 0.0167 (1.67%)
#在第11个产品第一次发现次品的概率: 0.0163 (1.63%)
#在第12个产品第一次发现次品的概率: 0.0160 (1.60%)
#在第13个产品第一次发现次品的概率: 0.0157 (1.57%)
#在第14个产品第一次发现次品的概率: 0.0154 (1.54%)
#在第15个产品第一次发现次品的概率: 0.0151 (1.51%)
#在第16个产品第一次发现次品的概率: 0.0148 (1.48%)
#在第17个产品第一次发现次品的概率: 0.0145 (1.45%)
#在第18个产品第一次发现次品的概率: 0.0142 (1.42%)
#在第19个产品第一次发现次品的概率: 0.0139 (1.39%)
#在第20个产品第一次发现次品的概率: 0.0136 (1.36%)
找到最可能第一次发现次品的位置
max_prob_index = np.argmax(probabilities)
print(f"\n最可能第一次发现次品的位置: 第{max_prob_index + 1}个产品")
print(f"对应的概率: {probabilities[max_prob_index]:.4f}")
#最可能第一次发现次品的位置: 第1个产品
#对应的概率: 0.0200
可视化概率分布
try:
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 6))
plt.bar(k_values, probabilities, alpha=0.7, color='skyblue')
plt.xlabel('检测的产品序号')
plt.ylabel('概率')
plt.title('第一次发现次品的概率分布')
plt.grid(True, alpha=0.3)
plt.show()
except ImportError:
print("(需要matplotlib来显示图形)")
示例2. 网络传输 - 数据包成功传输分析
分析数据包第k次尝试才成功传输的概率分布
假设单次传输成功率为80%
success_rate = 0.80
print(f"单次数据包传输成功率: {success_rate * 100}%")
#单次数据包传输成功率: 80.0%
第k次尝试才成功传输的概率分布:
分析前10次尝试的概率分布
attempt_values = range(1, 11)
pdf_values = [geopdf(k, success_rate) for k in attempt_values]
for attempt, prob in zip(attempt_values, pdf_values):
print(f"第{attempt}次尝试才成功的概率: {prob:.4f} ({prob * 100:.2f}%)")
#第1次尝试才成功的概率: 0.8000 (80.00%)
#第2次尝试才成功的概率: 0.1600 (16.00%)
#第3次尝试才成功的概率: 0.0320 (3.20%)
#第4次尝试才成功的概率: 0.0064 (0.64%)
#第5次尝试才成功的概率: 0.0013 (0.13%)
#第6次尝试才成功的概率: 0.0003 (0.03%)
#第7次尝试才成功的概率: 0.0001 (0.01%)
#第8次尝试才成功的概率: 0.0000 (0.00%)
#第9次尝试才成功的概率: 0.0000 (0.00%)
#第10次尝试才成功的概率: 0.0000 (0.00%)
计算期望传输次数
expected_attempts = 1 / success_rate
print(f"\n期望传输次数: {expected_attempts:.2f} 次")
#期望传输次数: 1.25 次
不同网络质量下的首次成功传输概率分布
network_conditions = {
"优秀网络 (成功率95%)": 0.95,
"良好网络 (成功率80%)": 0.80,
"一般网络 (成功率60%)": 0.60,
"较差网络 (成功率30%)": 0.30
}
for desc, rate in network_conditions.items():
prob_1st = geopdf(1, rate)
prob_3rd = geopdf(3, rate)
print(f"{desc}: 第1次成功概率={prob_1st:.4f}, 第3次成功概率={prob_3rd:.4f}")
#优秀网络 (成功率95%): 第1次成功概率=0.9500, 第3次成功概率=0.0024
#良好网络 (成功率80%): 第1次成功概率=0.8000, 第3次成功概率=0.0320
#一般网络 (成功率60%): 第1次成功概率=0.6000, 第3次成功概率=0.0960
#较差网络 (成功率30%): 第1次成功概率=0.3000, 第3次成功概率=0.1470
示例3. 销售漏斗 - 客户转化时机分析
分析客户在第几次接触时才转化的概率分布
假设每次营销接触的转化率为3%
conversion_rate = 0.03
print(f"单次营销接触转化率: {conversion_rate * 100}%")
#单次营销接触转化率: 3.0%
客户在第k次接触时才转化的概率分布:
分析前15次接触的概率分布
contact_values = range(1, 16)
pdf_values = [geopdf(k, conversion_rate) for k in contact_values]
print("接触次数 | 概率 | 累计概率")
print("-" * 35)
cumulative_prob = 0
for contact, prob in zip(contact_values, pdf_values):
cumulative_prob += prob
print(f"第{contact:2d}次 | {prob:.4f} | {cumulative_prob:.4f}")
#接触次数 | 概率 | 累计概率
#-----------------------------------
#第 1次 | 0.0300 | 0.0300
#第 2次 | 0.0291 | 0.0591
#第 3次 | 0.0282 | 0.0873
#第 4次 | 0.0274 | 0.1147
#第 5次 | 0.0266 | 0.1413
#第 6次 | 0.0258 | 0.1670
#第 7次 | 0.0250 | 0.1920
#第 8次 | 0.0242 | 0.2163
#第 9次 | 0.0235 | 0.2398
#第10次 | 0.0228 | 0.2626
#第11次 | 0.0221 | 0.2847
#第12次 | 0.0215 | 0.3062
#第13次 | 0.0208 | 0.3270
#第14次 | 0.0202 | 0.3472
#第15次 | 0.0196 | 0.3667
找到转化概率最高的接触点
max_prob_contact = np.argmax(pdf_values) + 1
print(f"\n转化概率最高的接触点: 第{max_prob_contact}次接触")
#转化概率最高的接触点: 第1次接触
基于概率分布的营销策略建议:
high_prob_contacts = [i + 1 for i, p in enumerate(pdf_values) if p > 0.02]
print(f"应重点优化的接触点: {high_prob_contacts}")
#应重点优化的接触点: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
示例4. 游戏设计 - 稀有物品掉落分析
分析玩家在第几次尝试才获得稀有物品的概率分布
假设稀有物品掉落率为1%
drop_rate = 0.01
print(f"稀有物品掉落率: {drop_rate * 100}%")
#稀有物品掉落率: 1.0%
玩家在第k次尝试才获得稀有物品的概率分布:
分析概率分布
attempt_values = range(1, 201, 10) # 每10次尝试采样一次
pdf_values = [geopdf(k, drop_rate) for k in attempt_values]
for attempt, prob in zip(attempt_values, pdf_values):
if prob > 0.0001: # 只显示概率大于0.01%的情况
print(f"第{attempt:3d}次尝试才获得的概率: {prob:.6f}")
#第 1次尝试才获得的概率: 0.010000
#第 11次尝试才获得的概率: 0.009044
#第 21次尝试才获得的概率: 0.008179
#第 31次尝试才获得的概率: 0.007397
#第 41次尝试才获得的概率: 0.006690
#第 51次尝试才获得的概率: 0.006050
#第 61次尝试才获得的概率: 0.005472
#第 71次尝试才获得的概率: 0.004948
#第 81次尝试才获得的概率: 0.004475
#第 91次尝试才获得的概率: 0.004047
#第101次尝试才获得的概率: 0.003660
#第111次尝试才获得的概率: 0.003310
#第121次尝试才获得的概率: 0.002994
#第131次尝试才获得的概率: 0.002708
#第141次尝试才获得的概率: 0.002449
#第151次尝试才获得的概率: 0.002215
#第161次尝试才获得的概率: 0.002003
#第171次尝试才获得的概率: 0.001811
#第181次尝试才获得的概率: 0.001638
#第191次尝试才获得的概率: 0.001481
分析"幸运"和"不幸"玩家的概率
幸运玩家（前10次获得）
lucky_prob = sum([geopdf(k, drop_rate) for k in range(1, 11)])
print(f"幸运玩家（前10次获得）的概率: {lucky_prob:.4f}")
#幸运玩家（前10次获得）的概率: 0.0956
普通玩家（11-100次获得）
normal_prob = sum([geopdf(k, drop_rate) for k in range(11, 101)])
print(f"普通玩家（11-100次获得）的概率: {normal_prob:.4f}")
#普通玩家（11-100次获得）的概率: 0.5383
不幸玩家（100次后获得）
unlucky_prob = 1 - lucky_prob - normal_prob
print(f"不幸玩家（100次后获得）的概率: {unlucky_prob:.4f}")
#不幸玩家（100次后获得）的概率: 0.3660
保底机制设计
pity_threshold = 100
natural_prob = sum([geopdf(k, drop_rate) for k in range(1, pity_threshold + 1)])
print(f"\n前{pity_threshold}次自然获得的概率: {natural_prob:.4f}")
print(f"需要保底机制的玩家比例: {(1 - natural_prob) * 100:.2f}%")
#前100次自然获得的概率: 0.6340
#需要保底机制的玩家比例: 36.60%
示例5. 医学诊断 - 疾病检测分析
假设单次检测的确诊率为85%
diagnosis_rate = 0.85
print(f"单次检测确诊率: {diagnosis_rate * 100}%")
#单次检测确诊率: 85.0%
在第k次检测才确诊的概率分布:
分析前5次检测的概率分布
test_values = range(1, 6)
pdf_values = [geopdf(k, diagnosis_rate) for k in test_values]
for test, prob in zip(test_values, pdf_values):
print(f"第{test}次检测才确诊的概率: {prob:.4f} ({prob * 100:.2f}%)")
#第1次检测才确诊的概率: 0.8500 (85.00%)
#第2次检测才确诊的概率: 0.1275 (12.75%)
#第3次检测才确诊的概率: 0.0191 (1.91%)
#第4次检测才确诊的概率: 0.0029 (0.29%)
#第5次检测才确诊的概率: 0.0004 (0.04%)
不同检测方法的效率比较:
test_methods = {
"高精度检测 (95%)": 0.95,
"标准检测 (85%)": 0.85,
"快速检测 (70%)": 0.70,
"筛查检测 (50%)": 0.50
}
for desc, rate in test_methods.items():
prob_1st = geopdf(1, rate)
prob_2nd = geopdf(2, rate)
expected_tests = 1 / rate
print(f"{desc}: 第1次确诊概率={prob_1st:.4f}, 期望检测次数={expected_tests:.2f}")
#高精度检测 (95%): 第1次确诊概率=0.9500, 期望检测次数=1.05
#标准检测 (85%): 第1次确诊概率=0.8500, 期望检测次数=1.18
#快速检测 (70%): 第1次确诊概率=0.7000, 期望检测次数=1.43
#筛查检测 (50%): 第1次确诊概率=0.5000, 期望检测次数=2.00
示例6. 客户服务 - 问题解决分析
假设单次联系解决问题的概率为40%
resolution_rate = 0.40
print(f"单次联系问题解决率: {resolution_rate * 100}%")
#单次联系问题解决率: 40.0%
客户问题在第k次联系才解决的概率分布:
分析前8次联系的概率分布
contact_values = range(1, 9)
pdf_values = [geopdf(k, resolution_rate) for k in contact_values]
cumulative_prob = 0
for contact, prob in zip(contact_values, pdf_values):
cumulative_prob += prob
print(f"第{contact}次联系解决: 概率={prob:.4f}, 累计概率={cumulative_prob:.4f}")
#第1次联系解决: 概率=0.4000, 累计概率=0.4000
#第2次联系解决: 概率=0.2400, 累计概率=0.6400
#第3次联系解决: 概率=0.1440, 累计概率=0.7840
#第4次联系解决: 概率=0.0864, 累计概率=0.8704
#第5次联系解决: 概率=0.0518, 累计概率=0.9222
#第6次联系解决: 概率=0.0311, 累计概率=0.9533
#第7次联系解决: 概率=0.0187, 累计概率=0.9720
#第8次联系解决: 概率=0.0112, 累计概率=0.9832
服务质量评估
一次联系解决的客户比例
first_contact_resolution = geopdf(1, resolution_rate)
print(f"一次联系解决率: {first_contact_resolution * 100:.1f}%")
#一次联系解决率: 40.0%
三次联系内解决的客户比例
three_contact_resolution = sum([geopdf(k, resolution_rate) for k in range(1, 4)])
print(f"三次联系内解决率: {three_contact_resolution * 100:.1f}%")
#三次联系内解决率: 78.4%
需要升级处理的客户比例
escalation_threshold = 5
escalation_prob = 1 - sum([geopdf(k, resolution_rate) for k in range(1, escalation_threshold + 1)])
print(f"需要升级处理的客户比例: {escalation_prob * 100:.1f}%")
#需要升级处理的客户比例: 7.8%
示例7. 科学研究 - 实验成功分析
假设单次实验成功的概率为20%
experiment_success_rate = 0.20
print(f"单次实验成功率: {experiment_success_rate * 100}%")
#单次实验成功率: 20.0%
在第k次实验才成功的概率分布:
分析概率分布
experiment_values = range(1, 16)
pdf_values = [geopdf(k, experiment_success_rate) for k in experiment_values]
for exp, prob in zip(experiment_values, pdf_values):
print(f"第{exp:2d}次实验才成功的概率: {prob:.4f}")
#第 1次实验才成功的概率: 0.2000
#第 2次实验才成功的概率: 0.1600
#第 3次实验才成功的概率: 0.1280
#第 4次实验才成功的概率: 0.1024
#第 5次实验才成功的概率: 0.0819
#第 6次实验才成功的概率: 0.0655
#第 7次实验才成功的概率: 0.0524
#第 8次实验才成功的概率: 0.0419
#第 9次实验才成功的概率: 0.0336
#第10次实验才成功的概率: 0.0268
#第11次实验才成功的概率: 0.0215
#第12次实验才成功的概率: 0.0172
#第13次实验才成功的概率: 0.0137
#第14次实验才成功的概率: 0.0110
#第15次实验才成功的概率: 0.0088
科研资源规划
期望实验次数
expected_experiments = 1 / experiment_success_rate
print(f"期望实验次数: {expected_experiments:.1f} 次")
#期望实验次数: 5.0 次
不同置信度下需要的最大实验次数
confidence_levels = [0.5, 0.8, 0.9, 0.95]
for confidence in confidence_levels:
max_experiments = np.ceil(np.log(1 - confidence) / np.log(1 - experiment_success_rate))
actual_prob = sum(
[geopdf(k, experiment_success_rate) for k in range(1, int(max_experiments) + 1)])
print(f"{confidence * 100}%置信度: 最多需要{int(max_experiments)}次实验 (实际概率: {actual_prob:.4f})")
#50.0%置信度: 最多需要4次实验 (实际概率: 0.5904)
#80.0%置信度: 最多需要8次实验 (实际概率: 0.8322)
#90.0%置信度: 最多需要11次实验 (实际概率: 0.9141)
#95.0%置信度: 最多需要14次实验 (实际概率: 0.9560)
不同难度实验的成功时机比较:
difficulty_levels = {
"简单实验 (50%成功率)": 0.50,
"中等实验 (20%成功率)": 0.20,
"困难实验 (5%成功率)": 0.05,
"极难实验 (1%成功率)": 0.01
}
for desc, rate in difficulty_levels.items():
prob_5th = geopdf(5, rate)
expected = 1 / rate
print(f"{desc}: 第5次成功概率={prob_5th:.4f}, 期望次数={expected:.1f}")
#简单实验 (50%成功率): 第5次成功概率=0.0312, 期望次数=2.0
#中等实验 (20%成功率): 第5次成功概率=0.0819, 期望次数=5.0
#困难实验 (5%成功率): 第5次成功概率=0.0407, 期望次数=20.0
#极难实验 (1%成功率): 第5次成功概率=0.0096, 期望次数=100.0
示例8. 符号计算与理论分析
计算几何分布的PMF符号表达式
symbolic_pmf = geopdf(k, p)
print(f"几何分布PMF的符号表达式: {symbolic_pmf}")
#几何分布PMF的符号表达式: p*(1.0 - p)**(k - 1.0)
验证概率和为1
理论上：∑_{k=1}^∞ (1-p)^{k-1} * p = 1
theoretical_sum = sp.Sum(p * (1 - p) ** (k - 1), (k, 1, sp.oo))
sum_value = theoretical_sum.doit()
print(f"理论概率和: {sum_value}")
#理论概率和: p*Piecewise((1/p, Abs(p - 1) < 1), (Sum((1 - p)**(k - 1), (k, 1, oo)), True))
计算期望和方差
期望值 E[X] = 1/p
expected_value = 1 / p
print(f"期望值 E[X] = {expected_value}")
#期望值 E[X] = 1/p
方差 Var[X] = (1-p)/p²
variance = (1 - p) / (p ** 2)
print(f"方差 Var[X] = {variance}")
#方差 Var[X] = (1 - p)/p**2
不同参数下的概率分布可视化
try:
import matplotlib.pyplot as plt
# 数值计算示例
k_values = np.arange(1, 21)
p_values = [0.1, 0.3, 0.5, 0.7, 0.9]
plt.figure(figsize=(12, 8))
for p_val in p_values:
pmf_values = [float(geopdf(k, p_val)) for k in k_values]
plt.plot(k_values, pmf_values, marker='o', label=f'p = {p_val}')
plt.xlabel('试验次数 k')
plt.ylabel('概率 P(X=k)')
plt.title('几何分布概率质量函数')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
# 累积分布可视化
plt.figure(figsize=(12, 8))
for p_val in p_values:
pmf_values = [float(geopdf(k, p_val)) for k in k_values]
cdf_values = np.cumsum(pmf_values)
plt.plot(k_values, cdf_values, marker='s', label=f'p = {p_val}')
plt.xlabel('试验次数 k')
plt.ylabel('累积概率 P(X≤k)')
plt.title('几何分布累积分布函数')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
except ImportError:
print("(需要matplotlib来显示图形)")
无记忆性验证
几何分布具有无记忆性：P(X>m+n | X>m) = P(X>n)
conditional_prob = (1 - p) ** (m + n) / (1 - p) ** m
unconditional_prob = (1 - p) ** n
print(f"条件概率 P(X>{m + n}|X>{m}) = {conditional_prob}")
print(f"无条件概率 P(X>{n}) = {unconditional_prob}")
print(f"无记忆性成立: {sp.simplify(conditional_prob - unconditional_prob) == 0}")
#条件概率 P(X>m + n|X>m) = (1 - p)**(m + n)/(1 - p)**m
#无条件概率 P(X>n) = (1 - p)**n
#无记忆性成立: True
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import scipy.stats as st
import sympy as sp
from sympy.stats import Geometric, density
# 几何分布的概率质量函数
def geometric_pdf(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
# 计算卡方分布的概率密度函数 (PDF)
# 转换为 SymPy 数组（仅在符号分支使用）
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
k = int(expr[0]) # 实验的次数，标量或数组
p = float(expr[1]) # 成功的概率
# 计算概率值，处理k为负值的情况（设置概率为0）
result = np.where(k >= 0, st.geom.pmf(k, p), 0)
elif any(e.free_symbols for e in expr):
k, p = expr
X = Geometric('geom', p)
density_expr = density(X)(k)
result = density_expr.evalf()
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例使用:
# 标量输入
print(geometric_pdf("3, 0.5"))
# 0.125
# 符号变量输入
print(geometric_pdf("x, 0.5"))
# 0.5*0.5**(x - 1.0)
平面几何图
GeoPlot(Point2d,Segment,Circle,Ellipse,Polygon，Fill,Legend)绘制点，线，圆，椭圆及多边形。
Point2d[a,b]绘制(a,b)坐标的点。
Segment(a1,b1,a2,b2)绘制(a1,b1)到(a2,b2)的线条。
Circle(a,b,r)绘制以(a,b)为圆心，r为半径的圆。
Ellipse(a,b,r,vr)绘制以(a,b)为圆心，r为半径, vr为垂直半径的椭圆。
Polygon(a1,b1,a2,b2,a3,b3....)绘制多点之间的多边形。
Fill: 是否在几何体内填充颜色。
Legend: 是否有注释
立体几何图
GeoPlot3D(Point3d,Line3d,Plane,Legend)绘制点，线，平面。
Point3d[x,y,z]绘制(x,y,z)空间坐标的点。
Line3d(x1,y1,z1,x2,y2,z2)绘制空间(x1,y1,z1)到(x2,y2,z2)的线条。
Plane(x1,y1,z1,x2,y2,z2,(x_range),(y_range),(z_range))绘制空间平面。
Legend: 是否有注释
几何分布的随机采样函数
r=geornd（p）从具有概率参数p的几何分布中生成随机数。p可以是向量、矩阵或多维数组。r的大小等于p的大小。p中的参数必须位于[0,1]区间内。
r=geornd（p，m，n，…）或r=georn d（p、[m、n，…]）生成多维m-by-n-by-。..包含来自具有概率参数p的几何分布的随机数的数组。p可以是标量或与r大小相同的数组。
几何分布有助于在一系列独立试验中模拟一次成功之前的失败次数，其中每次试验都会导致成功或失败，任何单个试验的成功概率都是常数p。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy.stats import geom
def geornd_scipy(p, size=None):
"""
几何分布随机采样函数（使用 SciPy）
返回首次成功前的失败次数（k=0,1,2,...）
参数:
p: 成功概率（标量或矩阵）
size: 输出形状（整数或元组）
返回:
随机数或随机数矩阵
"""
# 处理标量概率
if np.isscalar(p):
if not 0 < p <= 1:
raise ValueError("概率 p 必须在 (0,1] 区间")
# 使用 SciPy 的几何分布生成器（p 参数对应失败概率）
return geom.rvs(p, size=size) - 1
# 处理矩阵概率
p_arr = np.asarray(p)
if p_arr.ndim == 0:
return geom.rvs(p_arr, size=size) - 1
# 检查所有概率值是否有效
if np.any((p_arr <= 0) | (p_arr > 1)):
raise ValueError("所有概率值必须在 (0,1] 区间")
# 为每个概率值生成随机数
return np.vectorize(lambda p_val: geom.rvs(p_val) - 1)(p_arr)
# 示例用法
if __name__ == "__main__":
np.random.seed(42) # 设置随机种子
# 示例1：标量概率（单个随机数）
print("单个随机数:", geornd_scipy(0.3))
# 1
# 示例2：指定维度（2x3矩阵）
print("\n2x3矩阵:")
print(geornd_scipy(0.5, size=(2, 3)))
# [[4 1 1]
# [0 0 0]]
# 示例3：矩阵概率
prob_matrix = [[0.2, 0.5], [0.3, 0.7]]
print("\n矩阵概率:")
print(geornd_scipy(prob_matrix))
# [[4 1]
# [0 2]]
# 示例4：不同形状
print("\n一维数组:")
print(geornd_scipy(0.6, size=5))
# [1 0 0 0 0]
冈珀茨拟合
gompertz: Y = d+(a-d)*exp(-exp(-b*(x-c)))
返回冈珀茨拟合方程式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.optimize import curve_fit
def gompertz_fit_nonlinear(input_str):
try:
expr = sp.sympify(input_str)
error = False
expression = None
maxfev = 100000
x_sym = sp.symbols('x') # 定义符号变量x
# 定义 Gompertz 函数
def gompertz(x, a, b, c, d):
return d + (a - d) * np.exp(-np.exp(-b * (x - c)))
# 检查输入是否为包含两个元素的元组
if isinstance(expr, tuple) and len(expr) == 2:
x_part, y_part = expr[0], expr[1]
# 转换为矩阵并验证
if not isinstance(x_part, list) or not isinstance(y_part, list):
error = True
else:
x_data = np.array(x_part, dtype=float).ravel()
y_data = np.array(y_part, dtype=float).ravel()
if len(x_data) != len(y_data):
error = True
else:
# 初始参数猜测
initial_guess = [max(y_data), 1, np.mean(x_data), min(y_data)]
# 非线性拟合
params, _ = curve_fit(
gompertz,
x_data,
y_data,
p0=initial_guess,
maxfev=maxfev
)
# 提取参数并构造Sympy表达式
a, b, c, d = params
expression = f"{d:.4f} + ({a:.4f} - {d:.4f}) * exp(-exp(-{b:.4f} * (x - {c:.4f})))"
else:
error = True # 输入格式不正确
if error:
return f"输入错误: {input_str}"
else:
return sp.simplify(expression)
except Exception as e:
return f"错误: {e}"
# 示例输入：x和y数据点
input_data = "([1,2,3,4,5], [2.0, 4.5, 6.5, 8.0, 8.5])"
result = gompertz_fit_nonlinear(input_data)
print(result)
# 0.4148 + 8.7322*exp(-3.77555321883638*exp(-0.7963*x))
符号标量场的梯度向量(数值梯度)
FX = gradient(F) 返回向量F的一维数值梯度.输出FX对应于∂F/∂x,即x(水平)方向上的差分.点之间的间距假定为1
FX = gradient(F,[vars],[points]) 返回符号表达式列表的符号标量场的梯度向量.
[FX,FY] = gradient(F) 返回矩阵F的二维数值梯度的x和y分量.附加输出FY对应于∂F/∂y,即y(垂直)方向上的差分.每个方向上的点之间的间距假定为1
F — 输入数组,向量,矩阵,符号表达式
vars -- 符号变量列表.
points -- 数值列表.
FX, FY — 数值梯度,数组
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def get_result_value(symbols_set, points=None, result=None):
"""
将符号表达式代入具体数值点求值。
参数:
symbols_set: 符号集合
points: 数值点列表，与符号一一对应
result: 待代入的表达式
返回:
代入后的数值结果或符号表达式
"""
if points is not None and result is not None:
if len(points) != len(symbols_set):
raise ValueError(f"提供的点数（{len(points)}）与符号数量（{len(symbols_set)}）不匹配。")
subs_dict = {symbol: value for symbol, value in zip(symbols_set, points)}
try:
# 尝试数值求值
result = result.subs(subs_dict).evalf()
except Exception as e:
# 非数值情况直接代入
result = result.subs(subs_dict)
return result
def gradient_by_function(input_str):
"""
对标 MATLAB 的 gradient 函数，计算数值或符号梯度。
MATLAB梯度规则:
- 对于向量：中心差分（边界使用单边差分）
- 对于矩阵：对每个维度分别计算梯度
参数:
input_str (str): 输入表达式，可以是以下形式：
1. 纯矩阵 (如 "Matrix([[1,2],[3,4]])")
2. 符号表达式 (如 "x**2 + y**3")
3. 元组结构 (如 "(f, [x,y], [1,2])"）
其中最后一个元素为可选的代入点
返回:
SymPy 矩阵/表达式 或错误信息字符串
"""
try:
# 将输入的字符串转换为 SymPy 表达式
expr = sp.sympify(input_str)
# 初始化错误标志为 False
error = False
# 初始化结果变量
result = None
def remove_last_from_tuple(tup):
"""
移除元组的最后一个元素。
参数:
tup (tuple): 待处理的元组
返回:
tuple: 移除最后一个元素后的元组
"""
return tup[:-1] if len(tup) > 1 else tup
# 检查输入是否为元组，且最后一个元素为数值列表
if isinstance(expr, tuple) and isinstance(expr[-1], list) and all(item.is_number for item in expr[-1]):
# 提取代入点
points = expr[-1]
# 移除元组的最后一个元素
expr = remove_last_from_tuple(expr)
else:
# 如果不满足条件，代入点设为 None
points = None
# 检查输入是否为包含两个元素的元组
if isinstance(expr, tuple) and len(expr) == 2:
# 提取函数和变量
func, variables = expr
# 计算函数关于变量的梯度
grad = sp.derive_by_array(func, variables)
# 如果梯度结果长度大于 1，转换为 SymPy 矩阵
result = sp.Matrix(grad) if len(grad) > 1 else grad
# 将变量代入梯度结果
result = get_result_value(variables, points, result)
# 检查输入是否为 SymPy 矩阵
elif isinstance(expr, list):
# 检查矩阵是否为数值矩阵
sp_matrix = sp.Matrix(expr)
if sp_matrix.is_symbolic() is False:
# 检查矩阵是否为向量（行数或列数为 1）
is_vector = sp_matrix.rows == 1 or sp_matrix.cols == 1
if is_vector:
# 将向量转换为一维 NumPy 数组
A = np.ravel(np.array(sp_matrix, dtype=float))
# 计算向量的梯度
np_grad = np.gradient(A)
else:
# 将矩阵转换为 NumPy 数组
np_array = np.array(sp_matrix.tolist(), dtype=float)
# 计算矩阵的梯度
np_grad = np.gradient(np_array)
# 将 NumPy 梯度结果转换为 SymPy 矩阵
if isinstance(np_grad, list):
result = [sp.Matrix(g) for g in np_grad]
else:
result = sp.Matrix(np_grad)
else:
# 提取矩阵中的自由符号
symbols = sp_matrix.free_symbols
# 计算矩阵关于符号的梯度
result = sp.derive_by_array(sp_matrix, list(symbols))
# 将梯度结果转换为 SymPy 矩阵
result = sp.Matrix(result)
# 将变量代入梯度结果
result = get_result_value(symbols, points, result)
# 检查输入表达式是否包含自由符号
elif expr.free_symbols:
# 提取表达式中的自由符号
symbols_set = list(expr.free_symbols)
# 计算表达式关于符号的梯度
result = sp.derive_by_array(expr, symbols_set)
# 将梯度结果转换为 SymPy 矩阵
result = sp.Matrix(result)
# 将变量代入梯度结果
result = get_result_value(symbols_set, points, result)
else:
# 如果输入表达式不包含自由符号，设置错误标志为 True
error = True
# 如果没有错误，返回结果；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"输入错误:{e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 符号表达式梯度
print("\n示例1: 符号梯度")
x, y = sp.symbols('x y')
input_expr = "x**2 + y**3"
grad = gradient_by_function(input_expr)
print(f"\n表达式: {input_expr}")
print("符号梯度:")
print(grad)
# Matrix([[2*x],
# [3*y**2]])
# 示例2: 带代入点的符号梯度
print("\n示例2: 带代入点的梯度")
input_tuple = "(x**2 + y**3, [x, y], [1, 2])"
grad = gradient_by_function(input_tuple)
print(f"\n输入：{input_tuple}")
print("代入点后的数值梯度:")
print(grad)
# Matrix([[2.00000000000000],
# [12.0000000000000]])
具有敏感特征值的托普利茨矩阵
A = grcar(n,k) 返回 n×n 托普利茨矩阵, 其中下对角线上的元素为 -1, 主对角线上的元素为 1, 主对角线上方的 k 个对角线上的元素为 1. k 必须为整数, 默认值为 k = 3. A 具有对扰动敏感的特征值.
n,k — 输入,标量
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def create_grcar_matrix(n, k=3):
"""
生成一个n×n的Grcar矩阵。
Grcar矩阵是一种具有特殊结构的测试矩阵，主对角线为1，下方第一条对角线为-1，
上方k条对角线为1，其余位置为0。
参数：
n (int): 矩阵的维度。
k (int, 可选): 主对角线上方的对角线数量，默认为3。
返回：
numpy.ndarray: 生成的Grcar矩阵。
"""
# 创建n×n的零矩阵
A = np.zeros((n, n), dtype=int)
# 设置主对角线为1
np.fill_diagonal(A, 1)
# 设置下方第一条对角线为-1
# A[1:]表示从第二行开始的所有行，填充其主对角线（即原矩阵的下第一条对角线）
np.fill_diagonal(A[1:], -1)
# 设置主对角线上方的k条对角线为1
for diag in range(1, k + 1):
# 对每行，从第diag列开始，填充其主对角线（即原矩阵右上方的第diag条对角线）
np.fill_diagonal(A[:, diag:], 1)
return A
def grcar_matrix(input_str):
"""
根据输入字符串生成Grcar矩阵。
输入可以是单个数字n（如"5"）或元组(n, k)（如"(5, 3)"）。
参数：
input_str (str): 描述矩阵参数的字符串。
返回：
sympy.Matrix or str: 生成的矩阵或错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple):
# 处理元组输入，例如"(5, 3)"
if len(expr) != 2:
error = True
else:
n = int(expr[0])
k = int(expr[1])
result = create_grcar_matrix(n, k)
elif expr.is_number:
# 处理单个数字输入，例如"5"
n = int(expr)
result = create_grcar_matrix(n)
else:
error = True
if error:
return f"输入格式错误: {input_str}"
else:
return sp.Matrix(result)
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：生成5x5，k=3的Grcar矩阵
print("示例1：n=5, k=3")
print(grcar_matrix("5"))
# Matrix([[1, 1, 1, 1, 0],
# [-1, 1, 1, 1, 1],
# [0, -1, 1, 1, 1],
# [0, 0, -1, 1, 1],
# [0, 0, 0, -1, 1]])
# 示例2：生成5x5，k=2的Grcar矩阵
print("\n示例2：n=5, k=2")
print(grcar_matrix("(5, 2)"))
# Matrix([[1, 1, 1, 0, 0],
# [-1, 1, 1, 1, 0],
# [0, -1, 1, 1, 1],
# [0, 0, -1, 1, 1],
# [0, 0, 0, -1, 1]])
广义奇异值分解
[U,V,S] = gsvd(A,B) 执行矩阵A和B的广义奇异值分解,并返回酉矩阵U和V以及非负对角矩阵S
A, B — 输入矩阵(以单独参量指定),矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.linalg import qr, svd
def generalized_singular_value_decomposition(input_str):
"""
执行广义奇异值分解（GSVD），对标MATLAB的gsvd函数。
参数:
input_str: 输入字符串，表示两个矩阵的元组，例如 "([[1,2],[3,4]], [[5,6],[7,8]])"
返回:
返回包含U, V, X, C, S的元组（SymPy矩阵），或错误信息。
"""
try:
# 解析输入字符串为SymPy表达式
expr = sp.sympify(input_str)
if not isinstance(expr, tuple) or len(expr) != 2:
return "输入错误：请输入包含两个矩阵的元组。"
if any(isinstance(e, list) is False for e in expr):
return "输入错误：无法解析为两个有效矩阵。"
# 转换为NumPy数组
A_np = np.array(expr[0], dtype=float)
B_np = np.array(expr[1], dtype=float)
# 检查矩阵列数是否相同
if A_np.shape[1] != B_np.shape[1]:
return f"错误：矩阵A和B的列数必须相同。A的列数：{A_np.shape[1]}, B的列数：{B_np.shape[1]}"
# 计算矩阵B的QR分解（经济模式）
# B^T = Q_B * R_B → B = R_B^T @ Q_B^T
Q_B, R_B = qr(B_np.T, mode='economic')
R_B = R_B.T # 调整后B = Q_B @ R_B
# 计算A * Q_B^T
AQ_B = A_np @ Q_B.T
# 对AQ_B进行奇异值分解（SVD）
U, D, Zt = svd(AQ_B, full_matrices=False)
Z = Zt.T
# 计算X矩阵：X = Q_B^T @ Z
X = Q_B.T @ Z
# 计算C和S的对角元素
C_diag = D
S_diag = np.sqrt(1 - D ** 2)
# 构造对角矩阵C和S
C = np.diag(C_diag)
S = np.diag(S_diag)
# 计算V矩阵：B = V @ S @ X^T → V = B @ X @ inv(S)
# 由于B = Q_B @ R_B，代入得 Q_B @ R_B = V @ S @ X^T
# 结合X = Q_B^T @ Z → Q_B = X @ Z^T
V = (R_B @ Z) @ np.diag(1.0 / S_diag)
# 转换为SymPy矩阵
U_sp = sp.Matrix(U)
V_sp = sp.Matrix(V)
X_sp = sp.Matrix(X.T) # X在分解式中转置
C_sp = sp.Matrix(C)
S_sp = sp.Matrix(S)
return (U_sp, V_sp, X_sp, C_sp, S_sp)
except Exception as e:
return f"错误：{e}"
# 示例使用
input_str = "([[1,6,11],[2,7,12],[3,8,13],[4,9,14],[5,10,15]], [[8,1,6],[3,5,7],[4,9,2]])"
result = generalized_singular_value_decomposition(input_str)
if isinstance(result, tuple):
U, V, X, C, S = result
print("U =")
print(U)
# Matrix([[-0.354557057037681, 0.688686643768252, -0.100420070696462],
# [-0.398696369998832, 0.375554529395871, 0.401020580634150],
# [-0.442835682959984, 0.0624224150234908, -0.685649669614349],
# [-0.486974995921135, -0.250709699348890, 0.569917880112097],
# [-0.531114308882287, -0.563841813721270, -0.184868720435436]])
print("\nV =")
print(V)
# Matrix([[nan, nan, 8.64112621464805],
# [nan, nan, 4.34326123905226],
# [nan, nan, 0.200526501807073]])
print("\nX =")
print(X)
# Matrix([[-0.201664911192694, -0.516830501392304, -0.831996091591915],
# [-0.890317132783019, -0.257331626824051, 0.375653879134918],
# [0.408248290463863, -0.816496580927726, 0.408248290463863]])
print("\nC =")
print(C)
# Matrix([[35.1272233335747, 0, 0],
# [0, 2.46539669691652, 0],
# [0, 0, 1.17309315627448e-15]])
print("\nS =")
print(S)
# Matrix([[nan, 0, 0],
# [0, nan, 0],
# [0, 0, 1.00000000000000]])
else:
print(result)
古德曼函数
Y = Gudermannian(x)得到古德曼函数结果,将三角函数和双曲函数连系起来.
x -- 输入,数值,符号函数,向量,矩阵.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def gudermannian_function(input_str):
"""
计算输入表达式的 Gudermannian 函数值，支持标量、矩阵和符号表达式。
Gudermannian 函数定义为：gd(x) = 2 * arctan(exp(x)) - π/2
参数:
input_str (str): 输入的数学表达式字符串，可以是标量、矩阵或符号。
返回:
计算结果 (SymPy 表达式/矩阵) 或错误信息字符串。
"""
try:
# 将输入字符串解析为 SymPy 表达式
expr = sp.sympify(input_str)
error = False
result = None
# 处理元组类型输入（通常为无效输入）
if expr.free_symbols:
# 符号表达式处理：返回符号表达式
result = sp.simplify(2 * sp.atan(sp.exp(expr)) - sp.pi / 2)
elif expr.is_number:
z = complex(expr)
result = np.arcsin(np.tanh(z))
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"输入错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例 1: 标量数值输入
print("示例 1:")
print(gudermannian_function("1+2j"))
# (1.2304675059389543-0.6765532587260855j)
# 示例 2: 符号输入
print("\n示例 2:")
x = sp.symbols('x')
print(gudermannian_function("x"))
# 2*atan(exp(x)) - pi/2