负二项式累积分布函数
y=nbincdf(x,R,p)使用相应的成功次数,R和单次试验的成功概率计算x中每个值的负二项式cdf.
x,R和p可以是向量或矩阵,它们都具有相同的大小,也就是y的大小.
x,R或p的标量输入被扩展为与其他输入具有相同维度的常数数组.
import matplotlib.pyplot as plt
import numpy as np
示例1. 质量控制 - 缺陷产品检测
在质量控制中，负二项分布用于建模在找到第R个缺陷产品前需要检查的产品数量
场景：在找到第3个缺陷产品前，检查不超过10个产品的概率
result = nbincdf(10, 3, 0.05) # 缺陷率5%
print(f"在找到第3个缺陷品前检查≤10个产品的概率: {result:.4f}")
#在找到第3个缺陷品前检查≤10个产品的概率: 0.0245
不同缺陷率下的比较
defect_rates = [0.02, 0.05, 0.1, 0.15]
for rate in defect_rates:
prob = nbincdf(10, 3, rate)
print(f"缺陷率{rate:.0%}时概率: {prob:.4f}")
#缺陷率2%时概率: 0.0020
#缺陷率5%时概率: 0.0245
#缺陷率10%时概率: 0.1339
#缺陷率15%时概率: 0.3080
可视化
x_values = list(range(1, 31))
probs_5pct = [nbincdf(f"({x}, 3, 0.05)") for x in x_values]
probs_10pct = [nbincdf(f"({x}, 3, 0.1)") for x in x_values]
plt.figure(figsize=(10, 6))
plt.plot(x_values, probs_5pct, 'b-', label='缺陷率 5%', linewidth=2)
plt.plot(x_values, probs_10pct, 'r-', label='缺陷率 10%', linewidth=2)
plt.axhline(0.95, color='gray', linestyle='--', alpha=0.7, label='95%置信水平')
plt.xlabel('检查的产品数量')
plt.ylabel('累积概率')
plt.title('找到第3个缺陷品前的累积分布')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
示例2. 销售与市场营销 - 客户转化分析
在销售中，负二项分布用于建模获得第R个成功转化前需要的接触次数
场景：在获得第5个订单前，进行不超过50次销售接触的概率
conversion_rate = 0.1 # 10%转化率
result = nbincdf(50, 5, conversion_rate)
print(f"获得第5个订单前进行≤50次接触的概率: {result:.4f}")
#获得第5个订单前进行≤50次接触的概率: 0.6549
计算不同转化率下的效率
conversion_rates = [0.05, 0.08, 0.1, 0.12, 0.15]
target_contacts = 50
不同转化率下的效率比较:
for rate in conversion_rates:
prob = nbincdf(target_contacts, 5, rate)
print(f"转化率{rate:.1%}: {prob:.1%}概率在{target_contacts}次接触内获得5个订单")
#转化率5.0%: 14.0%概率在50次接触内获得5个订单
#转化率8.0%: 45.2%概率在50次接触内获得5个订单
#转化率10.0%: 65.5%概率在50次接触内获得5个订单
#转化率12.0%: 80.5%概率在50次接触内获得5个订单
#转化率15.0%: 93.0%概率在50次接触内获得5个订单
资源规划分析
required_successes = [3, 5, 8, 10]
confidence_level = 0.9
print(f"\n达到{confidence_level:.0%}置信水平所需接触次数:")
for r in required_successes:
contacts = 1
while True:
prob = nbincdf(contacts, r, conversion_rate)
if prob >= confidence_level:
print(f"获得{r}个订单: {contacts}次接触")
break
contacts += 1
#达到90%置信水平所需接触次数:
#获得3个订单: 49次接触
#获得5个订单: 73次接触
#获得8个订单: 108次接触
#获得10个订单: 130次接触
示例3. 医学研究 - 临床试验设计
在临床试验中，负二项分布用于建模观察到第R个有效病例前需要筛选的患者数量
场景：某种疾病的有效治疗响应率为30%，需要找到10个有效病例
response_rate = 0.3
required_successes = 10
计算在不同患者数量下的累积概率
patient_counts = [20, 30, 40, 50, 60]
print(f"找到{required_successes}个有效病例的累积概率:")
for n in patient_counts:
prob = nbincdf(n, required_successes, response_rate)
print(f"筛选{n}名患者: {prob:.1%}")
#找到10个有效病例的累积概率:
#筛选20名患者: 41.1%
#筛选30名患者: 80.4%
#筛选40名患者: 96.0%
#筛选50名患者: 99.4%
#筛选60名患者: 99.9%
样本量确定
confidence_levels = [0.8, 0.9, 0.95]
不同置信水平下的样本量要求:
for conf in confidence_levels:
n = required_successes
while True:
prob = nbincdf(n, required_successes, response_rate)
if prob >= conf:
print(f"{conf:.0%}置信水平: 需要筛选{n}名患者")
break
n += 1
#80%置信水平: 需要筛选30名患者
#90%置信水平: 需要筛选35名患者
#95%置信水平: 需要筛选39名患者
示例4. 生态学 - 物种观察研究
在生态学中，负二项分布用于建模观察到第R个稀有物种个体前需要的调查次数
场景：稀有物种的观察概率为2%，需要观察到5个个体
observation_prob = 0.02
target_observations = 5
计算在不同调查次数下的成功概率
survey_counts = [100, 200, 300, 400, 500]
print(f"观察到{target_observations}个个体的累积概率:")
for n in survey_counts:
prob = nbincdf(n, target_observations, observation_prob)
print(f"进行{n}次调查: {prob:.1%}")
#观察到5个个体的累积概率:
#进行100次调查: 6.0%
#进行200次调查: 39.1%
#进行300次调查: 73.1%
#进行400次调查: 90.8%
#进行500次调查: 97.4%
资源优化分析
budget_constraints = [200, 300, 400]
预算约束下的观察策略:
for budget in budget_constraints:
best_r = 1
best_prob = 0
for r in range(1, 10):
prob = nbincdf(budget, r, observation_prob)
if prob > best_prob:
best_prob = prob
best_r = r
print(f"预算{budget}次调查: 建议目标{best_r}个个体 (成功概率{best_prob:.1%})")
#预算200次调查: 建议目标1个个体 (成功概率98.3%)
#预算300次调查: 建议目标1个个体 (成功概率99.8%)
#预算400次调查: 建议目标1个个体 (成功概率100.0%)
示例5. 网络可靠性 - 故障分析
在网络工程中，负二项分布用于建模第R次网络故障前运行的时间段数
场景：每周网络故障概率为5%，分析第3次故障前的运行时间
failure_prob = 0.05
target_failures = 3
计算在不同时间周期内的累积概率
time_periods = [20, 40, 60, 80, 100] # 周数
print(f"第{target_failures}次故障前的运行时间概率:")
for weeks in time_periods:
prob = nbincdf(weeks, target_failures, failure_prob)
print(f"运行{weeks}周: {prob:.1%}概率发生不超过{target_failures}次故障")
#第3次故障前的运行时间概率:
#运行20周: 10.5%概率发生不超过3次故障
#运行40周: 36.5%概率发生不超过3次故障
#运行60周: 61.6%概率发生不超过3次故障
#运行80周: 79.1%概率发生不超过3次故障
#运行100周: 89.4%概率发生不超过3次故障
维护策略优化
maintenance_intervals = [26, 52, 104] # 半年、1年、2年
不同维护周期下的风险分析:
for interval in maintenance_intervals:
prob_1_failure = nbincdf(interval, 1, failure_prob)
prob_2_failures = nbincdf(interval, 2, failure_prob)
prob_3_failures = nbincdf(interval, 3, failure_prob)
print(f"\n{interval}周({interval // 52}年)维护周期:")
print(f" 至少1次故障: {1 - prob_1_failure:.1%}")
print(f" 至少2次故障: {1 - prob_2_failures:.1%}")
print(f" 至少3次故障: {1 - prob_3_failures:.1%}")
#26周(0年)维护周期:
# 至少1次故障: 25.0%
# 至少2次故障: 58.8%
# 至少3次故障: 82.5%
#52周(1年)维护周期:
# 至少1次故障: 6.6%
# 至少2次故障: 24.1%
# 至少3次故障: 47.7%
#104周(2年)维护周期:
# 至少1次故障: 0.5%
# 至少2次故障: 2.9%
# 至少3次故障: 9.2%
示例6. 客户服务 - 投诉处理
在客户服务中，负二项分布用于建模收到第R个严重投诉前处理的客户交互次数
场景：严重投诉率为1%，分析服务质量
complaint_rate = 0.01
critical_complaints = 2
不同交互数量下的风险分析
interaction_levels = [100, 200, 500, 1000]
print(f"收到第{critical_complaints}个严重投诉前的累积概率:")
for interactions in interaction_levels:
prob = nbincdf(interactions, critical_complaints, complaint_rate)
risk = 1 - prob
print(f"处理{interactions}次交互: 风险率{risk:.2%}")
#收到第2个严重投诉前的累积概率:
#处理100次交互: 风险率72.84%
#处理200次交互: 风险率39.92%
#处理500次交互: 风险率3.91%
#处理1000次交互: 风险率0.05%
服务质量标准设定
risk_tolerance = 0.05 # 5%风险容忍度
print(f"\n达到{risk_tolerance:.0%}风险容忍度的服务质量标准:")
for r in range(1, 4):
interactions = 1
while True:
prob = nbincdf(interactions, r, complaint_rate)
risk = 1 - prob
if risk <= risk_tolerance:
print(f"允许{r}个严重投诉: 每{interactions}次交互")
break
interactions += 1
#达到5%风险容忍度的服务质量标准:
#允许1个严重投诉: 每298次交互
#允许2个严重投诉: 每471次交互
#允许3个严重投诉: 每625次交互
示例7. 制造业 - 设备维护
在制造业中，负二项分布用于建模设备第R次故障前的生产批次数量
场景：设备每批次的故障概率为3%，制定维护计划
batch_failure_prob = 0.03
maintenance_trigger = 3 # 第3次故障时进行维护
不同生产批次下的维护概率
batch_counts = [50, 100, 150, 200]
print(f"第{maintenance_trigger}次故障前的生产批次分析:")
for batches in batch_counts:
prob = nbincdf(batches, maintenance_trigger, batch_failure_prob)
print(f"生产{batches}批次: {prob:.1%}概率需要维护")
#第3次故障前的生产批次分析:
#生产50批次: 21.2%概率需要维护
#生产100批次: 60.0%概率需要维护
#生产150批次: 84.0%概率需要维护
#生产200批次: 94.5%概率需要维护
预防性维护策略比较:
preventive_intervals = [50, 100, 150]
for interval in preventive_intervals:
# 预防性维护：在固定间隔维护，不考虑故障次数
preventive_cost = interval * 100 # 假设每批次生产成本
# corrective维护：基于故障的维护
expected_batches = 0
for batches in range(1, 1000):
prob = nbincdf(batches, maintenance_trigger, batch_failure_prob)
if prob >= 0.5: # 中位数
expected_batches = batches
break
corrective_cost = expected_batches * 100 + 5000 # 额外维护成本
print(f"{interval}批次预防性维护: 成本{preventive_cost} vs 纠正性维护: 成本{corrective_cost}")
#50批次预防性维护: 成本5000 vs 纠正性维护: 成本13600
#100批次预防性维护: 成本10000 vs 纠正性维护: 成本13600
#150批次预防性维护: 成本15000 vs 纠正性维护: 成本13600
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import NegativeBinomial, cdf
from sympy import Expr
from scipy.stats import nbinom
def negative_binomial_cdf(input_str):
"""
对标 MATLAB 的 nbincdf 函数，计算负二项分布的累积分布函数
参数:
input_str (str): 输入表达式，支持三种形式：
- 标量模式："(2, 3, 0.5)" 表示 x=2, R=3, p=0.5
- 矩阵模式："([[1,2],[3,4]], 5, 0.2)" 表示 x矩阵, R=5, p=0.2
返回:
Union[float, Matrix, str]: 计算结果或错误信息
MATLAB 参数对照：
nbincdf(x, R, p) 其中：
- x: 失败次数（标量或矩阵）
- R: 要求成功的次数（正整数）
- p: 单次试验成功概率（0 ≤ p ≤ 1）
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def nbincdf_sci(k, n, p):
"""
计算负二项分布的累积分布函数（CDF）
对标 MATLAB 的 nbincdf 函数
参数:
k (int/float/array): 成功次数
n (int): 需要的失败次数（正整数）
p (float): 每次试验成功的概率 (0 ≤ p ≤ 1)
返回:
float/array: 负二项分布的CDF值
"""
# 参数验证
if not isinstance(n, int) or n <= 0:
raise ValueError(f"参数 n 必须为正整数，当前 n={n}")
if not (0 <= p <= 1):
raise ValueError(f"概率 p 必须在 [0,1] 之间，当前 p={p}")
# 计算CDF（注意scipy的n参数表示成功次数，与MATLAB相反）
return nbinom.cdf(k, n, p)
# 计算逻辑 ======================================================
def nbincdf_sym(x_val, R, p):
# 验证 R 是正整数
if not (R.is_integer and R > 0):
return f"参数 R 必须为正整数，当前 R={R}"
# 验证 p 在概率范围内
if not (0 <= p <= 1):
return f"概率 p 必须在 [0,1] 之间，当前 p={p.evalf(4)}"
# 分布定义 ======================================================
X = NegativeBinomial("X", R, p) # 定义分布：X ~ NegativeBinomial(R, p)
"""计算单个值的 CDF，自动处理符号和数值"""
cdf_expr = cdf(X)(x_val)
# 如果是符号表达式直接返回，数值则求值
if isinstance(cdf_expr, Expr):
return cdf_expr
else:
return cdf_expr.evalf()
# 标量处理
# 参数解析 ======================================================
if isinstance(expr, tuple) and len(expr) == 3:
if all(e.is_number for e in expr):
x, R, p = float(expr[0]), int(expr[1]), float(expr[2])
result = nbincdf_sci(x, R, p)
else:
result = nbincdf_sym(*expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1：标量计算（数值）
print(negative_binomial_cdf("(2, 3, 0.5)"))
# 0.5
# 示例3：符号计算
x = sp.symbols('x')
print(negative_binomial_cdf(f"({x}, 2, 0.5)"))
# Piecewise((-0.25*0.5**(floor(x) + 1)*(floor(x) + 2)*hyper((1, 1.0*floor(x) + 3.0),
# (floor(x) + 2,), 0.5) + 1.0, x >= 0), (0, True))
负二项式参数估计
[parmhat，parmci]=nbinfit（data，alpha）返回MLE和100（1-alpha）百分比置信区间。默认情况下，alpha=0.05，对应于95%的置信区间。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.stats as stats
import scipy.optimize as optimize
import ast
def parse_input(input_str):
"""
解析输入字符串并转换为数值数组。
参数:
input_str (str): 输入的字符串，可以是列表、元组或SymPy矩阵的字符串形式。
返回:
numpy.ndarray: 转换后的数值数组。若解析失败则返回None。
"""
try:
# 尝试解析为Python字面量（列表、元组或数值）
data = ast.literal_eval(input_str)
if isinstance(data, (list, tuple)):
return np.array(data, dtype=float).flatten()
else:
return np.array([float(data)], dtype=float)
except:
try:
# 尝试作为SymPy表达式解析
expr = sp.sympify(input_str)
if isinstance(expr, sp.MatrixBase):
# 转换SymPy矩阵为numpy数组并展平
return np.array(expr.tolist(), dtype=float).flatten()
elif isinstance(expr, (sp.List, list, tuple)):
# 转换SymPy列表或类似结构为数组
return np.array([float(e) for e in expr], dtype=float)
else:
# 处理单个数值的情况
return np.array([float(expr)], dtype=float)
except:
return None
def negative_binomial_parameter_estimates(input_str):
"""
估计负二项分布的参数r和p，对标Matlab的nbinfit函数。
参数:
input_str (str): 输入数据的字符串形式，如"[3,5,7,2,4]"或"Matrix([1,2,3])"。
返回:
tuple or str: 成功时返回(r, p)的估计值，失败时返回错误信息。
"""
try:
data = parse_input(input_str)
if data is None:
return "错误：无法解析输入数据"
if len(data) == 0:
return "错误：输入数据为空"
# 计算矩估计作为初始猜测
mean = np.mean(data)
var = np.var(data, ddof=0) # 使用样本方差（非无偏估计）
if var > mean:
p_initial = mean / var
r_initial = mean ** 2 / (var - mean)
initial_guess = [r_initial, p_initial]
else:
# 方差不足时使用默认初始值
initial_guess = [1.0, 0.5]
# 定义负对数似然函数
def negative_binomial_nll(params):
r, p = params
# 处理参数超出范围的情况
if r <= 1e-6 or p <= 1e-6 or p >= 1 - 1e-6:
return np.inf
# 计算对数似然，避免无效值
try:
log_pmf = stats.nbinom.logpmf(data, r, p)
return -np.sum(log_pmf)
except:
return np.inf
# 参数约束：r > 0，0 < p < 1
bounds = [(1e-6, None), (1e-6, 1 - 1e-6)]
# 使用L-BFGS-B方法进行优化
res = optimize.minimize(
negative_binomial_nll,
initial_guess,
method='L-BFGS-B',
bounds=bounds
)
if res.success:
r_hat, p_hat = res.x
return (r_hat, p_hat)
else:
return f"优化失败：{res.message}"
except Exception as e:
return f"错误：{str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例输入数据（字符串形式）
test_cases = [
"[3,5,7,2,4]",
# 估计参数: r=202.9085, p=0.9798
"Matrix([3,5,7,2,4])",
# 估计参数: r=202.9085, p=0.9798
"5,3,7,8,2",
# 估计参数: r=125.0000, p=0.9615
]
for input_str in test_cases:
print(f"输入: {input_str}")
result = negative_binomial_parameter_estimates(input_str)
if isinstance(result, tuple):
r, p = result
print(f"估计参数: r={r:.4f}, p={p:.4f}\n")
else:
print(f"结果: {result}\n")
负二项式逆累积分布函数
X=nbininv(Y,R,P)返回负二项式cdf的倒数,其中包含相应的成功次数,R和单次试验中的成功概率P.
由于二项式分布是离散的,nbininv返回最小整数X,使得在X处计算的负二项cdf等于或超过Y. Y,R和P可以是向量或矩阵, 它们都具有相同的大小, 也就是X的大小.
将Y,R或P的标量输入扩展为与其他输入具有相同维度的常数数组.
示例1. 制造业 - 质量检测计划制定
制定质量检测计划：确定在给定置信水平下需要检查的产品数量
场景：缺陷率为5%，要找到第3个缺陷产品，置信水平为90%
confidence_level = 0.90
r = 3 # 目标缺陷产品数
p = 0.05 # 缺陷率
required_inspections = nbininv(confidence_level, r, p)
total_products = required_inspections + r # 总检查产品数 = 失败次数 + 成功次数
print(f"目标: 找到{r}个缺陷产品 (缺陷率{p:.1%})")
print(f"置信水平: {confidence_level:.0%}")
print(f"需要检查的产品数量: {total_products:.0f}")
print(f"其中: {required_inspections:.0f}个合格产品 + {r}个缺陷产品")
#目标: 找到3个缺陷产品 (缺陷率5.0%)
#置信水平: 90%
#需要检查的产品数量: 105
#其中: 102个合格产品 + 3个缺陷产品
不同置信水平下的比较
confidence_levels = [0.80, 0.90, 0.95, 0.99]
不同置信水平下的需求:
for conf in confidence_levels:
insp = nbininv(conf, r, p)
total = insp + r
print(f"{conf:.0%}置信水平: {total:.0f}个产品")
#80%置信水平: 85个产品
#90%置信水平: 105个产品
#95%置信水平: 124个产品
#99%置信水平: 165个产品
示例2. 销售管理 - 销售目标规划
规划达成销售目标所需的客户接触次数
场景：转化率8%，要获得5个订单，置信水平85%
conversion_rate = 0.08
target_orders = 5
confidence = 0.85
required_failures = nbininv(confidence, target_orders, conversion_rate)
total_contacts = required_failures + target_orders
print(f"销售目标: {target_orders}个订单")
print(f"转化率: {conversion_rate:.1%}")
print(f"置信水平: {confidence:.0%}")
print(f"需要接触的客户数: {total_contacts:.0f}")
print(f"其中: {required_failures:.0f}次未转化 + {target_orders}次转化")
#销售目标: 5个订单
#转化率: 8.0%
#置信水平: 85%
#需要接触的客户数: 90
#其中: 85次未转化 + 5次转化
敏感性分析：不同转化率下的需求
conversion_rates = [0.05, 0.08, 0.10, 0.12]
不同转化率下的客户接触需求:
for rate in conversion_rates:
failures = nbininv(confidence, target_orders, rate)
contacts = failures + target_orders
print(f"{rate:.1%}转化率: {contacts:.0f}次接触")
#5.0%转化率: 144次接触
#8.0%转化率: 90次接触
#10.0%转化率: 72次接触
#12.0%转化率: 59次接触
示例3. 医疗研究 - 临床试验样本量确定
确定临床试验中需要筛选的患者数量
场景：治疗响应率25%，需要10个有效病例，置信水平90%
response_rate = 0.25
required_successes = 10
confidence = 0.90
screening_failures = nbininv(confidence, required_successes, response_rate)
total_patients = screening_failures + required_successes
print(f"研究目标: {required_successes}个有效病例")
print(f"预期响应率: {response_rate:.1%}")
print(f"置信水平: {confidence:.0%}")
print(f"需要筛选的患者数: {total_patients:.0f}")
print(f"其中: {screening_failures:.0f}例无效 + {required_successes}例有效")
#研究目标: 10个有效病例
#预期响应率: 25.0%
#置信水平: 90%
#需要筛选的患者数: 55
#其中: 45例无效 + 10例有效
预算约束分析
budget_constraints = [100, 200, 300, 400]
预算约束下的可行性分析:
for budget in budget_constraints:
max_successes = 1
while True:
failures = nbininv(confidence, max_successes, response_rate)
if (failures + max_successes) <= budget:
max_successes += 1
else:
break
print(f"预算{budget}患者: 最多可获得{max_successes - 1}个有效病例")
#预算100患者: 最多可获得20个有效病例
#预算200患者: 最多可获得42个有效病例
#预算300患者: 最多可获得65个有效病例
#预算400患者: 最多可获得89个有效病例
示例4. 生态研究 - 野外调查规划
规划野外调查中需要的观察次数
场景：稀有物种观察概率2%，要观察到5次，置信水平95%
observation_prob = 0.02
target_observations = 5
confidence = 0.95
unsuccessful_surveys = nbininv(confidence, target_observations, observation_prob)
total_surveys = unsuccessful_surveys + target_observations
print(f"调查目标: 观察到{target_observations}次目标物种")
print(f"单次观察概率: {observation_prob:.1%}")
print(f"置信水平: {confidence:.0%}")
print(f"需要进行的调查次数: {total_surveys:.0f}")
print(f"其中: {unsuccessful_surveys:.0f}次未观察到 + {target_observations}次成功观察")
#调查目标: 观察到5次目标物种
#单次观察概率: 2.0%
#置信水平: 95%
#需要进行的调查次数: 456
#其中: 451次未观察到 + 5次成功观察
资源优化策略
if total_surveys > 500:
print("🔴 调查成本过高，建议:")
print(" - 增加单次调查效率")
print(" - 降低置信水平要求")
print(" - 调整观察目标")
else:
print("🟢 调查计划在可行范围内")
#🟢 调查计划在可行范围内
示例5. 网络运维 - 故障预测与容量规划
预测网络故障发生前的运行时间
场景：每周故障概率3%，预测第3次故障发生的时间
failure_prob = 0.03
target_failures = 3
confidence = 0.8
failure_free_weeks = nbininv(confidence, target_failures, failure_prob)
total_weeks = failure_free_weeks + target_failures
print(f"预测目标: 第{target_failures}次网络故障")
print(f"每周故障概率: {failure_prob:.1%}")
print(f"置信水平: {confidence:.0%}")
print(f"预测运行时间: {total_weeks:.0f}周")
print(f"其中: {failure_free_weeks:.0f}周无故障 + {target_failures}周有故障")
#预测目标: 第3次网络故障
#每周故障概率: 3.0%
#置信水平: 80%
#预测运行时间: 142周
#其中: 139周无故障 + 3周有故障
维护策略建议
maintenance_interval = 52 # 1年
if total_weeks <= maintenance_interval:
print(f"⚠️ 警告: 预计{total_weeks}周内会发生{target_failures}次故障")
print(" 建议缩短维护周期或改进系统可靠性")
else:
print(f"✅ 系统可靠性满足年度维护计划")
#✅ 系统可靠性满足年度维护计划
示例6. 客户服务 - 服务质量标准制定
制定基于风险容忍度的服务质量标准
场景：严重投诉率1%，制定服务交互次数标准
complaint_rate = 0.01
risk_tolerance = 0.05 # 5%风险容忍度
confidence = 1 - risk_tolerance
针对不同投诉次数的标准
complaint_thresholds = [1, 2, 3]
print(f"服务质量标准 (投诉率{complaint_rate:.1%}, 风险容忍度{risk_tolerance:.0%})")
print("=" * 50)
for complaints in complaint_thresholds:
safe_interactions = nbininv(confidence, complaints, complaint_rate)
total_interactions = safe_interactions + complaints
print(f"\n允许{complaints}次严重投诉:")
print(f" - 安全交互次数: {safe_interactions:.0f}")
print(f" - 总交互次数标准: {total_interactions:.0f}")
print(f" - 实际风险: {risk_tolerance:.1%}")
#服务质量标准 (投诉率1.0%, 风险容忍度5%)
#==================================================
#允许1次严重投诉:
# - 安全交互次数: 298
# - 总交互次数标准: 299
# - 实际风险: 5.0%
#允许2次严重投诉:
# - 安全交互次数: 471
# - 总交互次数标准: 473
# - 实际风险: 5.0%
#允许3次严重投诉:
# - 安全交互次数: 625
# - 总交互次数标准: 628
# - 实际风险: 5.0%
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import scipy.stats as stats
def negative_binomial_inverse_cdf(input_str):
"""
计算负二项分布的逆累积分布函数，支持标量和矩阵输入。
参数:
input_str: 输入字符串，格式为三元组 (Y, R, P)，例如 "([0.5], 5, 0.4)"
返回:
计算结果，SymPy矩阵或标量。错误时返回字符串描述。
"""
try:
# 解析输入字符串为Python对象
expr = sp.sympify(input_str)
if not isinstance(expr, (tuple, list)) or len(expr) != 3:
return f"输入错误: 需要三元组 (Y, R, P)，但得到的是 {input_str}"
def eval_element(y_val, r_val, p_val):
"""计算单个元素的逆CDF值"""
try:
y_val = float(y_val)
r_val = float(r_val)
p_val = float(p_val)
if not (0 <= y_val <= 1):
return np.nan
if not (0 < p_val < 1):
return np.nan
if r_val <= 0:
return np.nan
x = stats.nbinom.ppf(y_val, r_val, p_val)
return x
except Exception as e:
print(f"元素计算错误: {e}")
return np.nan
# 标量情况
result = eval_element(*expr)
return sp.sympify(result) if not np.isnan(result) else "无效输入值"
except Exception as e:
return f"错误: {str(e)}"
# 示例用法
if __name__ == "__main__":
# 标量输入
print(negative_binomial_inverse_cdf("(0.99, 10, 0.5)"))
# 23.0000000000000
负二项式概率密度函数
Y=nbinpdf(X,R,P)使用相应的成功次数,R和单次试验的成功概率P,在X中的每个值处返回负二项式pdf.
X,R和P可以是向量或矩阵,它们都具有相同的大小,也就是Y的大小.X,R或P的标量输入被扩展为与其他输入具有相同维度的常数数组.
示例1: 市场营销 - 客户转化分析
假设：需要获得3个成功转化，每次联系客户的转化概率为0.2")
k_values = list(range(11)) # 失败次数（未转化的客户数）
r = 3 # 需要的成功转化数
p = 0.2 # 单次转化概率
for k in k_values:
prob = nbinpdf(k, r, p)
print(f"在获得第{r}个转化前联系{k}个未转化客户的概率: {prob:.4f}")
#在获得第3个转化前联系0个未转化客户的概率: 0.0080
#在获得第3个转化前联系1个未转化客户的概率: 0.0192
#在获得第3个转化前联系2个未转化客户的概率: 0.0307
#在获得第3个转化前联系3个未转化客户的概率: 0.0410
#在获得第3个转化前联系4个未转化客户的概率: 0.0492
#在获得第3个转化前联系5个未转化客户的概率: 0.0551
#在获得第3个转化前联系6个未转化客户的概率: 0.0587
#在获得第3个转化前联系7个未转化客户的概率: 0.0604
#在获得第3个转化前联系8个未转化客户的概率: 0.0604
#在获得第3个转化前联系9个未转化客户的概率: 0.0591
#在获得第3个转化前联系10个未转化客户的概率: 0.0567
示例2: 质量控制 - 产品缺陷检测
假设：需要找到5个合格产品，每个产品合格概率为0.9
k_values = list(range(9)) # 缺陷产品数
r = 5 # 需要的合格产品数
p = 0.9 # 产品合格率
for k in k_values:
prob = nbinpdf(k, r, p)
print(f"找到第{r}个合格产品前发现{k}个缺陷产品的概率: {prob:.4f}")
#找到第5个合格产品前发现0个缺陷产品的概率: 0.5905
#找到第5个合格产品前发现1个缺陷产品的概率: 0.2952
#找到第5个合格产品前发现2个缺陷产品的概率: 0.0886
#找到第5个合格产品前发现3个缺陷产品的概率: 0.0207
#找到第5个合格产品前发现4个缺陷产品的概率: 0.0041
#找到第5个合格产品前发现5个缺陷产品的概率: 0.0007
#找到第5个合格产品前发现6个缺陷产品的概率: 0.0001
#找到第5个合格产品前发现7个缺陷产品的概率: 0.0000
#找到第5个合格产品前发现8个缺陷产品的概率: 0.0000
示例3: 医学研究 - 药物有效性试验
假设：需要观察到8个有效病例，每个患者有效的概率为0.3
k_values = [0, 5, 10, 15, 20, 25] # 无效病例数
r = 8 # 需要的有效病例数
p = 0.3 # 治疗有效性概率
for k in k_values:
prob = nbinpdf(k, r, p)
print(f"观察到第{r}个有效病例前治疗{k}个无效病例的概率: {prob:.6f}")
#观察到第8个有效病例前治疗0个无效病例的概率: 0.000066
#观察到第8个有效病例前治疗5个无效病例的概率: 0.008733
#观察到第8个有效病例前治疗10个无效病例的概率: 0.036043
#观察到第8个有效病例前治疗15个无效病例的概率: 0.053122
#观察到第8个有效病例前治疗20个无效病例的概率: 0.046490
#观察到第8个有效病例前治疗25个无效病例的概率: 0.029615
示例4: 销售预测 - 完成销售目标
不同销售能力下完成5笔销售前的失败次数分布：
r = 5 # 销售目标
probabilities = [0.1, 0.3, 0.5] # 不同销售能力的成功率
k_values = list(range(16)) # 失败次数
for p in probabilities:
print(f"\n销售成功率 {p * 100}%:")
for k in [0, 5, 10, 15]:
prob = nbinpdf(k, r, p)
print(f" 完成第{r}笔销售前失败{k}次的概率: {prob:.6f}")
#销售成功率 10.0%:
# 完成第5笔销售前失败0次的概率: 0.000010
# 完成第5笔销售前失败5次的概率: 0.000744
# 完成第5笔销售前失败10次的概率: 0.003490
# 完成第5笔销售前失败15次的概率: 0.007980
#销售成功率 30.0%:
# 完成第5笔销售前失败0次的概率: 0.002430
# 完成第5笔销售前失败5次的概率: 0.051460
# 完成第5笔销售前失败10次的概率: 0.068710
# 完成第5笔销售前失败15次的概率: 0.044716
#销售成功率 50.0%:
# 完成第5笔销售前失败0次的概率: 0.031250
# 完成第5笔销售前失败5次的概率: 0.123047
# 完成第5笔销售前失败10次的概率: 0.030548
# 完成第5笔销售前失败15次的概率: 0.003696
示例5: 向量输入 - 批量计算
同时计算多个(k, r, p)组合的概率：
多个失败次数，相同的r和p
k_vector = [0, 1, 2, 3, 4, 5]
r_val = 2
p_val = 0.4
result = nbinpdf(k_vector, r_val, p_val)
print(f"k={k_vector}, r={r_val}, p={p_val}")
print(f"概率分布: {result}")
#k=[0, 1, 2, 3, 4, 5], r=2, p=0.4
#概率分布: [0.16,0.192,0.1728,0.13824,0.10368,0.0746496]
示例6: 矩阵输入 - 多维数据分析
使用矩阵输入进行批量计算：
创建矩阵输入
k_matrix = [[0, 1, 2], [3, 4, 5]] #失败次数
r_matrix = [[2, 2, 2], [3, 3, 3]] #成功次数
p_matrix = [[0.3, 0.5, 0.7], [0.3, 0.5, 0.7]] #成功概率
result_matrix = nbinpdf(k_matrix, r_matrix, p_matrix)
对应的概率密度矩阵:
print(result_matrix)
#[[0.09,0.25,0.1323]
[0.09261,0.1171875,0.01750329]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import NegativeBinomial, density
def negative_binomial_pdf(input_str):
"""
计算负二项分布的概率密度函数（PDF）。
参数:
input_str: 输入字符串，格式为三元组(k, r, p)，其中：
k: 失败次数（标量、向量或矩阵）
r: 成功次数（标量、向量或矩阵）
p: 成功概率（标量、向量或矩阵）
返回:
PDF值组成的SymPy矩阵或标量。错误时返回错误信息字符串。
"""
try:
# 解析输入字符串为Python对象
expr = sp.sympify(input_str)
error = False
result = None
# 计算单个元素的PDF
def eval_element(k_val, r_val, p_val):
"""计算负二项分布PDF值"""
if r_val > 0 and (0 < p_val < 1):
X = NegativeBinomial('X', r_val, p_val)
pdf = density(X)(k_val)
return pdf.evalf() # 转换为数值形式
else:
sp.nan
if isinstance(expr, tuple) and len(expr) == 3:
result = eval_element(*expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误：{str(e)}"
# 示例用法
if __name__ == "__main__":
# 标量输入
print("标量输入:", negative_binomial_pdf("(2, 3, 0.5)"))
# 0.187500000000000
负二项式随机数
RND=nbinnd(R,P)是从负二项分布中选择的随机数矩阵,具有相应的成功数R和单次试验的成功概率P.
R和P可以是大小相同的向量或矩阵,也就是RND的大小.R或P的标量输入被扩展为与其他输入具有相同维数的常数数组.
示例1. 客户服务场景 - 电话呼叫中心
模拟客服在接到第5个投诉电话之前，成功处理的满意电话数量
print("客服场景:", nbinnd(5, 0.7))
#客服场景: 3
#意义：在接到第5个投诉前，预计能成功处理多少个客户问题
多个客服人员的情况
print("多个客服:", nbinnd([3,4,5], [0.6,0.7,0.8]))
#多个客服: [0 1 0]
示例2. 市场营销 - 客户转化
在获得第10个付费客户前，需要进行的营销接触次数
print("营销转化:", nbinnd(10, 0.05))
#营销转化: 215
#意义：获得10个付费客户需要接触多少潜在客户
不同产品线的营销效果比较
print("多产品线:", nbinnd([8,12,15], [0.03,0.04,0.06]))
#多产品线: [154 153 192]
示例3. 医疗研究 - 药物试验
在观察到第8个有效病例前，需要试验的患者数量
print("药物试验:", nbinnd(8, 0.15))
#药物试验: 24
#意义：需要多少患者才能观察到8个有效病例
多个研究中心的数据
print("多中心试验:", nbinnd(8, 0.15, (3, 5)))
#多中心试验:
#[[30 55 66 31 48]
[63 55 55 48 64]
[45 55 41 46 71]]
#生成3×5矩阵，模拟3个中心的5次试验
示例4. 销售预测 - 销售团队业绩
销售人员在达成第20笔交易前，需要进行的客户拜访次数
print("销售业绩:", nbinnd(20, 0.1))
#销售业绩: 137
#意义：达成20笔交易需要拜访多少客户
整个销售团队的表现
print("销售团队:", nbinnd(20, 0.1, 8))
#销售团队: [198 163 138 172 186 169 171 150]
#模拟8个销售人员的表现
示例5. 教育评估 - 考试通过率
学生在第3次不及格前，通过的考试科目数量
print("学业表现:", nbinnd(3, 0.85))
#学业表现: 0
#意义：在出现3次不及格前，预计能通过多少科目
多个班级的比较
print("多班级数据:", nbinnd([3,3,3], [0.8,0.85,0.9]))
#多班级数据: [1 0 0]
示例6. 分析销售团队绩效
参数：r=成功交易数，p=转化率
sales_data = nbinnd(15, 0.08, 50)
print(f"销售团队需要客户接触次数: {sales_data}")
print(f"平均需要接触次数: {np.mean(sales_data):.1f}")
print(f"最低接触次数: {np.min(sales_data)}")
print(f"最高接触次数: {np.max(sales_data)}")
#销售团队需要客户接触次数:
#[126 133 196 173 175 219 162 140 154 191 227 174 174 141 108 186 89 157
110 129 156 206 188 134 168 186 183 208 222 126 138 117 257 178 233 159
209 200 177 192 186 133 272 127 206 211 172 176 183 172]
#平均需要接触次数: 172.8
#最低接触次数: 89
#最高接触次数: 272
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import scipy.stats as stats
def negative_binomial_random_numbers(input_str):
"""
生成负二项分布的随机数，对标MATLAB的nbinrnd函数。
参数:
input_str: 输入字符串，格式为 (r, p) 或 (r, p, size)，
其中r为成功次数，p为成功概率，size为输出维度。
返回:
SymPy矩阵或数值。错误时返回错误信息字符串。
"""
try:
# 安全解析输入字符串
expr = sp.sympify(input_str)
error = False
result = None
def eval_random_number(r, p):
if r > 0 and (0 < p < 1):
random_number = stats.nbinom.rvs(float(r), float(p))
# 计算随机数
return random_number
else:
return sp.nan
if isinstance(expr, tuple):
# 解析参数
size = None
if len(expr) == 2:
result = eval_random_number(*expr)
elif len(expr) == 3:
n, x, size = expr
if (not size.is_integer) or size <= 0:
return "错误：size必须为正整数"
result = stats.nbinom.rvs(float(n), float(x), size=int(size))
result = sp.Matrix(result).T
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误：{str(e)}"
# 示例用法
if __name__ == "__main__":
# 标量输入
print("标量输出:", negative_binomial_random_numbers("(5, 0.4)"))
# 7
# 带size参数的向量输出
print("向量输出:", negative_binomial_random_numbers("(3, 0.5, 5)"))
# Matrix([[1, 11, 3, 4, 1]])
负二项均值和方差
M = nbinstat(R,P)返回负二项分布的平均值,其中R成功,单次试验成功的概率为P.
[M,V]=nbinstat(R,P)也返回分布的方差.
示例1.客户服务场景 - 呼叫中心预测
在接到第5个投诉前，预计成功处理的电话数量及其波动范围
mean, variance = nbinstat(5, 0.7)
print(f"客户服务预测: 均值={mean:.2f}, 方差={variance:.2f}")
#客户服务预测: 均值=2.14, 方差=3.06
#意义：平均在2.14个成功电话后会遇到第5个投诉，波动范围约±1.75个电话
示例2. 生产线质量控制分析
分析不同次品率下的生产表现
scenarios = [
(3, 0.95), # 高质量生产线：次品率5%
(3, 0.90), # 中等质量：次品率10%
(3, 0.80) # 低质量：次品率20%
]
for r, p in scenarios:
mean, var = nbinstat(r, p)
std_dev = np.sqrt(var)
print(f"次品率{1 - p:.0%}: 均值={mean:.1f}个合格品, 标准差={std_dev:.1f}")
#次品率5%: 均值=0.2个合格品, 标准差=0.4
#次品率10%: 均值=0.3个合格品, 标准差=0.6
#次品率20%: 均值=0.7个合格品, 标准差=1.0
示例3. 市场营销投入分析
分析不同转化率下获得10个客户的营销成本
conversion_rates = [0.02, 0.05, 0.08, 0.10]
获得10个客户需要的营销接触次数分析:
for p in conversion_rates:
mean, var = nbinstat(10, p)
std_dev = np.sqrt(var)
print(f"转化率{p:.1%}: 平均需要{mean:.0f}次接触 (±{std_dev:.0f}次)")
#转化率2.0%: 平均需要490次接触 (±157次)
#转化率5.0%: 平均需要190次接触 (±62次)
#转化率8.0%: 平均需要115次接触 (±38次)
#转化率10.0%: 平均需要90次接触 (±30次)
示例4. 医疗研究样本量规划
药物试验：观察8个有效病例前的患者数量预测
effectiveness_rates = [0.1, 0.15, 0.2, 0.25]
药物试验样本量预测 (需要8个有效病例):
for p in effectiveness_rates:
mean, var = nbinstat(8, p)
confidence_interval = 1.96 * np.sqrt(var) # 95%置信区间
print(f"有效率{p:.0%}: 平均需要{mean:.0f}名患者, 95%CI: ±{confidence_interval:.0f}名")
#有效率10%: 平均需要72名患者, 95%CI: ±53名
#有效率15%: 平均需要45名患者, 95%CI: ±34名
#有效率20%: 平均需要32名患者, 95%CI: ±25名
#有效率25%: 平均需要24名患者, 95%CI: ±19名
示例5. 网络可靠性工程
服务器在5次崩溃前的正常运行时间分析
reliability_levels = [0.97, 0.98, 0.99, 0.995]
服务器可靠性分析 (5次崩溃前):
for p in reliability_levels:
mean, var = nbinstat(5, p)
print(f"可靠性{p:.1%}: 平均运行{mean:.0f}小时, 方差={var:.0f}")
#可靠性97.0%: 平均运行0小时, 方差=0
#可靠性98.0%: 平均运行0小时, 方差=0
#可靠性99.0%: 平均运行0小时, 方差=0
#可靠性99.5%: 平均运行0小时, 方差=0
示例6. 销售团队绩效评估
比较不同销售人员的转化效率
sales_team = [
("新手销售", 15, 0.05),
("资深销售", 15, 0.08),
("精英销售", 15, 0.12)
]
销售团队绩效分析 (达成15笔交易):
for name, r, p in sales_team:
mean, var = nbinstat(r, p)
efficiency = r / mean if mean > 0 else 0
print(f"{name}: 平均需要{mean:.0f}次拜访, 转化效率{efficiency:.2%}")
#新手销售: 平均需要285次拜访, 转化效率5.26%
#资深销售: 平均需要172次拜访, 转化效率8.70%
#精英销售: 平均需要110次拜访, 转化效率13.64%
示例7. 多维度业务分析
不同产品线在不同目标下的表现
success_targets = [[5, 10], [8, 15]] # 二维数组：行代表产品类型，列代表成功目标
conversion_rates = [[0.1, 0.1], [0.15, 0.15]] # 对应的转化率
多产品线业务分析:
for i in range(len(success_targets)):
for j in range(len(success_targets[i])):
r = success_targets[i][j]
p = conversion_rates[i][j]
mean, var = nbinstat(r, p)
print(f"产品{i + 1}-目标{r}: 均值={mean:.1f}, 标准差={np.sqrt(var):.1f}")
#产品1-目标5: 均值=45.0, 标准差=21.2
#产品1-目标10: 均值=90.0, 标准差=30.0
#产品2-目标8: 均值=45.3, 标准差=17.4
#产品2-目标15: 均值=85.0, 标准差=23.8
示例8. 教育资源分配规划
教育机构：学生在第3次不及格前通过的科目数预测
pass_rates = [0.7, 0.8, 0.9]
学生学业表现预测 (3次不及格前):
for p in pass_rates:
mean, var = nbinstat(3, p)
print(f"通过率{p:.0%}: 平均通过{mean:.1f}门课, 波动范围±{np.sqrt(var):.1f}门")
#通过率70%: 平均通过1.3门课, 波动范围±1.4门
#通过率80%: 平均通过0.7门课, 波动范围±1.0门
#通过率90%: 平均通过0.3门课, 波动范围±0.6门
基于负二项分布的营销预算决策
target_customers = 20
estimated_conversion_rate = 0.06
cost_per_contact = 50 # 每次营销接触的成本
mean, var = nbinstat(target_customers, estimated_conversion_rate)
total_cost = mean * cost_per_contact
cost_uncertainty = np.sqrt(var) * cost_per_contact
print(f"\n营销预算分析:")
print(f"目标: 获得{target_customers}个客户")
print(f"预计转化率: {estimated_conversion_rate:.1%}")
print(f"平均营销预算: ¥{total_cost:,.0f}")
print(f"预算波动范围: ±¥{cost_uncertainty:,.0f} (95%置信区间)")
#营销预算分析:
#目标: 获得20个客户
#预计转化率: 6.0%
#平均营销预算: ¥15,667
#预算波动范围: ±¥3,613 (95%置信区间)
风险评估
if cost_uncertainty / total_cost > 0.3:
print("⚠️ 高风险: 预算不确定性较高")
else:
print("✅ 低风险: 预算相对确定")
#✅ 低风险: 预算相对确定
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import NegativeBinomial, E, variance
def negative_binomial_mean_variance(input_str):
"""
对标Matlab的nbinstat函数，计算负二项分布的均值和方差。
参数:
input_str: 输入表达式字符串，格式为元组 (r, p)，其中r和p可以是标量、矩阵或符号。
返回:
均值矩阵和方差矩阵。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
M, V = None, None
import numpy as np
def nbinstat_simple(r, p):
"""
计算负二项分布的均值和方差（标量输入版本）
负二项分布定义：在获得r次成功前的失败次数
参数:
r (float): 需要的成功次数（正数）
p (float): 每次试验成功的概率 (0 ≤ p ≤ 1)
返回:
tuple: (mean, variance) 两个浮点数
"""
# 参数验证
if r <= 0:
raise ValueError(f"r 必须为正数，当前值: {r}")
if not (0 <= p <= 1):
raise ValueError(f"p 必须在 [0,1] 范围内，当前值: {p}")
# 特殊情况处理
if p == 0:
return np.inf, np.inf # 不可能成功，无限次失败
if p == 1:
return 0.0, 0.0 # 立即成功，无失败
# 计算公式
q = 1 - p
mean = r * q / p
variance = r * q / (p ** 2)
return mean, variance
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
r = int(expr[0])
p = float(expr[1])
M, V = nbinstat_simple(r, p)
else:
error = True
else:
error = True
return (M, V) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 标量输入
print(negative_binomial_mean_variance("(2, 0.3)"))
# (4.666666666666667, 15.555555555555555)
二项式系数或所有组合
b = nchoosek(n,k) 返回二项式系数, 这是从 n 项中一次取 k 项的组合的数目。n 和 k 必须为非负整数
C = nchoosek(v,k) 返回一个矩阵，其中包含了从向量 v 的元素中一次取 k 个元素的所有组合。矩阵 C 有 k 列和 m!/((m–k)! k!) 行，其中 m 为 length(v)
n — 可选项的数目,非负整数标量
k — 选中项的数目,非负整数标量
b — 二项式系数,非负标量值
C — v 中的所有组合,矩阵
示例1. 彩票与概率计算
双色球：从33个红球中选6个的组合数
print("双色球红球组合数:", nchoosek(33, 6))
#双色球红球组合数: 1107568.0
大乐透：前区35选5，后区12选2
front_comb = nchoosek(35, 5)
back_comb = nchoosek(12, 2)
total_comb = front_comb * back_comb
print(f"大乐透总组合数: {total_comb:,}")
#大乐透总组合数: 21,425,712.0
示例2. 团队建设与分组
从10人中选出3人组成项目团队的所有可能组合
team_combinations = nchoosek([1,2,3,4,5,6,7,8,9,10], 3)
print("3人团队组合数:", len(team_combinations))
print("前5个团队组合示例:")
for i in range(5):
print(f" 团队{i + 1}: {team_combinations[i]}")
#3人团队组合数: 360
#前5个团队组合示例:
# 团队1: [1.0, 2.0, 3.0]
# 团队2: [1.0, 2.0, 4.0]
# 团队3: [1.0, 2.0, 5.0]
# 团队4: [1.0, 2.0, 6.0]
# 团队5: [1.0, 2.0, 7.0]
部门分组：20人分成4个5人小组
group_combinations = nchoosek(20, 5)
print(f"\n5人小组的组合数: {group_combinations:,}")
#5人小组的组合数: 15,504.0
示例3. 产品组合与套餐设计
餐厅菜单：从8道主菜中选3道组成套餐
menu_items = ["牛排", "鱼排", "鸡排", "猪排", "羊排", "虾", "蟹", "素菜"]
menu_combinations = nchoosek(len(menu_items), 3)
print(f"3道菜套餐组合数: {menu_combinations}")
#3道菜套餐组合数: 56.0
实际生成套餐组合
menu_indices = nchoosek([1,2,3,4,5,6,7,8], 3)
部分套餐组合:
for i, combo in enumerate(menu_indices[:5]):
dishes = [menu_items[int(idx) - 1] for idx in combo]
print(f" 套餐{i + 1}: {', '.join(dishes)}")
# 套餐1: 牛排, 鱼排, 鸡排
# 套餐2: 牛排, 鱼排, 猪排
# 套餐3: 牛排, 鱼排, 羊排
# 套餐4: 牛排, 鱼排, 虾
# 套餐5: 牛排, 鱼排, 蟹
示例4. 网络与路由选择
网络节点间的路径选择
nodes = ["A", "B", "C", "D", "E", "F"]
print("网络节点间的双向连接数:", nchoosek(f"({len(nodes)}, 2)"))
#网络节点间的双向连接数: 15.0
生成所有可能的网络连接
connections = nchoosek([1,2,3,4,5,6], 2)
所有网络连接:
for conn in connections:
node1, node2 = nodes[int(conn[0]) - 1], nodes[int(conn[1]) - 1]
print(f" {node1} - {node2}")
#A - B
#A - C
#A - D
#A - E
#A - F
#B - C
#B - D
#B - E
#B - F
#C - D
#C - E
#C - F
#D - E
#D - F
#E - F
示例5. 质量控制与抽样检测
从100个产品中随机抽取5个进行质量检查
sample_combinations = nchoosek(100, 5)
print(f"质量抽检的样本组合数: {sample_combinations:,}")
#质量抽检的样本组合数: 75,287,520.0
如果知道有3个次品，检测到至少1个次品的概率
defective_prob = 1 - nchoosek(97, 5) / sample_combinations
print(f"抽到至少1个次品的概率: {defective_prob:.2%}")
#抽到至少1个次品的概率: 14.40%
示例6. 投资组合管理
从20只股票中构建投资组合
stocks = ["AAPL", "GOOGL", "MSFT", "AMZN", "TSLA", "META", "NVDA", "NFLX"]
for k in [2, 3, 4]:
comb_count = nchoosek(len(stocks), k)
print(f"{k}只股票的投资组合数: {comb_count}")
#2只股票的投资组合数: 28.0
#3只股票的投资组合数: 56.0
#4只股票的投资组合数: 70.0
生成具体的投资组合
portfolios = nchoosek([1,2,3,4,5,6,7,8], 3)
部分3股投资组合:
for i, port in enumerate(portfolios[:5]):
stock_combo = [stocks[int(idx) - 1] for idx in port]
print(f" 组合{i + 1}: {', '.join(stock_combo)}")
# 组合1: AAPL, GOOGL, MSFT
# 组合2: AAPL, GOOGL, AMZN
# 组合3: AAPL, GOOGL, TSLA
# 组合4: AAPL, GOOGL, META
# 组合5: AAPL, GOOGL, NVDA
示例7. 科学研究与实验设计
药物试验：从8种候选药物中选择3种进行组合试验
drug_combinations = nchoosek(8, 3)
print(f"3药组合试验数: {drug_combinations}")
#3药组合试验数: 56.0
基因研究：检测5个基因位点的组合效应
gene_combinations = nchoosek(5, 2) + nchoosek(5, 3) + nchoosek(5, 4)
print(f"多基因相互作用研究组合数: {gene_combinations}")
#3药组合试验数: 56.0000000000000
#多基因相互作用研究组合数: 25.0
示例8. 旅游路线规划
从6个城市中选择4个设计旅游路线
cities = ["北京", "上海", "广州", "深圳", "杭州", "成都"]
route_combinations = nchoosek(len(cities), 4)
print(f"4城市旅游路线组合数: {route_combinations}")
#4城市旅游路线组合数: 15.0
生成具体路线组合
routes = nchoosek([1,2,3,4,5,6], 4)
部分旅游路线:
for i, route in enumerate(routes[:5]):
city_combo = [cities[int(idx) - 1] for idx in route]
print(f" 路线{i + 1}: {' → '.join(city_combo)}")
# 路线1: 北京 → 上海 → 广州 → 深圳
# 路线2: 北京 → 上海 → 广州 → 杭州
# 路线3: 北京 → 上海 → 广州 → 成都
# 路线4: 北京 → 上海 → 深圳 → 杭州
# 路线5: 北京 → 上海 → 深圳 → 成都
示例9. 委员会与组织结构
公司从15个部门经理中选出5人组成执行委员会
committee_combinations = nchoosek(15, 5)
print(f"执行委员会组合数: {committee_combinations:,}")
#执行委员会组合数: 3,003.0
不同规模的委员会选择
committee_sizes = [3, 5, 7]
for size in committee_sizes:
comb = nchoosek(15, size)
print(f"{size}人委员会组合数: {comb:,}")
#3人委员会组合数: 455.0
#5人委员会组合数: 3,003.0
#7人委员会组合数: 6,435.0
示例10. 符号计算的数学应用
符号计算：二项式定理展开系数
for k in range(6):
coeff = nchoosek(n, k)
print(f" C(n, {k}) = {coeff}")
#C(n, 0) = 1.00000000000000
#C(n, 1) = n
#C(n, 2) = n*(n - 1)/2
#C(n, 3) = n*(n - 2)*(n - 1)/6
#C(n, 4) = n*(n - 3)*(n - 2)*(n - 1)/24
#C(n, 5) = n*(n - 4)*(n - 3)*(n - 2)*(n - 1)/120
帕斯卡三角形性质验证
left = nchoosek(n, k)
right = nchoosek(n-1, k) + nchoosek(n-1, k-1)
print(f"C(n, k) = C(n-1, k) + C(n-1, k-1)")
print(f"{left} = {right}")
#C(n, k) = C(n-1, k) + C(n-1, k-1)
#gamma(n + 1)/(gamma(k + 1)*gamma(-k + n + 1)) = gamma(n)/(gamma(-k + n)*gamma(k + 1)) + gamma(n)/(gamma(k)*gamma(-k + n + 1))
示例11. 项目团队组合分析
available_engineers = 12
needed_engineers = 4
available_designers = 8
needed_designers = 2
计算可能的团队组合数
engineer_combs = nchoosek(available_engineers, needed_engineers)
designer_combs = nchoosek(available_designers, needed_designers)
total_combs = engineer_combs * designer_combs
项目团队组合分析:
print(f"工程师组合数: {engineer_combs:,}")
print(f"设计师组合数: {designer_combs:,}")
print(f"总团队组合数: {total_combs:,}")
#工程师组合数: 495.0
#设计师组合数: 28.0
#总团队组合数: 13,860.0
资源利用效率
utilization = (needed_engineers + needed_designers) / (available_engineers + available_designers)
print(f"人力资源利用率: {utilization:.1%}")
#人力资源利用率: 30.0%
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import itertools
def nk_binomial_coefficient(input_str):
"""
对标Matlab的nchoosek函数，实现二项式系数计算和组合生成。
参数:
input_str: 输入表达式字符串，可以是元组 (n, k) 或 (vector, k)
返回:
二项式系数值 或 组合矩阵。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
v, n = expr
if isinstance(v, list) and n.is_number:
np_v = np.array(v, dtype=float).ravel()
combinations = list(itertools.combinations(np_v, n))
# 将组合转换为矩阵
result = sp.Matrix(combinations)
else:
b_result = sp.binomial(v, n).evalf()
# 简化表达式
result = sp.combsimp(b_result)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 计算二项式系数
print(nk_binomial_coefficient("(5, 2)"))
# 10.0000000000000
print(nk_binomial_coefficient("(sympify('n'), 2)"))
# n*(n - 1)/2
# 示例2: 生成组合矩阵
print(nk_binomial_coefficient("([1, 2, 3], 2)"))
# Matrix([[1.00000000000000, 2.00000000000000],
# [1.00000000000000, 3.00000000000000],
# [2.00000000000000, 3.00000000000000]])
负对数似然
nll=negloglik（X）返回用于拟合概率分布pd的数据的负对数似然函数的值,负对数似然值越小，表示数据与分布的拟合越好。
mu — 向量X的中心值, 默认值是1.5, 用于计算每个数据点的概率密度
sigma — 向量X的标准差, 代表分布的离散程度, 默认值是0.5, 用于计算每个数据点的概率密度
normal — 正态分布, 数据对称分布在均值附近。
weibull — 威布尔分布, 描述故障，寿命等非对称数据。
X — 输入值, 向量
示例1. 产品质量寿命分析
产品寿命数据（小时）
lifetimes = [1200, 1500, 1800, 2100, 2400, 2700, 3000, 3300, 3600, 3900]
产品寿命数据负对数似然分析:
正态分布假设下的负对数似然
normal_nll = negloglik(lifetimes, normal)
print(f"正态分布负对数似然: {normal_nll}")
#正态分布负对数似然: 81.7782760803417
Weibull分布假设下的负对数似然
weibull_nll = negloglik(flifetimes, weibull)
print(f"Weibull分布负对数似然: {weibull_nll}")
#Weibull分布负对数似然: 81.5903993676933
比较模型拟合优度
if weibull_nll < normal_nll:
print("Weibull分布拟合更好")
else:
print("正态分布拟合更好")
#Weibull分布拟合更好
示例2. 金融收益率建模
股票日收益率数据（%）
returns = [1.2, -0.5, 2.3, -1.8, 0.7, 1.5, -0.9, 1.1, -2.1, 0.8,
1.6, -0.3, 2.1, -1.2, 0.9, 1.8, -0.7, 1.3, -1.5, 0.6]
股票收益率分布拟合:
returns_nll_normal = negloglik(returns, normal)
returns_nll_weibull = negloglik(returns[:10]}, weibull) # Weibull需要正数
print(f"收益率正态分布负对数似然: {returns_nll_normal}")
print(f"收益率绝对值Weibull分布负对数似然: {returns_nll_weibull}")
#收益率正态分布负对数似然: 33.9719280933577
#收益率绝对值Weibull分布负对数似然: 6.9010604226196 + 1.1886888165027*I
示例3. 设备故障时间预测
机械设备故障间隔时间（天）
failure_intervals = [45, 78, 32, 91, 56, 67, 83, 29, 71, 62, 88, 51]
设备故障间隔时间分析:
failure_nll_normal = negloglik(failure_intervals, normal)
failure_nll_weibull = negloglik(failure_intervals, weibull)
print(f"正态分布负对数似然: {failure_nll_normal}")
print(f"Weibull分布负对数似然: {failure_nll_weibull}")
#正态分布负对数似然: 52.9110109482821
#Weibull分布负对数似然: 52.6748040566129
模型选择
models = {
"正态分布": failure_nll_normal,
"Weibull分布": failure_nll_weibull
}
best_model = min(models, key=models.get)
print(f"最佳拟合分布: {best_model} (负对数似然最小)")
#最佳拟合分布: Weibull分布 (负对数似然最小)
示例4. 气候数据分析
月降水量数据（毫米）
precipitation = [85, 92, 78, 105, 88, 96, 112, 67, 89, 101, 95, 83]
降水量分布拟合:
precip_nll_normal = negloglik(precipitation, normal)
precip_nll_weibull = negloglik(precipitation, weibull)
print(f"正态分布负对数似然: {precip_nll_normal}")
print(f"Weibull分布负对数似然: {precip_nll_weibull}")
#正态分布负对数似然: 46.5212405194057
#Weibull分布负对数似然: 46.6173540303342
AIC准则比较（AIC = 2k + 2NLL，k为参数个数）
aic_normal = 2 * 2 + 2 * precip_nll_normal # 正态分布有2个参数
aic_weibull = 2 * 2 + 2 * precip_nll_weibull # Weibull分布有2个参数
print(f"正态分布AIC: {aic_normal}")
print(f"Weibull分布AIC: {aic_weibull}")
#正态分布AIC: 97.0424810388115
#Weibull分布AIC: 97.2347080606684
示例5. 医疗数据生存分析
患者生存时间（月）
survival_times = [12, 18, 24, 30, 36, 42, 48, 54, 60, 66, 72, 78, 84]
患者生存时间分布分析:
survival_nll_normal = negloglik(survival_times, normal)
survival_nll_weibull = negloglik(survival_times, weibull)
print(f"正态分布负对数似然: {survival_nll_normal}")
print(f"Weibull分布负对数似然: {survival_nll_weibull}")
#正态分布负对数似然: 58.8929466741246
#Weibull分布负对数似然: 58.5462173753973
if survival_nll_weibull < survival_nll_normal:
print("Weibull分布更适合生存时间建模")
#Weibull分布更适合生存时间建模
示例6. 生产制造过程控制
零件尺寸测量数据（毫米）
dimensions = [10.02, 10.01, 10.03, 9.99, 10.00, 10.02, 9.98, 10.01, 10.00, 10.02]
零件尺寸质量控制分析:
dims_nll_normal = negloglik(dimensions, normal)
print(f"尺寸数据正态分布负对数似然: {dims_nll_normal}")
#尺寸数据正态分布负对数似然: -28.0117754193544
计算过程能力指数所需参数
data_array = float(dimensions)
mu_est = mean(data_array)
sigma_est = std(data_array, ddof=1)
print(f"估计均值: {mu_est:.4f}")
print(f"估计标准差: {sigma_est:.4f}")
#估计均值: 10.0080
#估计标准差: 0.0155
示例7. 可靠性工程应用
电子元件寿命测试数据（千小时）
component_life = [1.5, 2.3, 3.1, 4.2, 5.0, 6.3, 7.1, 8.4, 9.2, 10.5]
电子元件可靠性分析:
comp_nll_weibull = negloglik(component_life, weibull)
print(f"Weibull分布负对数似然: {comp_nll_weibull}")
#Weibull分布负对数似然: 24.3940129491024
可靠性指标计算
A_est = 8.5 # 尺度参数估计
B_est = 2.1 # 形状参数估计
reliability_1000h = exp(-(1 / A_est) ** B_est)
print(f"1000小时可靠度估计: {reliability_1000h:.3f}")
#1000小时可靠度估计: 0.989
示例8. 符号计算与理论分析
使用符号数据进行理论分析
symbolic_data = [x1, x2, x3, x4]
symbolic_nll = negloglik(symbolic_data, normal)
符号正态分布负对数似然:
print(symbolic_nll)
#2*mu**2/sigma**2 + log(4*pi**2*sigma**4) - log(exp(x1*(mu - x1/2)/sigma**2)) - log(exp(x2*(mu - x2/2)/sigma**2)) - \
log(exp(x3*(mu - x3/2)/sigma**2)) - log(exp(x4*(mu - x4/2)/sigma**2))
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from sympy import symbols, log, exp, Piecewise
def matrix_to_vector(matrix):
"""
将矩阵转换为向量（列表形式）。若输入不是向量，返回None。
参数:
matrix: SymPy矩阵对象
返回:
列表形式的向量，若非向量则返回None。
"""
if not isinstance(matrix, sp.Matrix):
raise ValueError("输入必须是一个 sympy.Matrix 对象")
# 获取矩阵的形状
rows, cols = matrix.shape
# 如果只有一行或一列，则将其转换为列表
if rows == 1:
# 行向量：直接转换为列表
return list(matrix.row(0))
elif cols == 1:
# 列向量：提取每一行的第一个元素并转换为列表
return [matrix[i, 0] for i in range(rows)]
else:
# 如果不是向量，返回空列表
return []
def nll_dist(input_str, distribution='normal'):
"""
计算分布的负对数似然函数，支持正态分布和Weibull分布。
参数:
input_str: 输入数据（字符串形式的向量，如 "[1,2,3]"）
distribution: 分布类型，可选 'normal' 或 'weibull'
params: 已知参数（字典，例如 {'A': 26.5, 'B': 3.27}），若为None则从数据估计
返回:
数值结果（若输入为数值且参数已知）或符号表达式
"""
try:
# 解析输入数据
expr = sp.sympify(input_str)
matrix = sp.Matrix(expr) if isinstance(expr, list) else None
if matrix is None:
return "错误：输入无法转换为矩阵"
vector = matrix_to_vector(matrix)
if vector is None:
return "错误：输入不是向量"
# 定义符号变量和参数
if distribution == 'normal':
mu, sigma = symbols('mu sigma', real=True, positive=True)
z = symbols('z', real=True)
# 正态分布PDF
pdf = (1 / (sigma * sp.sqrt(2 * sp.pi))) * exp(-(z - mu) ** 2 / (2 * sigma ** 2))
param_symbols = {'mu': mu, 'sigma': sigma}
elif distribution == 'weibull':
A, B = symbols('A B', real=True, positive=True) # A:尺度参数, B:形状参数
z = symbols('z', real=True, nonnegative=True)
# Weibull分布PDF
pdf = (B / A) * (z / A) ** (B - 1) * exp(-(z / A) ** B)
param_symbols = {'A': A, 'B': B}
else:
return "错误：不支持的分布类型"
# 负对数似然函数
neg_log_likelihood = -log(pdf)
# 总负对数似然
total_nll = sum([neg_log_likelihood.subs(z, val) for val in vector])
# 处理参数
# 从数据估计参数
if all(val.is_number for val in vector):
data = np.array([float(v) for v in vector], dtype=float)
if distribution == 'normal':
mu_est = np.mean(data)
sigma_est = np.std(data, ddof=0)
result = total_nll.subs({mu: mu_est, sigma: sigma_est}).evalf()
elif distribution == 'weibull':
# Weibull参数估计（需使用scipy等库进行更精确估计）
from scipy.stats import weibull_min
shape, loc, scale = weibull_min.fit(data, floc=0)
result = total_nll.subs({A: scale, B: shape}).evalf()
else:
result = total_nll # 返回符号表达式
return result
except Exception as e:
return f"错误：{str(e)}"
# 示例使用
if __name__ == "__main__":
# 示例1: 使用已知Weibull参数
data = "[2.5, 3.0, 5.1, 7.2, 9.5]"
params = {'A': 26.5079, 'B': 3.27193}
print(nll_dist(data, distribution='weibull'))
# 11.6157041432461
# 示例2：数值向量
input2 = "[1, 2, 3, 4]"
print(f"示例1输入: {input2}")
print("输出:", nll_dist(input2))
# 6.12204123544711
2的更高次幂的指数
P = nextpowTwo(A)返回对A中每个元素满足的最小的2的幂的指数. 按照惯例,nextpowTwo(0)返回零.
您可以使用nextpowTwo填充传递到fft的信号。当信号长度并非2次幂时，这样做可以加快FFT的运算速度
A — 输入值,由实数组成的标量,向量,数组
示例1. 信号处理与FFT优化
FFT计算时选择最合适的长度（2的幂次）
signal_lengths = [1000, 2047, 4096, 8191, 10000]
FFT长度优化:
for length in signal_lengths:
next_pow = nextpowTwo(length)
optimal_length = 2 ** next_pow
print(f"信号长度{length:5d} → 2^{next_pow} = {optimal_length:5d} (最优FFT长度)")
#信号长度 1000 → 2^10 = 1024 (最优FFT长度)
#信号长度 2047 → 2^11 = 2048 (最优FFT长度)
#信号长度 4096 → 2^12 = 4096 (最优FFT长度)
#信号长度 8191 → 2^13 = 8192 (最优FFT长度)
#信号长度10000 → 2^14 = 16384 (最优FFT长度)
实际信号处理应用
sample_rates = [44100, 48000, 96000, 192000]
采样率对应的FFT长度:
for rate in sample_rates:
# 假设我们想要约1秒的FFT窗口
window_samples = rate
next_pow = nextpowTwo(window_samples)
fft_length = 2 ** next_pow
print(f"采样率{rate:6d}Hz → FFT长度: {fft_length}")
#采样率 44100Hz → FFT长度: 65536
#采样率 48000Hz → FFT长度: 65536
#采样率 96000Hz → FFT长度: 131072
#采样率192000Hz → FFT长度: 262144
示例2. 图像处理与纹理映射
图像尺寸调整为2的幂次（OpenGL纹理要求）
image_sizes = [
[640, 480], # VGA
[1280, 720], # HD
[1920, 1080], # Full HD
[3840, 2160] # 4K
]
图像纹理尺寸优化:
for width, height in image_sizes:
next_pow_w = nextpowTwo(width)
next_pow_h = nextpowTwo(height)
tex_width = 2 ** next_pow_w
tex_height = 2 ** next_pow_h
print(f"原始: {width:4d}×{height:4d} → 纹理: {tex_width:4d}×{tex_height:4d}")
#原始: 640× 480 → 纹理: 1024× 512
#原始: 1280× 720 → 纹理: 2048×1024
#原始: 1920×1080 → 纹理: 2048×2048
#原始: 3840×2160 → 纹理: 4096×4096
示例3. 内存分配与数据结构
动态数组扩容策略
current_sizes = [100, 500, 1000, 5000, 10000]
动态数组扩容策略:
for size in current_sizes:
next_pow = nextpowTwo(size)
new_size = 2 ** next_pow
growth_factor = new_size / size
print(f"当前大小{size:5d} → 扩容到{new_size:5d} (增长{growth_factor:.1f}倍)")
#当前大小 100 → 扩容到 128 (增长1.3倍)
#当前大小 500 → 扩容到 512 (增长1.0倍)
#当前大小 1000 → 扩容到 1024 (增长1.0倍)
#当前大小 5000 → 扩容到 8192 (增长1.6倍)
#当前大小10000 → 扩容到16384 (增长1.6倍)
哈希表桶大小优化
hash_table_sizes = [127, 255, 511, 1023, 2047]
哈希表桶大小优化:
for size in hash_table_sizes:
next_pow = nextpowTwo(size)
optimal_size = 2 ** next_pow
print(f"建议大小{size:4d} → 优化为{optimal_size:4d} (2^{next_pow})")
#建议大小 127 → 优化为 128 (2^7)
#建议大小 255 → 优化为 256 (2^8)
#建议大小 511 → 优化为 512 (2^9)
#建议大小1023 → 优化为1024 (2^10)
#建议大小2047 → 优化为2048 (2^11)
示例4. 音频处理与缓冲区设计
音频缓冲区大小优化
buffer_sizes_ms = [10, 20, 50, 100, 200] # 毫秒
sample_rate = 44100 # Hz
音频缓冲区大小计算:
for buffer_ms in buffer_sizes_ms:
samples_needed = int(sample_rate * buffer_ms / 1000)
next_pow = nextpowTwo(samples_needed)
optimal_buffer = 2 ** next_pow
actual_ms = optimal_buffer / sample_rate * 1000
print(f"{buffer_ms:3d}ms需求: {samples_needed:5d}样本 → 实际: {optimal_buffer:5d}样本 ({actual_ms:.1f}ms)")
# 10ms需求: 441样本 → 实际: 512样本 (11.6ms)
# 20ms需求: 882样本 → 实际: 1024样本 (23.2ms)
# 50ms需求: 2205样本 → 实际: 4096样本 (92.9ms)
#100ms需求: 4410样本 → 实际: 8192样本 (185.8ms)
#200ms需求: 8820样本 → 实际: 16384样本 (371.5ms)
示例5. 网络数据包与帧大小
网络协议帧大小优化
packet_sizes = [512, 1024, 1500, 4096, 8192] # 字节
网络数据包大小优化:
for size in packet_sizes:
next_pow = nextpowTwo(size)
optimal_size = 2 ** next_pow
efficiency = size / optimal_size * 100
print(f"原始{size:4d}字节 → 对齐{optimal_size:4d}字节 (效率{efficiency:.1f}%)")
#原始 512字节 → 对齐 512字节 (效率100.0%)
#原始1024字节 → 对齐1024字节 (效率100.0%)
#原始1500字节 → 对齐2048字节 (效率73.2%)
#原始4096字节 → 对齐4096字节 (效率100.0%)
#原始8192字节 → 对齐8192字节 (效率100.0%)
视频帧压缩块大小
frame_blocks = [8, 16, 32, 64]
视频编码块大小:
for block in frame_blocks:
next_pow = nextpowTwo(block)
print(f"块大小{block:2d} → 2^{next_pow} = {2 ** next_pow}")
#块大小 8 → 2^3 = 8
#块大小16 → 2^4 = 16
#块大小32 → 2^5 = 32
#块大小64 → 2^6 = 64
示例6. 游戏开发与资源管理
游戏资源内存对齐
resource_sizes = [
("角色模型", 15360),
("纹理贴图", 524288),
("音频文件", 1048576),
("关卡数据", 4194304)
]
游戏资源内存对齐:
for name, size in resource_sizes:
next_pow = nextpowTwo(size)
aligned_size = 2 ** next_pow
memory_kb = aligned_size / 1024
print(f"{name:8}: {size:8,d}字节 → {aligned_size:8,d}字节 ({memory_kb:,.0f}KB)")
#角色模型 : 15,360字节 → 16,384字节 (16KB)
#纹理贴图 : 524,288字节 → 524,288字节 (512KB)
#音频文件 : 1,048,576字节 → 1,048,576字节 (1,024KB)
#关卡数据 : 4,194,304字节 → 4,194,304字节 (4,096KB)
示例7. 数据库与索引优化
数据库页面大小优化
page_sizes_kb = [2, 4, 8, 16, 32, 64]
数据库页面大小优化:
for size_kb in page_sizes_kb:
size_bytes = size_kb * 1024
next_pow = nextpowTwo(size_bytes)
optimal_bytes = 2 ** next_pow
optimal_kb = optimal_bytes / 1024
print(f"建议{size_kb:2d}KB → 实际{optimal_kb:2.0f}KB (2^{next_pow})")
#建议 2KB → 实际 2KB (2^11)
#建议 4KB → 实际 4KB (2^12)
#建议 8KB → 实际 8KB (2^13)
#建议16KB → 实际16KB (2^14)
#建议32KB → 实际32KB (2^15)
#建议64KB → 实际64KB (2^16)
示例8. 科学计算与数值分析
数值积分网格大小
grid_sizes = [100, 500, 1000, 2000, 5000]
数值积分网格优化:
for size in grid_sizes:
next_pow = nextpowTwo(size)
optimal_grid = 2 ** next_pow
improvement = (optimal_grid - size) / size * 100
print(f"网格{size:4d} → 优化{optimal_grid:4d} (增加{improvement:+.1f}%)")
#网格 100 → 优化 128 (增加+28.0%)
#网格 500 → 优化 512 (增加+2.4%)
#网格1000 → 优化1024 (增加+2.4%)
#网格2000 → 优化2048 (增加+2.4%)
#网格5000 → 优化8192 (增加+63.8%)
矩阵运算分块大小
matrix_dims = [256, 512, 1024, 2048]
矩阵分块计算优化:
for dim in matrix_dims:
next_pow = nextpowTwo(dim)
block_size = 2 ** next_pow
print(f"矩阵维度{dim:4d} → 分块大小{block_size:4d}")
#矩阵维度 256 → 分块大小 256
#矩阵维度 512 → 分块大小 512
#矩阵维度1024 → 分块大小1024
#矩阵维度2048 → 分块大小2048
示例9. 嵌入式系统与硬件优化
微控制器内存分配
memory_chunks = [256, 512, 1024, 2048, 4096] # 字节
嵌入式系统内存分配:
for chunk in memory_chunks:
next_pow = nextpowTwo(chunk)
aligned_chunk = 2 ** next_pow
print(f"需求{chunk:4d}字节 → 分配{aligned_chunk:4d}字节")
#需求 256字节 → 分配 256字节
#需求 512字节 → 分配 512字节
#需求1024字节 → 分配1024字节
#需求2048字节 → 分配2048字节
#需求4096字节 → 分配4096字节
DMA传输缓冲区
dma_buffers = [128, 256, 512, 1024]
DMA缓冲区优化:
for buffer in dma_buffers:
next_pow = nextpowTwo(buffer)
optimal_buffer = 2 ** next_pow
print(f"缓冲区{buffer:4d} → 优化{optimal_buffer:4d}")
#缓冲区 128 → 优化 128
#缓冲区 256 → 优化 256
#缓冲区 512 → 优化 512
#缓冲区1024 → 优化1024
示例10. 符号计算与理论分析
符号表达式分析
symbolic_inputs = [n, n+1, 2n, n^2]
符号输入的2的幂次分析:
for expr in symbolic_inputs:
result = nextpowTwo(expr)
print(f"nextpow2({expr:5}) = {result}")
#nextpow2(n ) = ceiling(log(Abs(n))/log(2))
#nextpow2(n+1 ) = ceiling(log(Abs(n + 1))/log(2))
#nextpow2(2n ) = ceiling(log(2*Abs(n))/log(2))
#nextpow2(n^2 ) = ceiling(log(Abs(n**2))/log(2))
复杂表达式
complex_exprs = [
x + y,
2x + 1,
x^2 + y^2",
sqrt(x)
]
复杂符号表达式:
for expr in complex_exprs:
try:
result = nextpowTwo(expr)
print(f"nextpow2({expr:10}) = {result}")
except:
print(f"nextpow2({expr:10}) = 无法计算")
#nextpow2(x + y ) = ceiling(log(Abs(x + y))/log(2))
#nextpow2(2*x + 1 ) = ceiling(log(Abs(2*x + 1))/log(2))
#nextpow2(x^2 + y^2) = ceiling(log(Abs(x**2 + y**2))/log(2))
#nextpow2(sqrt(x) ) = ceiling(log(Abs(sqrt(x)))/log(2))
示例11. 系统参数优化分析
各种系统参数
parameters = {
"音频缓冲区": 4410, # 100ms @ 44.1kHz
"视频帧缓冲区": 1920 * 1080 * 3, # 1080p RGB
"网络发送窗口": 65535, # TCP最大窗口
"数据库缓存": 1048576, # 1MB
"图像处理块": 512 # 处理块大小
}
for name, value in parameters.items():
next_pow = nextpowTwo(value)
optimized = 2 ** next_pow
ratio = optimized / value
print(f"{name:12}: {value:8,d} → {optimized:8,d} (x{ratio:.2f})")
#音频缓冲区 : 4,410 → 8,192 (x1.86)
#视频帧缓冲区 : 6,220,800 → 8,388,608 (x1.35)
#网络发送窗口 : 65,535 → 65,536 (x1.00)
#数据库缓存 : 1,048,576 → 1,048,576 (x1.00)
#图像处理块 : 512 → 512 (x1.00)
示例12. 性能测试与基准比较
sizes = [100, 1000, 10000, 100000]
不同数据规模的2的幂次对齐:
for size in sizes:
next_pow = nextpowTwo(size)
aligned_size = 2 ** next_pow
overhead = (aligned_size - size) / size * 100
efficiency = "优秀" if overhead < 10 else "良好" if overhead < 25 else "一般"
print(f"大小{size:6,d} → 对齐{aligned_size:6,d} (开销{overhead:5.1f}%) - {efficiency}")
#大小 100 → 对齐 128 (开销 28.0%) - 一般
#大小 1,000 → 对齐 1,024 (开销 2.4%) - 优秀
#大小10,000 → 对齐16,384 (开销 63.8%) - 一般
#大小100,000 → 对齐131,072 (开销 31.1%) - 一般
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def next_pow_2_e(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
def evaluation_new_pow2(x_val):
"""核心计算函数"""
if x_val == 0:
return 0
else:
return sp.ceiling(sp.log(abs(x_val), 2))
# 处理单个数值
if expr.is_number or expr.free_symbols:
result = evaluation_new_pow2(expr)
# 无法处理的类型（如符号表达式）
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: next_pow_2_e('5') =", next_pow_2_e('5'))
# 3
# 示例2: 零值
print("示例2: next_pow_2_e('0') =", next_pow_2_e('0'))
# 0
# 示例3: 符号输入
print("示例3: next_pow_2_e('x') =", next_pow_2_e('x'))
# ceiling(log(Abs(x))/log(2))
尼柯尔斯图
nicholsplot(tf)
尼柯尔斯图是控制系统频域分析的重要工具，以 分贝（dB）为纵轴（增益）、度（°）为横轴（相位），绘制开环传递函数的频率响应曲线。
电机位置控制系统:
相位裕度 ≈ 45°（曲线离 -180° 线的垂直距离）。
幅值裕度 ≈ 10 dB（曲线离 0 dB 的水平距离）。
谐振峰值 ≈ 2 dB → 瞬态响应轻微振荡。
nicholsplot(10,s*(0.5s+1)*(0.1s+1))
无人机姿态控制器:
相位裕度 > 60° → 高稳定性。
带宽 ≈ 3 rad/s（0 dB 交点频率）→ 快速响应。
无谐振峰值（曲线避开 +3 dB 区域）→ 平滑跟踪。
nicholsplot(2*(s+5),s*(s+1)*(s+10))
机械臂力控制系统:
相位裕度 < 30° → 需增加补偿器。
幅值裕度 ≈ 6 dB → 临界稳定。
谐振峰值 +5 dB → 操作时存在抖动。
nicholsplot(0.8*(s+3),s^2*(s+4))
温度调节系统:
相位裕度 ≈ 70° → 强鲁棒性。
带宽 ≈ 0.5 rad/s → 缓慢调节，避免超调。
无谐振 → 温度变化平稳。
nicholsplot(1.2,s*(2s+1)*(0.2s+1))
一阶惯性系统
分析温度传感器或 RC 电路的动态响应
检查相位裕度（曲线在 -90° 附近是否远离临界点）
nicholsplot([1], [1, 2])
欠阻尼二阶系统
弹簧-质量阻尼系统或 LCR 振荡电路
观察谐振峰（曲线在 -120°~-150° 处的凸起是否接触 +3dB 等幅线）
nicholsplot([1], [1, 0.6, 1])
含微分项的系统
带速度反馈的电机位置控制
验证零点（s=-0.5s=−0.5）是否提升相位裕度（曲线向左移动）
nicholsplot([1, 0.5], [1, 2, 4])
非最小相位系统
飞机俯仰动力学或化学过程控制
右半平面零点（s=1s=1）导致高频段相位滞后（曲线向右卷绕）
nicholsplot([-1,1],[1,4,3])
三阶系统
伺服电机或机械臂关节控制
评估稳定性：曲线在 -180° 时幅值需 <0dB（避免包围临界点）
nicholsplot([2],[1,6,5,0])
叠加尼柯尔斯图
nicholsplot([tf1],[tf2]) 是在同一坐标系中绘制多个开环传递函数的尼柯尔斯曲线，用于直观比较不同系统的稳定性、鲁棒性和动态性能。
汽车悬架系统优化对比: 比较被动悬架 vs 主动悬架
主动悬架曲线（红色）相位滞后减少 25° → 过弯稳定性提升
谐振峰值从 +4dB 降至 +1.2dB → 颠簸路面振动减弱
带宽从 2 rad/s 增至 8 rad/s → 响应速度提高 3 倍
nicholsplot([0.8,(s^2+1.2s+4)*(0.02s+1)], # 被动悬架 (传统弹簧-阻尼)
[3.5*(s+5), s*(s^2+2s+6)*(0.01*s+1)]) # 主动悬架 (带液压作动器)
航空发动机燃油控制系统: 不同飞行高度下的鲁棒性验证
需增加增益调度补偿器补偿高度变化影响
nicholsplot([50,s*(s+1)*(0.1*s^2+0.5s+2)], # 海平面模式 (h=0)
[28,s*(s+0.7)*(0.15*s^2+0.3s+1.8)]) # 巡航高度模式 (h=10km)
医疗输液泵安全验证: 正常状态 vs 管路堵塞故障
故障曲线（橙色）逼近 -180°线 → 相位裕度从 65° 降至 32°
增益交点左移 → 带宽缩减 60%（触发堵塞报警阈值）
曲线进入 +6dB 危险区 → 需启动压力保护机制
nicholsplot([2.4,s*(0.5s+1)*(0.1s+1)], # 正常流动 (标称模型)
[1.2,s*(1.2s+1)*(0.15s+1)]) # 管路堵塞 (流量下降40%)
机器人关节控制器迭代: PID 控制器 vs 自适应模糊控制器
相位穿越频率从 5 rad/s → 12 rad/s → 动作更敏捷
避开 +3dB 谐振区 → 关节运动超调量从 15% 降至 3%
曲线形态更平坦 → 负载扰动敏感度降低
nicholsplot([3*(s+0.5),s^2*(0.02s+1)], # 传统PID
[4.5*(s^2+2s+1), s^2*(0.01*s^2+0.1s+1)]) # 模糊自适应控制器
光伏逆变器并网性能: 晴天 vs 阴天发电工况
阴天曲线下移 6dB → 谐波失真率增加 2.3%
相位裕度保持 >50° → 满足并网稳定标准
0dB 交点频率变化 0.5 rad/s → 需调整锁相环参数
nicholsplot([8.5,(0.05*s^2+0.3s+1)*(s+10)], # 最大功率输出 (晴天)
[3.2,(0.08*s^2+0.4s+1.2)*(s+8)]) # 降额运行 (阴天)
向量范数和矩阵范数
n = norm(v) 返回向量 v 的欧几里德范数. 此范数也称为 2-范数,向量模或欧几里德长度.
n = norm(v,p) 返回广义向量 p 范数.
n = norm(X) 返回矩阵 X 的 2-范数或最大奇异值，该值近似于 max(svd(X)).
n = norm(X,p) 返回矩阵 X 的 p-范数，其中 p 为 1,2 或 Inf
如果 p = 1, 则 n 是矩阵的最大绝对列之和.
如果 p = 2, 则 n 近似于 max(svd(X)). 此值等效于 norm(X)
如果 p = Inf, 则 n 是矩阵的最大绝对行之和.
n = norm(X,"fro") 返回矩阵或数组X的Frobenius范数
v — 输入向量, 向量
X — 输入数组, 矩阵, 数组
示例1. 机器学习与特征工程
特征向量归一化
feature_vectors = [
[1.2, -0.8, 2.1, -1.5], # 样本1
[0.9, 1.3, -0.7, 2.2], # 样本2
[-1.1, 0.5, 1.8, -0.9] # 样本3
]
特征向量范数分析:
for i, features in enumerate(feature_vectors):
l2_norm = norm(features)
l1_norm = norm(features, 1)
inf_norm = norm(features, inf)
print(f"样本{i + 1}: L2={l2_norm:.3f}, L1={l1_norm:.3f}, L∞={inf_norm:.3f}")
#样本1: L2=2.956, L1=5.600, L∞=2.100
#样本2: L2=2.798, L1=5.100, L∞=2.200
#样本3: L2=2.347, L1=4.300, L∞=1.800
特征归一化应用
for i, features in enumerate(feature_vectors):
l2_norm = float(norm(features))
normalized = [x / l2_norm for x in features]
print(f"样本{i + 1}: {[f'{x:.3f}' for x in normalized]}")
#样本1: ['0.406', '-0.271', '0.710', '-0.507']
#样本2: ['0.322', '0.465', '-0.250', '0.786']
#样本3: ['-0.469', '0.213', '0.767', '-0.383']
示例2. 图像处理与计算机视觉
图像梯度矩阵分析
gradient_matrix = [
[12, -8, 15, -6],
[-5, 20, -12, 8],
[18, -10, 25, -15],
[-7, 14, -9, 11]
]
图像梯度矩阵范数分析:
fro_norm = norm(gradient_matrix, fro)
l2_norm = norm(gradient_matrix)
l1_norm = norm(gradient_matrix, 1)
inf_norm = norm(gradient_matrix, @inf)
print(f"Frobenius范数: {fro_norm:.3f} (总梯度强度)")
print(f"谱范数 (L2): {l2_norm:.3f} (最大奇异值)")
print(f"列和范数 (L1): {l1_norm:.3f} (最大列和)")
print(f"行和范数 (L∞): {inf_norm:.3f} (最大行和)")
#Frobenius范数: 53.132 (总梯度强度)
#谱范数 (L2): 53.132 (最大奇异值)
#列和范数 (L1): 61.000 (最大列和)
#行和范数 (L∞): 68.000 (最大行和)
图像相似性比较
def image_similarity(img1, img2):
"""计算两幅图像的差异范数"""
diff = [[img1[i][j] - img2[i][j] for j in range(len(img1[0]))]
for i in range(len(img1))]
return norm(diff, fro)
img_a = [[100, 150], [200, 50]]
img_b = [[105, 145], [195, 55]]
similarity = image_similarity(img_a, img_b)
print(f"\n图像差异Frobenius范数: {similarity:.3f}")
#图像差异Frobenius范数: 10.000
示例3. 控制系统与稳定性分析
系统矩阵范数分析
system_matrix = [
[0.8, 0.1, -0.2],
[0.05, 0.9, 0.1],
[-0.1, 0.05, 0.85]
]
控制系统矩阵范数分析:
matrix_norms = {
"谱范数 (L2)": norm(system_matrix),
"列和范数 (L1)": norm(system_matrix, 1),
"行和范数 (L∞)": norm(system_matrix, inf),
"Frobenius范数": norm(system_matrix, fro)
}
for name, norm_val in matrix_norms.items():
stability = "稳定" if float(norm_val) < 1 else "可能不稳定"
print(f"{name}: {norm_val:.4f} - {stability}")
#谱范数 (L2): 1.4992 - 可能不稳定
#列和范数 (L1): 1.1500 - 可能不稳定
#行和范数 (L∞): 1.1000 - 可能不稳定
#Frobenius范数: 1.4992 - 可能不稳定
状态向量分析
state_vector = [1.5, -0.8, 2.1]
state_norm = norm(state_vector)
print(f"\n状态向量L2范数: {state_norm:.3f} (系统能量)")
#状态向量L2范数: 2.702 (系统能量)
示例4. 数值分析与误差估计
线性方程组误差分析
A = [[4, 1], [1, 3]] # 系数矩阵
b = [1, 2] # 右侧向量
x_approx = [0.2, 0.6] # 近似解
计算残差范数
residual = [b[0] - (A[0][0] * x_approx[0] + A[0][1] * x_approx[1]),
b[1] - (A[1][0] * x_approx[0] + A[1][1] * x_approx[1])]
residual_norm = norm(residual)
print(f"残差L2范数: {residual_norm:.6f}")
#残差L2范数: 0.400000
条件数估计（使用谱范数）
A_norm = norm(A)
A_inv_norm = norm([[0.2727, -0.0909], [-0.0909, 0.3636]]) # 近似逆
cond_number = float(A_norm) * float(A_inv_norm)
print(f"矩阵条件数估计: {cond_number:.4f}")
#矩阵条件数估计: 2.4543
示例5. 信号处理与音频分析
音频信号能量计算
audio_frames = [
[0.12, -0.08, 0.15, -0.06, 0.09],
[-0.05, 0.11, -0.07, 0.13, -0.04],
[0.08, -0.09, 0.14, -0.05, 0.10]
]
音频帧能量分析 (L2范数平方):
for i, frame in enumerate(audio_frames):
frame_energy = float(norm(frame)) ** 2
print(f"帧{i + 1}: 能量 = {frame_energy:.6f}")
#帧1: 能量 = 0.055000
#帧2: 能量 = 0.038000
#帧3: 能量 = 0.046600
信号相似性（相关性矩阵范数）
correlation_matrix = [
[1.0, 0.8, 0.3],
[0.8, 1.0, 0.5],
[0.3, 0.5, 1.0]
]
corr_norm = norm(correlation_matrix, 2)
print(f"\n相关矩阵谱范数: {corr_norm:.3f} (最大特征值)")
#相关矩阵谱范数: 2.095 (最大特征值)
示例6. 量子计算与状态向量
量子态向量范数（概率幅）
quantum_state = [0.6, 0.8] # |ψ⟩ = 0.6|0⟩ + 0.8|1⟩
state_norm = norm(quantum_state)
print(f"量子态向量L2范数: {state_norm:.6f}")
#量子态向量L2范数: 1.0
验证归一化条件
probability_sum = sum(x ** 2 for x in quantum_state)
print(f"概率和: {probability_sum:.6f}")
print(f"归一化检查: {'是' if abs(state_norm - 1.0) < 1e-10 else '否'}")
#概率和: 1.000000
#归一化检查: 是
量子门矩阵范数
hadamard_gate = [[1 / sqrt(2), 1 / sqrt(2)],
[1 / sqrt(2), -1 / sqrt(2)]]
gate_norm = norm(hadamard_gate)
print(f"Hadamard门谱范数: {gate_norm:.6f}")
#Hadamard门谱范数: 1.414214
示例7. 金融风险分析与投资组合
投资组合权重向量
portfolio_weights = [0.3, 0.4, 0.2, 0.1] # 四个资产的权重
weight_norm = norm(portfolio_weights)
print(f"投资组合权重L2范数: {weight_norm:.6f}")
#投资组合权重L2范数: 0.547723
协方差矩阵范数分析
covariance_matrix = [
[0.04, 0.02, 0.01, 0.015],
[0.02, 0.09, 0.03, 0.025],
[0.01, 0.03, 0.16, 0.020],
[0.015, 0.025, 0.020, 0.25]
]
cov_norm = norm(covariance_matrix, 2)
print(f"协方差矩阵谱范数: {cov_norm:.4f} (最大风险特征值)")
#协方差矩阵谱范数: 0.2619 (最大风险特征值)
风险向量分析
risk_contributions = [0.12, 0.18, 0.25, 0.15]
risk_norm = norm(risk_contributions)
print(f"风险贡献L1范数: {norm(risk_contributions, 1):.4f} (总风险)")
print(f"风险贡献L∞范数: {norm(risk_contributions, inf):.4f} (最大单个风险)")
#风险贡献L1范数: 0.7000 (总风险)
#风险贡献L∞范数: 0.2500 (最大单个风险)
示例8. 物理系统与能量计算
力学系统速度向量
velocity_vector = [3.2, -1.5, 4.1] # m/s
speed = norm(velocity_vector)
print(f"速度向量L2范数: {speed:.3f} m/s (速率)")
#速度向量L2范数: 5.413 m/s (速率)
假设质量1kg，计算动能
mass = 1.0
kinetic_energy = 0.5 * mass * float(speed) ** 2
print(f"动能: {kinetic_energy:.3f} J")
#动能: 14.650 J
应力张量范数
stress_tensor = [
[100, -20, 15],
[-20, 80, 25],
[15, 25, 120]
]
stress_norm = norm(stress_tensor, fro)
print(f"应力张量Frobenius范数: {stress_norm:.3f} MPa")
#应力张量Frobenius范数: 182.483 MPa
示例9. 网络分析与图论
邻接矩阵范数分析
adjacency_matrix = [
[0, 1, 1, 0],
[1, 0, 1, 1],
[1, 1, 0, 1],
[0, 1, 1, 0]
]
网络邻接矩阵范数分析:
adj_norms = {
"谱范数": norm(adjacency_matrix),
"L1范数": norm(adjacency_matrix, 1),
"L∞范数": norm(adjacency_matrix, inf)
}
for name, norm_val in adj_norms.items():
print(f"{name}: {norm_val}")
#谱范数: 3.16227766016838
#L1范数: 3.00000000000000
#L∞范数: 3.00000000000000
节点重要性向量
centrality_scores = [0.8, 1.2, 0.9, 1.1]
centrality_norm = norm(centrality_scores)
print(f"中心性得分L2范数: {centrality_norm:.3f}")
#中心性得分L2范数: 2.025
示例10. 数据压缩与降维
主成分分析中的奇异值
singular_values = [45.2, 23.1, 8.7, 3.2, 1.1]
sv_norm = norm(singular_values)
print(f"奇异值向量L2范数: {sv_norm:.3f}")
#奇异值向量L2范数: 51.612
计算能量保留比例
total_energy = float(sv_norm) ** 2
for k in range(1, len(singular_values) + 1):
partial_norm = norm(singular_values[:k])
energy_ratio = (float(partial_norm) ** 2) / total_energy
print(f"前{k}个主成分保留能量: {energy_ratio:.1%}")
#前1个主成分保留能量: 76.7%
#前2个主成分保留能量: 96.7%
#前3个主成分保留能量: 99.6%
#前4个主成分保留能量: 100.0%
#前5个主成分保留能量: 100.0%
数据重构误差
original_data_norm = 156.7
reconstructed_norm = 154.2
reconstruction_error = abs(original_data_norm - reconstructed_norm)
print(f"重构误差: {reconstruction_error:.3f}")
#重构误差: 2.500
示例11. 系统鲁棒性分析
system_matrix = [[0.5, 0.2], [0.1, 0.6]]
disturbance_vector = [0.3, -0.4]
系统增益分析
system_gain = norm(system_matrix)
disturbance_magnitude = norm(disturbance_vector)
最大可能输出
max_output = float(system_gain) * float(disturbance_magnitude)
print(f"系统增益 (谱范数): {system_gain:.4f}")
print(f"扰动幅度 (L2范数): {disturbance_magnitude:.4f}")
print(f"最大可能输出: {max_output:.4f}")
#系统增益 (谱范数): 0.8124
#扰动幅度 (L2范数): 0.5000
#最大可能输出: 0.4062
稳定性裕度
if float(system_gain) < 1:
robustness = "高鲁棒性"
elif float(system_gain) < 2:
robustness = "中等鲁棒性"
else:
robustness = "低鲁棒性"
print(f"系统鲁棒性: {robustness}")
#系统鲁棒性: 高鲁棒性
示例12. 比较不同范数的特性
data_matrix = [[1, -2, 3], [-4, 5, -6], [7, -8, 9]]
norms = {
"Frobenius": ('"fro"', "所有元素的总体大小"),
"谱范数 (L2)": ("2", "最大缩放能力"),
"列和范数 (L1)": ("1", "最大列和"),
"行和范数 (L∞)": ("inf", "最大行和")
}
for name, (norm_type, description) in norms.items():
norm_val = norm(data_matrix, norm_type)
print(f"{name:15}: {norm_val:8.4f} - {description}")
#Frobenius : 16.8819 - 所有元素的总体大小
#谱范数 (L2) : 16.8481 - 最大缩放能力
#列和范数 (L1) : 18.0000 - 最大列和
#行和范数 (L∞) : 24.0000 - 最大行和
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def norm_vector_matrix(input_str):
"""
对标MATLAB的norm函数，计算向量或矩阵的范数。
参数:
input_str: 输入字符串，格式示例:
- 矩阵: "[[1,2],[3,4]]"
- 带范数类型: "([1,2,3], 2)" 或 "([[1,2],[3,4]], 'fro')"
返回:
数值结果（浮点数） 或 错误信息字符串
"""
try:
# 解析输入字符串为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
v = expr[0]
p = str(expr[1])
if isinstance(v, list) and p in ['1', '2', 'inf', 'fro']:
matrix = sp.Matrix(v)
if p == '1':
result = matrix.norm(1).evalf()
elif p == '2':
result = matrix.norm(2).evalf()
elif p == 'inf':
result = matrix.norm(sp.oo).evalf()
elif p == 'fro':
result = matrix.norm('fro').evalf()
else:
error = True
elif isinstance(expr, list):
A = sp.Matrix(expr)
result = A.norm().evalf()
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误：{str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1: 向量默认范数（2范数）
print("示例1: norm('[1,2,3]') =", norm_vector_matrix('[1,2,3]'))
# 3.74165738677394
# 示例2: 矩阵默认范数（2范数，最大奇异值）
print("示例2: norm('[[1,2],[3,4]]') =", norm_vector_matrix('[[1,2],[3,4]]'))
# 5.47722557505166
# 示例3: 矩阵Frobenius范数
print("示例3: norm('([[1,2],[3,4]], \"fro\")') =",
norm_vector_matrix('([[1,2],[3,4]], "fro")'))
# 5.47722557505166
# 示例4: 向量无穷范数
print("示例4: norm('([1,-3,2], inf)') =", norm_vector_matrix('([1,-3,2], inf)'))
# 3.00000000000000
正态累积分布函数
p = normcdf(x) 返回标准正态分布的累积分布函数 (cdf)，在 x 中的值处计算函数值。
p = normcdf(x,mu) 返回具有均值 mu 和单位标准差的正态分布的 cdf，在 x 中的值处计算函数值。
p = normcdf(x,mu,sigma) 返回具有均值 mu 和标准差 sigma 的正态分布的 cdf，在 x 中的值处计算函数值。
x — 用于计算 cdf 的值,标量值,标量值组成的数组
mu — 均值,0 (默认),标量值,标量值组成的数组
sigma — 标准差,1 (默认),非负标量值,非负标量值组成的数组
p — cdf 值,标量值,标量值组成的数组
示例1. 质量控制与六西格玛
假设产品长度均值为100mm，标准差为2mm
spec_limits = [96, 98, 100, 102, 104]
mu, sigma = 100, 2
for limit in spec_limits:
cdf_value = normcdf(limit, mu, sigma)
defect_prob = 1 - float(cdf_value) if limit > mu else float(cdf_value)
print(f"尺寸≤{limit}mm的概率: {cdf_value:.4f} ({defect_prob:.2%}缺陷率)")
#尺寸≤96mm的概率: 0.0228 (2.28%缺陷率)
#尺寸≤98mm的概率: 0.1587 (15.87%缺陷率)
#尺寸≤100mm的概率: 0.5000 (50.00%缺陷率)
#尺寸≤102mm的概率: 0.8413 (15.87%缺陷率)
#尺寸≤104mm的概率: 0.9772 (2.28%缺陷率)
计算合格率（在98mm-102mm之间）
lower_cdf = normcdf(98, mu, sigma)
upper_cdf = normcdf(102, mu, sigma)
yield_rate = float(upper_cdf) - float(lower_cdf)
print(f"\n合格率(98-102mm): {yield_rate:.2%}")
#合格率(98-102mm): 68.27%
示例2. 考试成绩分析
假设考试均分72，标准差15
exam_mean = 72
exam_std = 15
grade_boundaries = [
(90, "A"), (80, "B"), (70, "C"),
(60, "D"), (0, "F")
]
各分数段概率分布:
prev_cdf = 0
for score, grade in sorted(grade_boundaries, reverse=True):
cdf_value = normcdf(score, exam_mean, exam_std)
prob_segment = float(cdf_value) - prev_cdf
print(f"{grade}等 (> {score}分): {prob_segment:.2%}")
prev_cdf = float(cdf_value)
#A等 (> 90分): 88.49%
#B等 (> 80分): -18.18%
#C等 (> 70分): -25.61%
#D等 (> 60分): -23.51%
#F等 (> 0分): -21.19%
录取分数线分析
admission_scores = [85, 90, 95]
录取分数线分析:
for score in admission_scores:
cdf_value = normcdf(score, exam_mean, exam_std)
admission_rate = 1 - float(cdf_value)
print(f"分数线{score}分以上的概率: {admission_rate:.2%}")
#分数线85分以上的概率: 19.31%
#分数线90分以上的概率: 11.51%
#分数线95分以上的概率: 6.26%
示例3. 医疗诊断与参考范围
假设血红蛋白正常值均值为14.5g/dL，标准差为1.2
hb_mean = 14.5
hb_std = 1.2
计算95%参考范围（均值±1.96标准差）
lower_ref = hb_mean - 1.96 * hb_std
upper_ref = hb_mean + 1.96 * hb_std
print(f"95%参考值范围: [{lower_ref:.1f}, {upper_ref:.1f}] g/dL")
#95%参考值范围: [12.1, 16.9] g/dL
异常值概率
abnormal_levels = [12.0, 13.0, 16.0, 17.0]
for level in abnormal_levels:
if level < hb_mean:
cdf_value = normcdf(level, hb_mean, hb_std)
prob = float(cdf_value)
status = "偏低"
else:
cdf_value = normcdf(f"({level}, {hb_mean}, {hb_std})")
prob = 1 - float(cdf_value)
status = "偏高"
print(f"血红蛋白{level}g/dL ({status})的概率: {prob:.4f} ({prob:.2%})")
#血红蛋白12.0g/dL (偏低)的概率: 0.0062 (0.62%)
#血红蛋白13.0g/dL (偏低)的概率: 0.0668 (6.68%)
#血红蛋白16.0g/dL (偏高)的概率: 0.0668 (6.68%)
#血红蛋白17.0g/dL (偏高)的概率: 0.0062 (0.62%)
示例4. 工程可靠性分析
假设设备寿命服从正态分布，均值5000小时，标准差800小时
life_mean = 5000
life_std = 800
warranty_periods = [1000, 2000, 3000, 4000, 5000]
不同保修期的故障概率:
for period in warranty_periods:
cdf_value = normcdf(period, life_mean, life_std)
failure_prob = float(cdf_value)
reliability = 1 - failure_prob
print(f"{period}小时保修期: 故障率{failure_prob:.2%}, 可靠度{reliability:.2%}")
#1000小时保修期: 故障率0.00%, 可靠度100.00%
#2000小时保修期: 故障率0.01%, 可靠度99.99%
#3000小时保修期: 故障率0.62%, 可靠度99.38%
#4000小时保修期: 故障率10.56%, 可靠度89.44%
#5000小时保修期: 故障率50.00%, 可靠度50.00%
计算可靠寿命（90%可靠度对应的寿命）
reliability_target = 0.90
from scipy.stats import norm
reliable_life = norm.ppf(1 - reliability_target, life_mean, life_std)
print(f"\n90%可靠度对应的寿命: {reliable_life:.0f}小时")
#90%可靠度对应的寿命: 3975小时
示例5. 气象数据分析
某地7月平均气温28°C，标准差3°C
temp_mean = 28
temp_std = 3
temperature_thresholds = [25, 30, 35, 40]
for temp in temperature_thresholds:
cdf_value = normcdf(temp, temp_mean, temp_std)
if temp > temp_mean:
prob = 1 - float(cdf_value)
desc = f"高于{temp}°C"
else:
prob = float(cdf_value)
desc = f"低于{temp}°C"
print(f"气温{desc}的概率: {prob:.4f} ({prob:.2%})")
#气温低于25°C的概率: 0.1587 (15.87%)
#气温高于30°C的概率: 0.2525 (25.25%)
#气温高于35°C的概率: 0.0098 (0.98%)
#气温高于40°C的概率: 0.0000 (0.00%)
极端天气事件
extreme_heat = 35 # 高温阈值
heat_prob = 1 - float(normcdf(extreme_heat, temp_mean, temp_std))
print(f"\n极端高温(≥{extreme_heat}°C)概率: {heat_prob:.2%}")
print(f"预计每年出现{heat_prob * 365:.1f}天")
#极端高温(≥35°C)概率: 0.98%
#预计每年出现3.6天
示例6. 生产调度与交付时间
交付时间均值为10天，标准差2天
delivery_mean = 10
delivery_std = 2
deadlines = [7, 8, 9, 10, 11, 12, 14]
按时交付概率分析:
for days in deadlines:
cdf_value = normcdf(days, delivery_mean, delivery_std)
on_time_prob = float(cdf_value)
late_prob = 1 - on_time_prob
print(f"{days}天交付: 按时概率{on_time_prob:.2%}, 延迟概率{late_prob:.2%}")
#7天交付: 按时概率6.68%, 延迟概率93.32%
#8天交付: 按时概率15.87%, 延迟概率84.13%
#9天交付: 按时概率30.85%, 延迟概率69.15%
#10天交付: 按时概率50.00%, 延迟概率50.00%
#11天交付: 按时概率69.15%, 延迟概率30.85%
#12天交付: 按时概率84.13%, 延迟概率15.87%
#14天交付: 按时概率97.72%, 延迟概率2.28%
服务水平协议(SLA)分析
sla_targets = [0.90, 0.95, 0.99]
SLA目标对应的交付时间:
for sla in sla_targets:
delivery_time = norm.ppf(sla, delivery_mean, delivery_std)
print(f"{sla:.0%} SLA: {delivery_time:.1f}天")
#90% SLA: 12.6天
#95% SLA: 13.3天
#99% SLA: 14.7天
示例7. 符号计算与理论分析
符号表达式
symbolic_expr = normcdf(x)
print(f"标准正态分布CDF: {symbolic_expr}")
#标准正态分布CDF: 0.5*erf(sqrt(2)*x/2) + 0.5
带参数的符号表达式
symbolic_with_params = normcdf(y, mu, sigma)
print(f"一般正态分布CDF: {symbolic_with_params}")
#一般正态分布CDF: 0.5*erf(sqrt(2)*(-mu + y)/(2*sigma)) + 0.5
特定点的符号计算
specific_point = normcdf(a+b, c, d)
print(f"线性变换CDF: {specific_point}")
#线性变换CDF: 0.5*erf(sqrt(2)*(a + b - c)/(2*d)) + 0.5
示例8. 批量数据处理
考试成绩矩阵
scores_matrix = [[65, 72, 80], [85, 90, 78], [92, 88, 95]]
exam_mean = 75
exam_std = 10
计算每个分数的累积概率
for row in scores_matrix:
prob_row = []
for score in row:
cdf_val = normcdf(score, exam_mean, exam_std)
prob_row.append(f"{float(cdf_val):.3f}")
print(f"[{', '.join(prob_row)}]")
#[0.159, 0.382, 0.691]
#[0.841, 0.933, 0.618]
#[0.955, 0.903, 0.977]
向量化计算
score_vector = [60, 70, 80, 90, 100]
vector_result = []
for score in score_vector:
cdf_val = normcdf(score, exam_mean, exam_std)
vector_result.append(float(cdf_val))
print(f"\n分数向量: {score_vector}")
print(f"累积概率: {[f'{x:.3f}' for x in vector_result]}")
#分数向量: [60, 70, 80, 90, 100]
#累积概率: ['0.067', '0.309', '0.691', '0.933', '0.994']
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy.stats import norm
import sympy as sp
from sympy.stats import Normal, cdf
def normal_distribution_cdf(input_str):
"""
对标Matlab的normcdf函数，计算正态分布的累积分布函数（CDF）。
参数:
input_str: 输入表达式字符串，可以是标量、矩阵（列表形式）、元组或符号。
返回:
正态分布的CDF值，类型可以是标量、矩阵或符号表达式。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def normal_cdf_sym(x_val, mu=0, sigma=1):
"""计算正态分布的CDF，处理数值和符号输入。"""
W = Normal("W", mu, sigma)
cdf_expr = cdf(W)(x_val)
# 仅当所有参数均为数值时求值
return cdf_expr.evalf()
if isinstance(expr, tuple):
if len(expr) == 3:
x, mu, sigma = expr
elif len(expr) == 2:
x, mu = expr
sigma = 1
else:
error = True
if all(e.is_number for e in expr):
result = norm.cdf(int(x), loc=float(mu), scale=float(sigma))
else:
result = normal_cdf_sym(x, mu, sigma)
elif expr.is_number:
x = expr
mu = 0
sigma = 1
result = norm.cdf(float(x), loc=mu, scale=sigma)
elif expr.free_symbols:
# 输入是标量或符号
result = normal_cdf_sym(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 标量输入
print(normal_distribution_cdf("0"))
# 0.5
# 示例2: 符号输入
print(normal_distribution_cdf("x"))
# 0.5*erf(sqrt(2)*x/2) + 0.5
# 数值计算
print(normal_distribution_cdf("0"))
# 0.5
print(normal_distribution_cdf("1, 0, 2"))
# 0.6914624612740131
2-范数估值
n = normest(S) 返回矩阵 S 的 2-范数估值.
本函数旨在主要用于稀疏矩阵,尽管也适用于大型满矩阵.
S — 输入矩阵,稀疏矩阵,满矩阵.
tol — 相对误差容限,1e-6 (默认),非负实数
示例1: 图像处理中的变换矩阵
45度旋转矩阵
theta = pi/4
rotation_matrix = [
[cos(theta), -sin(theta)],
[sin(theta), cos(theta)]
]
result = normest(rotation_matrix)
print(f"旋转矩阵范数: {result}")
#旋转矩阵范数: 1.0
#因为旋转不改变向量长度
示例2: 控制系统中的状态转移矩阵
离散系统的状态转移矩阵
system_matrix = [
[0.8, 0.2],
[-0.1, 0.9]
]
result = normest(system_matrix)
print(f"系统矩阵范数: {result}")
print(f"系统稳定性判断: {'稳定' if result < 1 else '可能不稳定'}")
#系统矩阵范数: 0.9338433171688098
#系统稳定性判断: 稳定
示例3: 机器学习中的权重矩阵
模拟一个神经网络的权重矩阵
weights = [
[0.1, -0.2, 0.3, -0.1],
[0.4, 0.2, -0.3, 0.5],
[-0.2, 0.1, 0.4, -0.3]
]
result = normest(weights)
print(f"权重矩阵范数: {result}")
#权重矩阵范数: 0.9046719563003617
#范数可用于梯度裁剪和正则化
示例4: 结构力学中的刚度矩阵
简单的弹簧系统刚度矩阵
stiffness_matrix = [
[2, -1, 0],
[-1, 2, -1],
[0, -1, 1]
]
result = normest(stiffness_matrix)
print(f"刚度矩阵范数: {result}")
#刚度矩阵范数: 3.246978619287783
#范数反映系统的刚性程度
示例5: 量子力学中的幺正矩阵
Hadamard 门
hadamard = [
[1/sqrt(2), 1/sqrt(2)],
[1/sqrt(2), -1/sqrt(2)]
]
result = normest(hadamard)
print(f"Hadamard门范数: {result}")
#Hadamard门范数: 0.9999999999999999
示例6: 大规模稀疏矩阵的近似（使用小规模模拟）
对角占优矩阵，模拟某些物理问题
large_matrix = [
[5, 1, 0, 0],
[1, 5, 1, 0],
[0, 1, 5, 1],
[0, 0, 1, 5]
]
exact_norm = np.linalg.norm(large_matrix, 2)
estimated_norm = normest(large_matrix)
print(f"精确范数: {exact_norm}")
print(f"估计范数: {estimated_norm}")
print(f"相对误差: {abs(exact_norm - estimated_norm)/exact_norm*100:.6f}%")
#精确范数: 6.618033988749895
#估计范数: 6.618032825998585
#相对误差: 0.000018%
示例7: 条件数估计
A = [[1, 2], [3, 4]]
A_norm = normest(A)
使用相同的算法估计逆矩阵的范数
A_inv = np.linalg.inv(np.array(A))
A_inv_norm = normest(A_inv)
condition_estimate = A_norm * A_inv_norm
exact_condition = np.linalg.cond(A, 2)
print(f"条件数估计: {condition_estimate}")
print(f"精确条件数: {exact_condition}")
#条件数估计: 14.933032778908295
#精确条件数: 14.933034373659268
示例8: 数值方法的收敛性分析
Jacobi 方法的迭代矩阵
A = [[4, 1, 1], [1, 4, 1], [1, 1, 4]]
D_inv = [[1/4, 0, 0], [0, 1/4, 0], [0, 0, 1/4]]
迭代矩阵 T = I - D^{-1}A
T = [[1 - 1/4*4, -1/4*1, -1/4*1],
[-1/4*1, 1 - 1/4*4, -1/4*1],
[-1/4*1, -1/4*1, 1 - 1/4*4]]
spectral_radius_est = normest(T)
print(f"迭代矩阵谱半径估计: {spectral_radius_est}")
print(f"收敛性: {'收敛' if spectral_radius_est < 1 else '发散'}")
#迭代矩阵谱半径估计: 0.4999999261708375
#收敛性: 收敛
示例9: 数据科学中的协方差矩阵
模拟三个变量的协方差矩阵
covariance = [
[2.5, 0.8, -0.3],
[0.8, 1.2, 0.5],
[-0.3, 0.5, 0.9]
]
result = normest(covariance)
print(f"协方差矩阵范数: {result}")
#协方差矩阵范数: 2.8824996891697734
#反映数据的主要变化方向
示例10: 对比不同容差的影响
test_matrix = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
exact = np.linalg.norm(test_matrix, 2)
for tol in [1e-2, 1e-4, 1e-6, 1e-8]:
estimated = normest(test_matrix, tol)
error = abs(exact - estimated)
print(f"容差 {tol:.0e}: 估计值={estimated:.8f}, 误差={error:.2e}")
#容差 1e-02: 估计值=16.84808732, 误差=1.60e-05
#容差 1e-04: 估计值=16.84810335, 误差=3.04e-09
#容差 1e-06: 估计值=16.84810335, 误差=1.23e-09
#容差 1e-08: 估计值=16.84810334, 误差=8.95e-09
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def norm_estimate_matrix(input_str):
"""
对标MATLAB的normest函数，估计矩阵的2-范数（最大奇异值）。
参数:
input_str: 输入字符串，支持以下格式:
- 仅矩阵: "[[1,2],[3,4]]"
- 矩阵+容差: "([[1,2],[3,4]], 1e-6)"
返回:
数值估计结果（浮点数） 或 错误信息字符串
"""
try:
def eval_normest(A_np, tol=1e-6, max_iter=100):
# 如果输入是向量，直接返回其2-范数
if A_np.ndim == 1 or A_np.shape[1] == 1:
return np.linalg.norm(A_np, 2)
# 以下是矩阵的幂迭代算法
b = np.random.randn(A_np.shape[1])
b /= np.linalg.norm(b)
current_norm = 0
for _ in range(max_iter):
Ab = A_np @ b
new_b = A_np.T @ Ab
new_norm = np.linalg.norm(Ab)
new_b_norm = np.linalg.norm(new_b)
if new_b_norm == 0:
break
b = new_b / new_b_norm
if abs(new_norm - current_norm) < tol:
break
current_norm = new_norm
return current_norm
# 解析输入
expr = sp.sympify(input_str)
A_sym = None
tol = 1e-6
# 处理输入格式：可能是 (矩阵, tol) 或 仅矩阵
if isinstance(expr, tuple) and len(expr) == 2:
A_sym = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
tol = float(expr[1])
else:
A_sym = sp.Matrix(expr) if isinstance(expr, list) else None
# 验证矩阵有效性
if A_sym is None:
return "错误：输入不是有效矩阵"
if A_sym.is_symbolic():
return "错误：矩阵包含非数值元素"
# 转换为NumPy数组
A_np = np.array(A_sym.tolist(), dtype=float)
# 调用核心算法
# 展平为向量（如果是列向量）
if A_np.shape[1] == 1:
A_np = A_np.flatten()
return eval_normest(A_np, tol)
except Exception as e:
return f"错误：{str(e)}"
# 示例代码
if __name__ == "__main__":
# 示例1: 小型稠密矩阵
print("示例1: normest([[3,0],[0,4]]) =",
norm_estimate_matrix("[[3,0],[0,4]]"))
# 3.999999438364604
# 示例2: 自定义容差
print("示例2: normest(([[1,2],[3,4]], 1e-8)) =",
norm_estimate_matrix("([[1,2],[3,4]], 1e-8)"))
# 5.4649857042187255
# 示例3: 一维输入
print("示例3: normest([1,2,3]) =",
norm_estimate_matrix("[1,2,3]"))
# 3.7416573867739413
y = normpdf(x) 返回标准正态分布的概率密度函数 (pdf)，在 x 中的值处计算函数值。
y = normpdf(x,mu) 返回具有均值 mu 和单位标准差的正态分布的 pdf，在 x 中的值处计算函数值。
y = normpdf(x,mu,sigma) 返回具有均值 mu 和标准差 sigma 的正态分布的 pdf，在 x 中的值处计算函数值。
x — 用于计算 pdf 的值,标量值 | 标量值组成的数组
mu — 均值, 0 (默认) | 标量值 | 标量值组成的数组
sigma — 标准差,1 (默认) | 正标量值 | 正标量值组成的数组
y — pdf 值,标量值 | 标量值组成的数组
示例1: 质量控制 - 产品尺寸分布
假设某零件长度服从正态分布，均值10mm，标准差0.2mm
lengths = [9.5, 9.8, 10.0, 10.2, 10.5]
for length in lengths:
prob = normpdf(length, 10, 0.2)
print(f"长度 {length}mm 的概率密度: {prob}")
#长度 9.5mm 的概率密度: 错误: 标准差sigma必须为正数
#长度 9.8mm 的概率密度: 错误: 标准差sigma必须为正数
#长度 10.0mm 的概率密度: 错误: 标准差sigma必须为正数
#长度 10.2mm 的概率密度: 错误: 标准差sigma必须为正数
#长度 10.5mm 的概率密度: 错误: 标准差sigma必须为正数
示例2: 考试成绩分析
某考试平均分75，标准差10
scores = [55, 65, 75, 85, 95]
for score in scores:
prob = normpdf(score, 75, 10)
print(f"分数 {score} 的概率密度: {prob:.6f}")
#分数 55 的概率密度: 0.005399
#分数 65 的概率密度: 0.024197
#分数 75 的概率密度: 0.039894
#分数 85 的概率密度: 0.024197
#分数 95 的概率密度: 0.005399
示例3: 金融风险管理 - 资产收益率
假设某股票日收益率服从正态分布，均值0.001，标准差0.02
returns = [-0.05, -0.02, 0.0, 0.02, 0.05]
for ret in returns:
prob = normpdf(ret, 0.001, 2)
print(f"收益率 {ret:.3f} 的概率密度: {prob}")
#收益率 -0.050 的概率密度: 0.1994062976877908
#收益率 -0.020 的概率密度: 0.19946014465718026
#收益率 0.000 的概率密度: 0.1994711152668254
#收益率 0.020 的概率密度: 0.19946213926859985
#收益率 0.050 的概率密度: 0.19941128290754798
示例4: 医学研究 - 生理指标
成人血压的舒张压，均值80mmHg，标准差10mmHg
blood_pressure = [60, 70, 80, 90, 100]
for bp in blood_pressure:
prob = normpdf(bp, 80, 10)
print(f"舒张压 {bp}mmHg 的概率密度: {prob:.6f}")
#舒张压 60mmHg 的概率密度: 0.005399
#舒张压 70mmHg 的概率密度: 0.024197
#舒张压 80mmHg 的概率密度: 0.039894
#舒张压 90mmHg 的概率密度: 0.024197
#舒张压 100mmHg 的概率密度: 0.005399
示例5: 工程测量 - 误差分析
测量误差服从均值为0的正态分布，标准差为测量精度
measurement_errors = [-0.3, -0.1, 0.0, 0.1, 0.3]
for error in measurement_errors:
prob = normpdf(error, 0, 1)
print(f"测量误差 {error} 的概率密度: {prob}")
#测量误差 -0.3 的概率密度: 0.38138781546052414
#测量误差 -0.1 的概率密度: 0.3969525474770118
#测量误差 0.0 的概率密度: 0.3989422804014327
#测量误差 0.1 的概率密度: 0.3969525474770118
#测量误差 0.3 的概率密度: 0.38138781546052414
示例6: 批量计算 - 多个点的PDF值
x_values = [-2, -1, 0, 1, 2]
result = normpdf(x_values)
print("标准正态分布在 x =", x_values, "处的PDF值:")
print([f"{p:.6f}" for p in result])
#标准正态分布在 x = [-2, -1, 0, 1, 2] 处的PDF值:
#['0.053991', '0.241971', '0.398942', '0.241971', '0.053991']
示例7: 不同参数的正态分布比较
x = 1.0
distributions = [
("N(0,1)", 0, 1),
("N(0,2)", 0, 2),
("N(1,1)", 1, 1),
("N(1,2)", 1, 2)
]
for name, mu, sigma in distributions:
prob = normpdf(x, mu, sigma)
print(f"{name} 在 x={x} 处的PDF: {prob:.6f}")
#N(0,1) 在 x=1.0 处的PDF: 0.241971
#N(0,2) 在 x=1.0 处的PDF: 0.176033
#N(1,1) 在 x=1.0 处的PDF: 0.398942
#N(1,2) 在 x=1.0 处的PDF: 0.199471
示例8: 符号计算 - 得到PDF的解析表达式
一般正态分布的PDF表达式
pdf_expr = normpdf(x, mu, sigma)
print(pdf_expr)
#0.398942280401433*exp(-(-mu + x)**2/(2*sigma**2))/sigma
示例9: 矩阵输入 - 处理二维数据
模拟二维测量数据
measurement_grid = [[0.5, 1.0, 1.5],
[2.0, 2.5, 3.0]]
result_matrix = normpdf(measurement_grid)
测量网格的PDF值矩阵:
print(result_matrix)
#[[0.35206533,0.24197072,0.1295176 ]
[0.05399097,0.0175283,0.00443185]]
示例10: 假设检验中的似然比计算
在均值检验中计算不同假设下的似然
observed_data = 1.96 # 常见的临界值
零假设 H0: mu=0, sigma=1
likelihood_h0 = normpdf(observed_data, 0, 1)
备择假设 H1: mu=2, sigma=1
likelihood_h1 = normpdf(observed_data, 2, 1)
likelihood_ratio = likelihood_h0 / likelihood_h1
print(f"H0下的似然: {likelihood_h0:.6f}")
print(f"H1下的似然: {likelihood_h1:.6f}")
print(f"似然比: {likelihood_ratio:.4f}")
#H0下的似然: 0.058441
#H1下的似然: 0.398623
#似然比: 0.1466
示例11: 贝叶斯统计中的先验分布
使用正态分布作为参数先验
parameter_values = [-1, 0, 1]
prior_mean = 0
prior_std = 1
for param in parameter_values:
prior_prob = normpdf(param, prior_mean, prior_std)
print(f"参数值 {param} 的先验概率密度: {prior_prob:.6f}")
#参数值 -1 的先验概率密度: 0.241971
#参数值 0 的先验概率密度: 0.398942
#参数值 1 的先验概率密度: 0.241971
示例12: 信号处理中的噪声模型
高斯白噪声的概率密度
noise_levels = [-3, -1.5, 0, 1.5, 3]
for noise in noise_levels:
prob = normpdf(noise, 0, 1)
print(f"噪声幅度 {noise} 的概率密度: {prob:.6f}")
#噪声幅度 -3 的概率密度: 0.004432
#噪声幅度 -1.5 的概率密度: 0.129518
#噪声幅度 0 的概率密度: 0.398942
#噪声幅度 1.5 的概率密度: 0.129518
#噪声幅度 3 的概率密度: 0.004432
示例13: 机器学习中的特征分布
假设某个特征服从正态分布
feature_values = [-2, -1, 0, 1, 2]
feature_mean = 0.5
feature_std = 1.2
for value in feature_values:
prob = normpdf(value, feature_mean, feature_std)
print(f"特征值 {value} 的概率密度: {prob:.6f}")
#特征值 -2 的概率密度: 0.017528
#特征值 -1 的概率密度: 0.129518
#特征值 0 的概率密度: 0.352065
#特征值 1 的概率密度: 0.352065
#特征值 2 的概率密度: 0.129518
示例14: 经济学中的收入分布
对数收入通常假设为正态分布
log_incomes = [10, 11, 12, 13, 14] # 对数收入
mean_log_income = 11.5
std_log_income = 1.0
for log_income in log_incomes:
prob = normpdf(log_income, mean_log_income, std_log_income)
print(f"对数收入 {log_income} 的概率密度: {prob:.6f}")
#对数收入 10 的概率密度: 0.129518
#对数收入 11 的概率密度: 0.352065
#对数收入 12 的概率密度: 0.352065
#对数收入 13 的概率密度: 0.129518
#对数收入 14 的概率密度: 0.017528
示例15: 验证数值精度
与精确值比较
test_points = [0, 1, 2, 3]
exact_values = [0.39894228, 0.24197072, 0.05399097, 0.00443185]
print("x\t我们的实现\t精确值\t\t误差")
print("-" * 50)
for x, exact in zip(test_points, exact_values):
our_value = normpdf(x)
error = abs(our_value - exact)
print(f"{x}\t{our_value:.8f}\t{exact:.8f}\t{error:.2e}")
#x 我们的实现 精确值 误差
#--------------------------------------------------
#0 0.39894228 0.39894228 4.01e-10
#1 0.24197072 0.24197072 4.52e-09
#2 0.05399097 0.05399097 3.49e-09
#3 0.00443185 0.00443185 1.59e-09
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
from scipy.stats import norm
import sympy as sp
from sympy.stats import Normal, density
def normal_distribution_pdf(input_str):
"""
对标Matlab的normcdf函数，计算正态分布的累积分布函数（CDF）。
参数:
input_str: 输入表达式字符串，可以是标量、矩阵（列表形式）、元组或符号。
返回:
正态分布的CDF值，类型可以是标量、矩阵或符号表达式。如果输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def normal_pdf_sym(x_val, mu=0, sigma=1):
"""计算正态分布的CDF，处理数值和符号输入。"""
W = Normal("W", mu, sigma)
pdf_expr = density(W)(x_val)
# 仅当所有参数均为数值时求值
return pdf_expr.evalf()
if isinstance(expr, tuple):
if len(expr) == 3:
x, mu, sigma = expr
elif len(expr) == 2:
x, mu = expr
sigma = 1
else:
error = True
if all(e.is_number for e in expr):
result = norm.pdf(int(x), loc=float(mu), scale=float(sigma))
else:
result = normal_pdf_sym(x, mu, sigma)
elif expr.is_number:
x = expr
mu = 0
sigma = 1
result = norm.cdf(float(x), loc=mu, scale=sigma)
elif expr.free_symbols:
# 输入是标量或符号
result = normal_pdf_sym(expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例用法
if __name__ == "__main__":
# 示例1: 标量输入
print(normal_distribution_pdf("0"))
# 0.5
# 示例2: 符号输入
print(normal_distribution_pdf("x"))
# 0.398942280401433*exp(-x**2/2)
# 数值计算
print(normal_distribution_pdf("0"))
# 0.5
print(normal_distribution_pdf("1, 0, 2"))
# 0.17603266338214976
质因子及其指数
nPrimes(n)返回质因子及对应的指数
n — 正整数
示例1: 密码学 - RSA加密
RSA模数通常是两个大质数的乘积
rsa_numbers = [15, 35, 77, 143] # 小数字用于演示
for num in rsa_numbers:
factors = nPrimes(num)
print(f"RSA模数 {num} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
if len(factors) == 2 and all(exp == 1 for exp in factors.values()):
print(" ✓ 是有效的RSA模数（两个质数乘积）")
#RSA模数 15 = 3^1 × 5^1
# ✓ 是有效的RSA模数（两个质数乘积）
#RSA模数 35 = 5^1 × 7^1
# ✓ 是有效的RSA模数（两个质数乘积）
#RSA模数 77 = 7^1 × 11^1
# ✓ 是有效的RSA模数（两个质数乘积）
#RSA模数 143 = 11^1 × 13^1
# ✓ 是有效的RSA模数（两个质数乘积）
示例2: 数论研究 - 完全数
完全数是等于其真因子和的数
perfect_numbers = [6, 28, 496, 8128]
for num in perfect_numbers:
factors = nPrimes(num)
print(f"完全数 {num} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
#完全数 6 = 2^1 × 3^1
#完全数 28 = 2^2 × 7^1
#完全数 496 = 2^4 × 31^1
#完全数 8128 = 2^6 × 127^1
示例3: 计算机科学 - 哈希表大小
哈希表大小通常选择质数或质数的乘积以减少冲突
hash_sizes = [16, 31, 64, 101, 1009]
for size in hash_sizes:
factors = nPrimes(size)
if len(factors) == 1 and list(factors.values())[0] == 1:
print(f"哈希表大小 {size} 是质数 ✓")
else:
print(f"哈希表大小 {size} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
#哈希表大小 16 = 2^4
#哈希表大小 31 是质数 ✓
#哈希表大小 64 = 2^6
#哈希表大小 101 是质数 ✓
#哈希表大小 1009 是质数 ✓
示例4: 音乐理论 - 音阶比例
音程频率比可以用质因数表示
musical_intervals = {
"八度": 2,
"纯五度": 3 / 2,
"纯四度": 4 / 3,
"大三度": 5 / 4
}
for interval, ratio in musical_intervals.items():
# 将分数转换为整数表示
numerator = ratio.as_integer_ratio()[0] if hasattr(ratio, 'as_integer_ratio') else ratio.numerator
denominator = ratio.as_integer_ratio()[1] if hasattr(ratio, 'as_integer_ratio') else ratio.denominator
num_factors = nPrimes(numerator)
den_factors = nPrimes(denominator)
print(f"{interval}: {numerator}/{denominator} = "
f"{' × '.join(f'{p}^{e}' for p, e in num_factors.items())} / "
f"{' × '.join(f'{p}^{e}' for p, e in den_factors.items())}")
#八度: 2/1 = 2^1 /
#纯五度: 3/2 = 3^1 / 2^1
#纯四度: 6004799503160661/4503599627370496 = 3^3 × 7^1 × 19^1 × 73^1 × 87211^1 × 262657^1 / 2^52
#大三度: 5/4 = 5^1 / 2^2
示例5: 金融数学 - 复利计算
年利率转换为日利率等
interest_periods = [365, 52, 12, 4, 2] # 年化周期
for period in interest_periods:
factors = nPrimes(period)
print(f"{period} 次复利/年 = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
#365 次复利/年 = 5^1 × 73^1
#52 次复利/年 = 2^2 × 13^1
#12 次复利/年 = 2^2 × 3^1
#4 次复利/年 = 2^2
#2 次复利/年 = 2^1
示例6: 编码理论 - 错误检测
CRC生成多项式通常用质因数表示
crc_polynomials = [3, 5, 9, 17, 33] # 简化的多项式值
for poly in crc_polynomials:
factors = nPrimes(poly)
print(f"CRC多项式 {poly} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
#CRC多项式 3 = 3^1
#CRC多项式 5 = 5^1
#CRC多项式 9 = 3^2
#CRC多项式 17 = 17^1
#CRC多项式 33 = 3^1 × 11^1
示例7: 数论趣味 - 平方数和立方数
special_numbers = [36, 64, 100, 144, 216, 1000]
for num in special_numbers:
factors = nPrimes(num)
is_square = all(exp % 2 == 0 for exp in factors.values())
is_cube = all(exp % 3 == 0 for exp in factors.values())
properties = []
if is_square:
properties.append("平方数")
if is_cube:
properties.append("立方数")
prop_str = " (" + ", ".join(properties) + ")" if properties else ""
print(f"{num} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}{prop_str}")
#36 = 2^2 × 3^2 (平方数)
#64 = 2^6 (平方数, 立方数)
#100 = 2^2 × 5^2 (平方数)
#144 = 2^4 × 3^2 (平方数)
#216 = 2^3 × 3^3 (立方数)
#1000 = 2^3 × 5^3 (立方数)
示例8: 物理常数比值的质因数
某些物理常数比值包含有趣的质因数
physical_ratios = {
"精细结构常数倒数": 137,
"圆周率近似": 355, # 355/113 ≈ π
"自然对数底数平方": 7 # e² ≈ 7.389
}
for name, value in physical_ratios.items():
factors = nPrimes(value)
print(f"{name} {value} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
#精细结构常数倒数 137 = 137^1
#圆周率近似 355 = 5^1 × 71^1
#自然对数底数平方 7 = 7^1
示例9: 生日悖论相关数字
combinatorial_numbers = [365, 52, 12, 7] # 年天数、周数、月数、星期数
for num in combinatorial_numbers:
factors = nPrimes(num)
print(f"{num} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
#365 = 5^1 × 73^1
#52 = 2^2 × 13^1
#12 = 2^2 × 3^1
#7 = 7^1
示例10: 验证函数对大数的处理
large_numbers = [1001, 1024, 123456, 999999]
for num in large_numbers:
factors = nPrimes(num)
print(f"{num} = {' × '.join(f'{p}^{e}' for p, e in factors.items())}")
#1001 = 7^1 × 11^1 × 13^1
#1024 = 2^10
#123456 = 2^6 × 3^1 × 643^1
#999999 = 3^3 × 7^1 × 11^1 × 13^1 × 37^1
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
def number_prime_factors(input_str):
"""
此函数用于将输入字符串解析为整数，并计算该整数的质因数分解结果。
参数:
input_str (str): 输入的字符串，该字符串需能被解析为正整数。
返回:
dict: 若输入合法，返回一个字典，字典的键为质因数，值为该质因数在分解中的指数。
str: 若输入不合法，返回包含错误信息的字符串。
"""
try:
# 解析输入字符串为数值
n = eval(input_str)
# 检查输入是否为整数
if isinstance(n, float):
if not n.is_integer():
return "错误：输入必须为整数"
n = int(n)
elif not isinstance(n, int):
return "错误：输入必须为整数"
# 检查输入是否为正整数
if n < 1:
return "错误：输入必须为正整数"
factors = {} # 存储质因数及其指数
# 处理所有2的因子，使n变为奇数
while n % 2 == 0:
factors[2] = factors.get(2, 0) + 1
n = n // 2
# 处理奇数因子，从3开始，每次递增2
i = 3
max_factor = int(n ** 0.5) # 初始最大因子为当前n的平方根
while i <= max_factor and n > 1:
while n % i == 0:
factors[i] = factors.get(i, 0) + 1
n = n // i
max_factor = int(n ** 0.5) # 更新最大因子为新的n的平方根
i += 2 # 跳过偶数
# 若剩余n为质数且大于2
if n > 1:
factors[n] = 1
return factors
except Exception as e:
return f"错误: {e}"
print(number_prime_factors("200"))
# {2: 2, 5: 2}
print(number_prime_factors("12345"))
# {3: 1, 5: 1, 823: 1}
实数的第n次实根
Y = nthroot(X,N)返回X的第n次实根.如果X中的元素为负数,则N必须为奇数.
X — 输入数组,标量,向量,矩阵
N — 要计算的根,标量,向量,矩阵
示例1: 基本的数值计算
test_cases = [
"-8, 3", # -8的立方根
"16, 4", # 16的4次方根
"27, 3", # 27的立方根
"-32, 5", # -32的5次方根
"9, 2", # 9的平方根
]
for case in test_cases:
result = nthroot(case)
print(f"nthroot{case} = {result}")
#nthroot(-8, 3) = -2.0
#nthroot(16, 4) = 2.0
#nthroot(27, 3) = 3.0
#nthroot(-32, 5) = -2.0
#nthroot(9, 2) = 3.0
示例2: 符号计算
symbolic_cases = [
"x, 3",
"y**2, 4"
]
for case in symbolic_cases:
result = nthroot(case)
print(f"nthroot{case} = {result}")
#nthroot(x, 3) = x**(1/3)
#nthroot(y**2, 4) = (y**2)**(1/4)
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from sympy.simplify.simplify import nthroot
def nth_root_surds(input_str):
"""
MATLAB nthroot函数的Python实现，支持矩阵广播机制
参数格式：
"(x_matrix, n_matrix)" 或 "x_matrix"
特殊行为：
- 当x是行向量(1×N)、n是列向量(M×1)时，生成M×N矩阵
- 自动处理标量与矩阵的混合运算
- 严格校验负数偶次根情况
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def nthroot_sci(x, n):
"""
计算实数的第N次实根（对标MATLAB的nthroot函数）- 标量优化版
参数:
x (float): 输入实数
n (float): 根的次数
返回:
float: 对应的实根，或NaN（当实根不存在时）
说明:
- 此版本专门优化标量输入，更简洁高效
- 对于数组输入，请使用向量化版本
"""
# 处理x=0的情况
if x == 0:
if n == 0:
return np.nan # 0^0未定义
return 0.0 if n > 0 else np.inf # 0的正数次根为0，负数次根为无穷大
# 处理n=0的情况
if n == 0:
return np.nan # 任何数的0次根未定义
# 处理x为正数的情况
if x > 0:
return x ** (1 / n)
# 处理x为负数的情况
if n % 2 == 1: # 奇数次根
return -(-x) ** (1 / n)
# 其他情况（负数开偶次根）无实数解
return np.nan
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e) for e in expr)
result = nthroot_sci(*params)
else:
result = nthroot(*expr)
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 测试用例 ------------------------------------------------------------------
if __name__ == "__main__":
# 原题测试案例
test_input = "(-5,3)"
print(f"输入: {test_input}")
print("输出:")
print(nth_root_surds(test_input))
# -1.7099759466766968
test_input = "x,3"
print(f"输入: {test_input}")
print("输出:")
print(nth_root_surds(test_input))
# x**(1/3)
n次方根
nRoot(z,n)返回底数的根指数次方根的符号化结果
z — 底数
n — 根指数
示例1. 几何学应用
立方体边长计算
result = nRoot(64, 3) # 体积为64的立方体
print(result)
#应用: 立方体体积求边长: V = a³ → a = ∛V
#数值结果: 4.0
球体半径计算
result = nRoot(36@pi, 3) # 体积为36π的球体
print(result)
#应用: 球体体积求半径: V = (4/3)πr³ → r = ∛(3V/4π)
#数值结果: 4.83597586204941
示例2. 物理学应用
简谐运动周期
result = nRoot(4*@pi**2, 2) # T = 2π√(m/k) 中的部分计算
print(result)
#应用: 简谐运动周期计算
#数值结果: 6.28318530717959
声强级计算
result = nRoot(1000, 10) # 声强级的对数关系
print(result)
#应用: 声学: 声强级计算
#数值结果: 1.99526231496888
示例3. 金融学应用
年化收益率计算
result = nRoot(1.331, 3) # 3年总收益1.331倍
print(result)
#应用: 年化收益率: (1 + r)³ = 1.331 → r = ∛1.331 - 1
#数值结果: 1.10
复利计算
result = nRoot(2, 10) # 10年翻倍所需的年化收益率
print(result)
#应用: 复利翻倍时间计算
#数值结果: 1.07177346253629
示例4. 工程学应用
圆管直径计算
result = nRoot(16/@pi, 2) # 截面积为4的圆管
print(result)
#应用: 圆管截面积求直径: A = πr² → r = √(A/π)
#数值结果: 2.25675833419103
梁的截面模量
result = nRoot(27, 3) # 立方关系
print(result)
#应用: 结构工程: 梁的强度计算
#数值结果: 3.0
示例5. 统计学应用
几何平均数
result = nRoot(8*27*64, 3) # 三个数的几何平均
print(result)
#应用: 几何平均数: ∛(8×27×64)
#数值结果: 24.0
标准差计算
result = nRoot(16, 2) # 方差为16
print(result)
#应用: 标准差: σ = √方差
#数值结果: 4.0
示例6. 化学应用
浓度计算
result = nRoot(0.001, 3) # 稀释计算
print(result)
#应用: 化学浓度稀释计算
#数值结果: 0.10
反应速率
result = nRoot(8, 3) # 立方关系
print(result)
#应用: 化学反应速率计算
#数值结果: 2.0
示例7. 电子工程应用
RMS电压计算
result = nRoot(100, 2) # 峰值电压10V的RMS值
print(result)
#应用: RMS电压: V_rms = V_peak/√2
#数值结果: 10.0
阻抗计算
result = nRoot(25, 2) # R² + X² = Z²
print(result)
#应用: 电路阻抗计算
#数值结果: 5.0
示例8. 生物学应用
生物体积计算
result = nRoot(125, 3) # 细胞体积
print(result)
#应用: 细胞体积与尺寸关系
#数值结果: 5.0
种群增长
result = nRoot(8, 3) # 三代后的种群大小
print(result)
#应用: 种群几何增长模型
#数值结果: 2.0
示例9. 符号计算示例
symbolic_examples = [
(x**2, 2), # 平方根
(x**6, 3), # 立方根
(x**3 * y**3, 3), # 乘积的立方根
(x**4, 4), # 四次方根
((x+1)**3, 3), # 表达式的立方根
]
descriptions = [
"变量的平方根化简",
"高次幂的方根化简",
"乘积的方根分配律",
"四次方根计算",
"表达式的方根计算"
]
for example, desc in zip(symbolic_examples, descriptions):
result = nRoot(example)
print(f"{desc}: {example} = {result}")
#变量的平方根化简: (x**2, 2) = (x**2)**0.5
#高次幂的方根化简: (x**6, 3) = (x**6)**0.333333333333333
#乘积的方根分配律: (x**3 * y**3, 3) = (x**3*y**3)**0.333333333333333
#四次方根计算: (x**4, 4) = (x**4)**0.25
#表达式的方根计算: ((x+1)**3, 3) = ((x + 1.0)**3)**0.333333333333333
示例10. 复数计算
complex_examples = [
(-1, 2), # 虚数单位
(-8, 3), # 复数的立方根
(I, 4), # 虚数单位的四次方根
]
for example in complex_examples:
result = nRoot(example)
print(f"{example} = {result}")
#(-1, 2) = 1.0*I
#(-8, 3) = 1.0 + 1.73205080756888*I
#(I, 4) = 0.923879532511287 + 0.38268343236509*I
示例11. 实际工程问题
engineering_problems = [
# 建筑: 混凝土柱尺寸设计
((27, 3), "建筑: 27m³混凝土柱的边长"),
# 机械: 轴径设计
((8, 3), "机械: 扭矩与轴径的立方关系"),
# 电气: 电缆尺寸
((16, 2), "电气: 电流与导线截面积关系"),
# 水利: 管道流速
((64, 3), "水利: 流量与管径关系"),
]
for problem, desc in engineering_problems:
result = nRoot(problem)
print(f"{desc}: {problem} = {result}")
#建筑: 27m³混凝土柱的边长: (27, 3) = 3.0
#机械: 扭矩与轴径的立方关系: (8, 3) = 2.0
#电气: 电流与导线截面积关系: (16, 2) = 4.0
#水利: 流量与管径关系: (64, 3) = 4.0
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def nth_root_symbolic(input_str):
"""
对输入的字符串表达式进行符号化处理，计算其n次方根。
参数:
input_str: 输入的字符串，格式应为包含两个元素的元组，如 "(z, n)"，其中z是底数，n是根指数
返回:
如果输入格式正确，返回计算得到的n次方根的符号化结果；否则返回错误信息
"""
try:
# 将输入字符串转换为SymPy表达式
expr = sp.sympify(input_str)
error = False
result = None
# 检查输入是否为包含两个元素的元组
if isinstance(expr, tuple) and len(expr) == 2:
# 计算z的n次方根，并将结果进行数值计算
result = sp.root(*expr).evalf()
else:
# 若输入格式不符合要求，标记错误
error = True
# 如果没有错误，返回计算结果；否则返回错误信息
return result if not error else f"输入错误: {input_str}"
except Exception as e:
# 捕获并返回计算过程中出现的错误信息
return f"错误: {e}"
# 示例输入，计算4的2次方根
input_example1 = "(4, 2)"
result1 = nth_root_symbolic(input_example1)
print(f"输入: {input_example1}, 结果: {result1}")
# 2.00000000000000
# 示例输入，使用变量进行符号化计算
x = sp.Symbol('x')
input_example3 = f"({x}, 3)"
result3 = nth_root_symbolic(input_example3)
print(f"输入: {input_example3}, 结果: {result3}")
# x**0.333333333333333
非均匀快速傅里叶变换
Y = nufft(X,t) 使用采样点 t 返回 X 的非均匀离散傅里叶变换 (NUDFT)。
如果 X 是向量，则 nufft 返回该向量的变换。
如果 X 是矩阵，则 nufft 将 X 的各列视为向量，并返回每列的变换。
如果 X 是一个多维数组，则 nufft 将沿大小不等于 1 的第一个数组维度的值视为向量，并返回每个向量的变换。
Y = nufft(X,t,f) 使用采样点 t 计算查询点 f 的 NUDFT。要指定 f 而不指定采样点，请使用 nufft(X,[],f)。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from pynufft import NUFFT
def nufft_response(input_str):
"""
对标MATLAB的nufft函数实现非均匀快速傅里叶变换。
参数:
input_str: 输入表达式，格式为元组，其中：
- 第一个元素为数据（多维列表、SymPy矩阵、numpy数组）
- 第二个元素为非均匀采样坐标（归一化到[0, 1)范围）
- 第三个元素（可选）为输出尺寸s（默认为数据形状）
- 第四个元素（可选）为是否执行正向变换（默认为True）
返回:
- numpy数组或错误信息字符串
示例:
>>> data = np.array([1.0, 2.0, 3.0])
>>> t = np.array([[0.1], [0.5], [0.9]]) # 注意二维坐标输入
>>> nufft_response((data, t, (3,)))
array([...])
"""
try:
expr = sp.sympify(input_str, evaluate=False)
if not isinstance(expr, tuple) or len(expr) < 2:
return "错误：输入必须为元组，至少包含数据和坐标"
# 解析输入参数
data_part = expr[0]
t_part = expr[1]
s_param = expr[2] if len(expr) >= 3 else None
forward = expr[3] if len(expr) >= 4 else True
# 转换数据
# 转换数据(复数)和坐标(浮点)
data = np.array(data_part, dtype=complex) if isinstance(data_part, list) else None
t = np.array(t_part, dtype=np.float64) if isinstance(t_part, list) else None
if data is None:
return "错误：数据格式不支持"
# 转换坐标并确保二维结构
if t is None:
return "错误：坐标格式不支持"
# 确保坐标是二维数组 (N, d)
if t.ndim == 1:
t = t.reshape(-1, 1) # 转换为(N, 1)
# 归一化坐标到[0,1)并转换到[-π, π)
om = 2 * np.pi * t - np.pi
# 确定输出尺寸
if s_param is None:
s_param = data.shape
else:
# 确保s_param是元组格式
if isinstance(s_param, (list, tuple)):
s_param = tuple(s_param)
else:
s_param = (int(s_param),) # 标量转换为单元素元组
# 初始化NUFFT参数
ndim = om.shape[1] # 现在om保证是二维数组
Kd = tuple(2 * s for s in s_param) # 默认网格尺寸为2倍输出尺寸
Jd = (6,) * ndim # 邻近点数
# 创建NUFFT对象
nufft = NUFFT()
nufft.plan(om, s_param, Kd, Jd)
# 执行变换
if forward:
result = nufft.forward(data.flatten())
else:
result = nufft.adjoint(data.flatten())
return result.reshape(s_param)
except Exception as e:
return f"错误：{str(e)}"
# 正确使用的测试用例
data = [1.0, 2.0, 3.0]
t = [[0.1], [0.5], [0.9]] # 二维坐标输入
result = nufft_response(f"{data}, {t}, 3")
print("NUFFT结果:", result)
# [0.73477517+1.53991759j 5.98673999+0.00310206j 0.73477517-1.53991759j]
矩阵的零空间
Z = null(A) 返回A的零空间的标准正交基.
Z = null(A,tol) 也指定容差. 小于或等于tol的A的奇异值被视为零, 这会影响Z中的列数.
Z = null(A,"rational") 返回 A 的零空间的有理基, 它通常不是正交基. 如果A是具有小整数元素的小矩阵,则Z的元素是小整数的比率.此方法在数值上不如 null(A) 准确.
A — 输入矩阵,矩阵
tol — 奇异值容差,标量
Z — 零空间基向量,矩阵
示例1: 线性方程组求解 - 齐次方程组的解空间
方程组: x + 2y + 3z = 0, 2x + 4y + 6z = 0
A1 = [[1, 2, 3], [2, 4, 6]]
result1 = null(A1, rational)
print(f"系数矩阵 A = {A1}")
print(f"零空间基向量: {result1}")
#系数矩阵 A = [[1, 2, 3], [2, 4, 6]]
#零空间基向量: [[-2],
[1],
[0],
[-3],
[0],
[1]]
#物理意义: 描述所有满足 Ax=0 的解向量
示例2: 结构力学 - 机构的自由度分析
简化的桁架结构约束矩阵
A2 = [[1, 0, -1, 0],
[0, 1, 0, -1],
[1, 1, 0, 0]]
result2 = null(A2, rational)
print(f"结构约束矩阵 A = {A2}")
print(f"机构自由度 (刚体运动模式): {result2}")
#结构约束矩阵 A = [[1, 0, -1, 0], [0, 1, 0, -1], [1, 1, 0, 0]]
#机构自由度 (刚体运动模式): [[-1],
[1],
[-1],
[1]]
#物理意义: 零空间向量表示可能的刚体位移模式
示例3: 电路理论 - 基尔霍夫电流定律
节点关联矩阵 (4个节点，5条支路)
A3 = [[1, 1, 0, 0, 0], # 节点1
[-1, 0, 1, 1, 0], # 节点2
[0, -1, -1, 0, 1], # 节点3
[0, 0, 0, -1, -1]] # 节点4
result3 = null(A3, rational)
print(f"电路节点关联矩阵 A = {A3}")
print(f"回路电流模式: {result3}")
#电路节点关联矩阵 A = [[1, 1, 0, 0, 0], [-1, 0, 1, 1, 0], [0, -1, -1, 0, 1], [0, 0, 0, -1, -1]]
#回路电流模式: [[1],
[-1],
[1],
[0],
[0],
[-1],
[1],
[0],
[-1],
[1]]
#物理意义: 零空间表示满足KCL的独立回路电流
示例4: 控制理论 - 系统的能控性分析
能控性矩阵的零空间
A4 = [[1, 0.5, 0.25],
[0, 1, 1],
[0, 0, 1]]
B4 = [[1], [0], [1]]
能控性矩阵: [B, AB, A²B]
C_ctrl = np.hstack([B4, np.dot(A4, B4), np.dot(np.dot(A4, A4), B4)])
result4 = null(C_ctrl, 1e-10)
print(f"能控性矩阵维度: {C_ctrl.shape}")
print(f"不能控子空间: {result4}")
#能控性矩阵维度: (3, 3)
#不能控子空间: Matrix(3, 0, [])
#物理意义: 零空间表示系统不能控的状态方向
示例5: 机器人学 - 机械臂雅可比矩阵的奇异性
简化2连杆机械臂的雅可比矩阵
J5 = [
[-sp.sin(theta1) - sp.sin(theta1 + theta2), -sp.sin(theta1 + theta2)],
[sp.cos(theta1) + sp.cos(theta1 + theta2), sp.cos(theta1 + theta2)]
]
在特定构型下求值
J5_num = J5.subs({theta1: 0, theta2: 0})
result5 = null(J5_num, rational)
print(f"机械臂雅可比矩阵在(0,0)位置: {J5_num}")
print(f"奇异性方向 (零空间): {result5}")
#机械臂雅可比矩阵在(0,0)位置: [[0, 0], [2, 1]]
#奇异性方向 (零空间): [[-1/2],
[1]]
#物理意义: 零空间表示末端无法运动的方向
示例6: 计算机图形学 - 投影变换的核
透视投影矩阵 (简化版本)
A6 = [[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 1],
[0, 0, 0, 0]]
result6 = null(A6, rational)
print(f"投影矩阵 A = {A6}")
print(f"投影消失点 (零空间): {result6}")
#投影矩阵 A = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 1], [0, 0, 0, 0]]
#投影消失点 (零空间): [[0],
[0],
[-1],
[1]]
#物理意义: 零空间向量映射到投影平面的消失点
示例7: 统计学 - 线性回归的多重共线性
存在多重共线性的设计矩阵
A7 = [[1, 2, 3, 5],
[1, 4, 6, 10],
[1, 1, 1.5, 2.5],
[1, 3, 4.5, 7.5]]
result7 = null(A7, 1e-10)
print(f"设计矩阵 A = {A7}")
print(f"共线性关系: {result7}")
#设计矩阵 A = [[1, 2, 3, 5], [1, 4, 6, 10], [1, 1, 1.5, 2.5], [1, 3, 4.5, 7.5]]
#共线性关系: [[0, 0],
[-0.550220614995865, -0.769411539385026],
[0.794593530209707, -0.363013796550100],
[-0.256667872127478, 0.525572893684070]]
#物理意义: 零空间揭示预测变量之间的线性依赖关系
示例8: 量子力学 - 哈密顿算符的零能量态
简单量子系统的哈密顿矩阵
A8 = [[1, -1, 0],
[-1, 2, -1],
[0, -1, 1]]
result8 = null(A8, 1e-10)
print(f"哈密顿矩阵 H = {A8}")
print(f"零能量量子态: {result8}")
#哈密顿矩阵 H = [[1, -1, 0], [-1, 2, -1], [0, -1, 1]]
#零能量量子态: [[-0.577350269189626],
[-0.577350269189626],
[-0.577350269189626]]
#物理意义: 零空间对应系统的零能量本征态
示例9: 经济学 - 列昂惕夫投入产出模型
投入产出系数矩阵 (3个部门)
A9 = [[0.2, 0.3, 0.1],
[0.4, 0.1, 0.2],
[0.1, 0.2, 0.3]]
I - A 矩阵
I_minus_A = np.eye(3) - np.array(A9)
result9 = null(I_minus_A, 1e-10)
print(f"投入产出矩阵 (I-A) = {I_minus_A}")
print(f"平衡增长路径: {result9}")
#投入产出矩阵 (I-A) =
#[[ 0.8 -0.3 -0.1]
[-0.4 0.9 -0.2]
[-0.1 -0.2 0.7]]
#平衡增长路径: Matrix(3, 0, [])
#物理意义: 零空间表示经济系统的平衡增长模式
示例10: 网络理论 - 图的拉普拉斯矩阵
简单图的拉普拉斯矩阵
A10 = [[2, -1, -1, 0],
[-1, 2, -1, 0],
[-1, -1, 3, -1],
[0, 0, -1, 1]]
result10 = null(A10, rational)
print(f"图拉普拉斯矩阵 L = {A10}")
print(f"连通分量指示向量: {result10}")
#图拉普拉斯矩阵 L = [[2, -1, -1, 0], [-1, 2, -1, 0], [-1, -1, 3, -1], [0, 0, -1, 1]]
#连通分量指示向量: [[1],
[1],
[1],
[1]]
#物理意义: 零空间维度等于图的连通分量个数
示例11: 结构静力学 - 超静定桁架的静不定力
3杆桁架，1个多余约束
S = [[1, 1, 0],
[0, 1, 1],
[1, 0, 1]]
result = null(S, rational)
print(f"平衡矩阵 S = {S}")
print(f"自应力模式 (静不定力): {result}")
#平衡矩阵 S = [[1, 1, 0], [0, 1, 1], [1, 0, 1]]
#自应力模式 (静不定力): Matrix(0, 0, [])
#工程意义: 零空间表示结构中的自平衡内力分布
示例12: 电力系统 - 潮流方程的零注入节点
简化的节点导纳矩阵
Y = [[2, -1, -1],
[-1, 2, -1],
[-1, -1, 2]]
result = null(Y, 1e-10)
print(f"节点导纳矩阵 Y = {Y}")
print(f"零注入模式: {result}")
#节点导纳矩阵 Y = [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]]
#零注入模式: [[0.577350269189626],
[0.577350269189626],
[0.577350269189626]]
#工程意义: 零空间对应系统中功率注入为零的运行状态
示例13: 化工过程 - 反应系统的原子平衡
化学反应计量矩阵 (CO + H2O -> CO2 + H2)
物种: CO, H2O, CO2, H2
元素: C, O, H
A = [[1, 0, 1, 0], # C原子平衡
[1, 1, 2, 0], # O原子平衡
[0, 2, 0, 2]] # H原子平衡
result = null(A, rational)
print(f"原子矩阵 A = {A}")
print(f"独立的化学反应: {result}")
#原子矩阵 A = [[1, 0, 1, 0], [1, 1, 2, 0], [0, 2, 0, 2]]
#独立的化学反应: [[-1],
[-1],
[1],
[1]]
#工程意义: 零空间给出系统可能的独立反应计量式
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def null_space_matrix(input_str):
"""
MATLAB null函数的Python实现，支持符号计算和数值计算
参数:
input_str: 输入字符串，支持格式：
- 纯矩阵（如"[[1,2],[3,6]]"）
- 矩阵+计算方式（如"(Matrix([[1,2],[3,6]]), 'rational')"）
- 矩阵+数值容差（如"(Matrix([[1,2],[3,6]]), 1e-5)"）
返回:
SymPy矩阵 或 错误信息字符串
示例:
>>> null_space_matrix("([[1,2],[3,6]], 'rational')")
Matrix([[-2], [1]])
>>> null_space_matrix("([[1,2],[3,6]], 1e-5)")
Matrix([[0.89442719], [-0.4472136]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def null_space(A, rtol=1e-5):
"""
使用SVD的数值方法计算零空间
参数:
A: numpy数组格式的矩阵
rtol: 奇异值判断阈值
返回:
零空间基向量组成的numpy数组（列向量形式）
"""
u, s, vh = np.linalg.svd(A)
# 识别接近零的奇异值（根据机器精度调整阈值）
null_mask = (s <= rtol)
# 提取对应的右奇异向量并转置为列向量
return np.compress(null_mask, vh, axis=0).T
# 处理带参数的输入，格式：(矩阵, 计算方式)
if isinstance(expr, tuple) and len(expr) == 2:
matrix_part = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
if matrix_part is None:
return f"输入错误: 无法解析矩阵部分 {input_str}"
# 符号计算模式
if str(expr[1]) == "rational":
result = matrix_part.nullspace()
# 数值计算模式
elif expr[1].is_number:
A = np.array(matrix_part.tolist(), dtype=float)
result = null_space(A, rtol=float(expr[1]))
# 处理无参数的纯矩阵输入
elif isinstance(expr, list):
A = np.array(expr, dtype=float)
result = null_space(A)
else:
error = True
# 结果格式处理
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"计算错误: {str(e)}"
# 示例测试代码
if __name__ == "__main__":
# 测试1：符号计算
print("符号计算测试:")
test_input1 = "([[1,2],[3,6]], 'rational')"
print(f"输入: {test_input1}")
print("输出:", null_space_matrix(test_input1))
# Matrix([[-2], [1]])
# 测试2：数值计算
print("\n数值计算测试:")
test_input2 = "([[1,2],[3,6]], 1e-5)"
print(f"输入: {test_input2}")
print("输出:", null_space_matrix(test_input2))
# Matrix([[-0.894427190999916], [0.447213595499958]])
提取分子和分母
[N，D]=numden(A)将A转换为有理形式,其中分子和分母是具有整数系数的相对素数多项式.
函数返回表达式有理形式的分子和分母.
如果A是符号矩阵或数字矩阵,那么N是分子的符号矩阵,D是分母的符号矩阵.N和D都是与A大小相同的矩阵.
示例1: 控制系统 - 传递函数分析
典型二阶系统传递函数
num1, den1 = numden((omega**2) / (s**2 + 2*alpha*omega*s + omega**2))
print(f"传递函数: G(s) = {G1}")
print(f"分子 (零点多项式): {num1}")
print(f"分母 (极点多项式): {den1}")
#传递函数: G(s) = (omega**2) / (s**2 + 2*alpha*omega*s + omega**2)
#分子 (零点多项式): omega**2
#分母 (极点多项式): 2*alpha*omega*s + omega**2 + s**2
#应用: 分析系统稳定性、频率响应特性
示例2: 电路分析 - 阻抗函数分解
RLC串联电路阻抗
Z2 = R + s*L + 1/(s*C)
Z2_rational = Z2/1 # 转换为有理式
num2, den2 = numden(Z2_rational)
print(f"阻抗函数: Z(s) = {Z2}")
print(f"分子: {num2}")
print(f"分母: {den2}")
#阻抗函数: Z(s) = R + s*L + 1/(s*C)
#分子: C*L*s**2 + C*R*s + 1
#分母: C*s
#应用: 确定谐振频率、分析频率特性
示例3: 信号处理 - 滤波器传递函数
IIR滤波器传递函数
num3, den3 = numden((b0 + b1*z**(-1) + b2*z**(-2)) / (1 + a1*z**(-1) + a2*z**(-2)))
print(f"滤波器传递函数: H(z) = {H3}")
print(f"分子系数 (前馈部分): {num3}")
print(f"分母系数 (反馈部分): {den3}")
#滤波器传递函数: H(z) = (b0 + b1*z**(-1) + b2*z**(-2)) / (1 + a1*z**(-1) + a2*z**(-2))
#分子系数 (前馈部分): b0*z**2 + b1*z + b2
#分母系数 (反馈部分): a1*z + a2 + z**2
#应用: 滤波器实现、稳定性分析
示例4: 机械振动 - 频响函数
单自由度系统频响函数
num4, den4 = numden(1 / (m*s**2 + c*s + k))
print(f"频响函数: H(s) = {H4}")
print(f"分子: {num4}")
print(f"分母: {den4}")
#频响函数: H(s) = 1 / (m*s**2 + c*s + k)
#分子: 1
#分母: c*s + k + m*s**2
#应用: 模态分析、共振频率计算
示例5: 热力学 - 传热方程
一维热传导传递函数
num5, den5 = numden(1 / (tau*s + 1)) # 一阶系统
print(f"热传递函数: G(s) = {G5}")
print(f"分子: {num5}")
print(f"分母: {den5}")
#热传递函数: G(s) = 1 / (tau*s + 1)
#分子: 1
#分母: s*tau + 1
#应用: 热惯性分析、温度控制设计
示例6: 电磁学 - 波阻抗计算
介质波阻抗
num9, den9 = numden(sqrt(mu/epsilon))
print(f"波阻抗: Z = {Z9}")
print(f"分子: {num9}")
print(f"分母: {den9}")
#波阻抗: Z = sqrt(mu/epsilon)
#分子: sqrt(mu/epsilon)
#分母: 1
#应用: 电磁波传播分析、阻抗匹配
示例7: 流体力学 - 雷诺数表达式
雷诺数
num10, den10 = numden((rho*v*L) / mu)
print(f"雷诺数: Re = {Re}")
print(f"分子 (惯性力): {num10}")
print(f"分母 (粘性力): {den10}")
#雷诺数: Re = (rho*v*L) / mu
#分子 (惯性力): L*rho*v
#分母 (粘性力): mu
#应用: 流动状态判断、相似理论
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy import Matrix, sympify
def numerator_denominator(input_str):
"""
MATLAB numden函数的Python实现，支持矩阵和符号表达式
参数:
input_str: 输入字符串，可以是：
- 符号表达式（如"x/y"）
- 矩阵（如"[[1/x, 2], [3, 4/y]]"）
返回:
(分子, 分母) 元组，类型可能是：
- 对于符号表达式：两个SymPy表达式
- 对于矩阵：两个SymPy矩阵
- 错误信息字符串
示例:
>>> numerator_denominator("(x**2 + 1)/y")
(x**2 + 1, y)
>>> numerator_denominator("[[1/x, 2], [3, 4/y]]")
(Matrix([
[1, 2],
[3, 4]]), Matrix([
[x, 1],
[1, y]]))
"""
try:
expr = sympify(input_str)
# 处理矩阵情况
if isinstance(expr, list):
sp_matrix = sp.Matrix(expr)
# 初始化分子分母矩阵
num_matrix = []
den_matrix = []
# 遍历矩阵每个元素
for element in sp_matrix:
num, den = sp.fraction(element)
num_matrix.append(num)
den_matrix.append(den)
# 重构矩阵保持原始形状
original_shape = sp_matrix.shape
return (
Matrix(num_matrix).reshape(*original_shape),
Matrix(den_matrix).reshape(*original_shape)
)
else:
# 对普通表达式直接分解
return sp.fraction(sp.cancel(expr))
except Exception as e:
return f"错误: {str(e)}"
# 示例测试
if __name__ == "__main__":
# 测试符号表达式
print("符号表达式测试:")
test1 = "(x**2 + 1)/y"
print(f"输入: {test1}")
num, den = numerator_denominator(test1)
print(f"分子: {num}")
# x**2 + 1
print(f"分母: {den}\n")
# y
# 测试矩阵输入
print("矩阵输入测试:")
test2 = "[[1/x, 2], [3, 4/y]]"
print(f"输入: {test2}")
m_num, m_den = numerator_denominator(test2)
print("分子矩阵:")
print(m_num)
# Matrix([[1, 2],
# [3, 4]])
print("分母矩阵:")
print(m_den)
# Matrix([[x, 1],
# [1, y]])
纳托尔的布莱克曼-哈里斯窗
w = nuttallwin(L,sym=1) 返回一个长度为L个点的修正的纳托尔的布莱克曼-哈里斯窗 默认sym=1
当sym=1, 生成一个对称窗口，用于滤波器设计.
当sym=0, 生成一个周期性窗口，用于光谱分析.
示例1: 频谱分析 - 高动态范围信号的精确测量
生成测试信号：强信号与弱信号共存
fs = 1000 # 采样率
t = np.linspace(0, 1, 1000)
signal_strong = 1.0 * np.sin(2 * np.pi * 50 * t) # 强信号
signal_weak = 0.001 * np.sin(2 * np.pi * 150 * t) # 弱信号（动态范围60dB）
test_signal = signal_strong + signal_weak
应用Nuttall窗
window_1 = nuttallwin(1000)
windowed_signal = test_signal * window_1
频谱分析
fft_result = np.fft.fft(windowed_signal)
freq = np.fft.fftfreq(len(fft_result), 1 / fs)
props = analyze_window_properties(window_1, 1000)
print(f"窗口长度: 1000")
print(f"主瓣宽度: {props['main_lobe_width']:.4f} Hz")
print(f"旁瓣衰减: {props['sidelobe_attenuation']:.1f} dB")
#窗口长度: 1000
#主瓣宽度: 0.0010 Hz
#旁瓣衰减: 0.0 dB
#应用: 检测强信号附近的弱信号成分
示例2: 雷达信号处理 - 脉冲压缩中的旁瓣抑制
线性调频信号（LFM）
pulse_width = 10e-6 # 脉冲宽度
bandwidth = 1e6 # 带宽
t_radar = np.linspace(-pulse_width / 2, pulse_width / 2, 256)
chirp_signal = np.exp(1j * np.pi * bandwidth / pulse_width * t_radar ** 2)
应用Nuttall窗进行加窗处理
window_2 = nuttallwin(256)
windowed_chirp = chirp_signal * window_2
脉冲压缩（匹配滤波）
compressed = np.correlate(windowed_chirp, chirp_signal, mode='same')
props = analyze_window_properties(window_2, 256)
print(f"窗口长度: 256")
print(f"旁瓣衰减: {props['sidelobe_attenuation']:.1f} dB")
#窗口长度: 256
#旁瓣衰减: 0.0 dB
#应用: 提高距离分辨率，减少虚假目标
示例3: 音频处理 - 音乐信号的谐波分析
模拟音乐信号（基频+谐波）
t_audio = np.linspace(0, 0.1, 4410) # 100ms @ 44.1kHz
fundamental = 440 # A4音
audio_signal = (0.8 * np.sin(2 * np.pi * fundamental * t_audio) +
0.5 * np.sin(2 * np.pi * 2 * fundamental * t_audio) +
0.3 * np.sin(2 * np.pi * 3 * fundamental * t_audio) +
0.1 * np.sin(2 * np.pi * 4 * fundamental * t_audio))
window_3 = nuttallwin(4410)
windowed_audio = audio_signal * window_3
props = analyze_window_properties(window_3, 4410)
print(f"窗口长度: 4410 (100ms @ 44.1kHz)")
print(f"主瓣宽度: {props['main_lobe_width'] * 44100:.1f} Hz")
#窗口长度: 4410 (100ms @ 44.1kHz)
#主瓣宽度: 10.8 Hz
#应用: 精确分离音乐信号的谐波成分
示例4: 通信系统 - OFDM子载波的频谱整形
OFDM符号（多个子载波）
n_subcarriers = 64
ofdm_symbol = np.sum([np.exp(1j * 2 * np.pi * k * np.arange(256) / 256)
for k in range(n_subcarriers)], axis=0)
window_4 = nuttallwin(256)
windowed_ofdm = ofdm_symbol * window_4
props = analyze_window_properties(window_4, 256)
print(f"窗口长度: 256")
print(f"频谱泄漏抑制: {props['sidelobe_attenuation']:.1f} dB")
#窗口长度: 256
#频谱泄漏抑制: 0.0 dB
#应用: 减少OFDM系统子载波间干扰
示例5: 振动分析 - 轴承故障特征提取
模拟轴承振动信号（冲击响应+噪声）
t_vib = np.linspace(0, 1, 5000)
bearing_freq = 100 # 故障特征频率
impact_times = np.arange(0, 1, 1 / bearing_freq)
vibration = np.random.normal(0, 0.1, len(t_vib)) # 背景噪声
添加周期性冲击
for impact_time in impact_times:
idx = int(impact_time * len(t_vib))
if idx < len(t_vib):
vibration[idx:idx + 50] += np.exp(-np.arange(50) / 10) * np.sin(2 * np.pi * 1000 * np.arange(50) / 5000)
window_5 = nuttallwin(5000)
windowed_vibration = vibration * window_5
props = analyze_window_properties(window_5, 5000)
print(f"窗口长度: 5000")
print(f"等效噪声带宽: {props['equivalent_noise_bandwidth']:.3f}")
#窗口长度: 5000
#等效噪声带宽: 1.977
#应用: 提高故障特征频率的检测灵敏度
示例6: 医学影像 - B模式成像的旁瓣抑制
模拟超声回波信号
t_ultrasound = np.linspace(0, 50e-6, 1000)
center_freq = 3.5e6 # 3.5MHz中心频率
ultrasound_signal = np.exp(-(t_ultrasound - 25e-6) ** 2 / (5e-6) ** 2) * np.sin(
2 * np.pi * center_freq * t_ultrasound)
window_6 = nuttallwin(1000)
windowed_ultrasound = ultrasound_signal * window_6
props = analyze_window_properties(window_6, 1000)
print(f"窗口长度: 1000")
print(f"旁瓣衰减: {props['sidelobe_attenuation']:.1f} dB")
#窗口长度: 1000
#旁瓣衰减: 0.0 dB
#应用: 提高医学超声图像的空间分辨率
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import scipy.signal as signal
import sympy as sp
from typing import Union, List, Tuple
def analyze_window_properties(window: List[float], M: int) -> dict:
"""分析窗口函数的特性"""
n = np.arange(M)
w = np.array(window)
# 计算频域特性
W = np.fft.fft(w, 4096)
freq = np.fft.fftfreq(4096)
mag = 20 * np.log10(np.abs(W) + 1e-10)
# 找到主瓣宽度和旁瓣衰减
main_lobe_width = None
sidelobe_attenuation = None
if len(mag) > 0:
peak_idx = np.argmax(mag)
# 主瓣宽度（-3dB点）
half_power = mag[peak_idx] - 3
left_idx = peak_idx
while left_idx > 0 and mag[left_idx] > half_power:
left_idx -= 1
right_idx = peak_idx
while right_idx < len(mag) - 1 and mag[right_idx] > half_power:
right_idx += 1
main_lobe_width = (right_idx - left_idx) * (freq[1] - freq[0])
# 旁瓣衰减
sidelobe_region = np.concatenate([mag[:peak_idx - 10], mag[peak_idx + 10:]])
sidelobe_attenuation = mag[peak_idx] - np.max(sidelobe_region) if len(sidelobe_region) > 0 else 0
return {
'main_lobe_width': main_lobe_width,
'sidelobe_attenuation': sidelobe_attenuation,
'coherent_gain': np.mean(w),
'processing_gain': np.sum(w ** 2) / M,
'equivalent_noise_bandwidth': M * np.sum(w ** 2) / (np.sum(w) ** 2)
}
def nuttallwin_window(input_str):
"""
MATLAB nuttallwin函数的Python实现
参数:
params: 输入参数，支持以下格式：
- 单个整数M（窗口长度）
- 元组(M, sym_flag)：
M: 窗口长度（正整数）
sym_flag: 对称性标志（0表示周期性窗口，1表示对称窗口）
返回:
list: 生成的Nuttall窗口系数列表
或错误信息字符串
示例:
>>> nuttallwin_window(5)
[0.0, 0.2235, 0.8705, 0.8705, 0.2235]
>>> nuttallwin_window((6, 0))
[0.0003628, 0.061334, 0.52923, 1.0, 0.52923, 0.061334]
"""
try:
expr = sp.sympify(input_str)
error = False
sym = True
if isinstance(expr, tuple) and len(expr) == 2:
M = int(expr[0])
sym_flag = expr[1]
# 验证sym_flag有效性
if sym_flag not in (0, 1):
raise ValueError("sym_flag 必须为 0 或 1")
sym = bool(sym_flag)
elif expr.is_number:
# 处理单个整数输入
M = int(expr)
else:
error = True
# 验证M有效性
if not isinstance(M, int):
raise TypeError("窗口长度M必须为整数")
if M <= 0:
raise ValueError("窗口长度必须为正整数")
# 生成窗口
window = signal.windows.nuttall(M, sym=sym)
return list(window) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例测试
if __name__ == "__main__":
# 测试基本功能
print("基本测试:")
print("nuttallwin_window(5) =>", nuttallwin_window(5))
# [0.0003628000000000381, 0.22698240000000006, 1.0, 0.22698240000000006, 0.0003628000000000381]
奈奎斯特图
nyquistplot(sys1,sys2...sysN)控制系统中用于分析系统稳定性的频域工具，它在复平面上绘制开环传递函数在频率从负无穷到正无穷变化时的轨迹。
该图以美国工程师哈里·奈奎斯特（Harry Nyquist）命名，是控制系统分析与设计中的重要工具。
直流电机位置控制系统
分析电机位置控制系统的稳定性
确定系统在不同增益下的稳定裕度
评估系统对负载变化的鲁棒性
预测超调和振荡行为
nyquistplot(5,s*(s+1)*(s+5))
温度控制系统
评估工业加热炉的温度控制性能
分析系统对温度扰动的响应特性
确定防止温度超调的最佳控制器参数
预测系统在环境温度变化时的稳定性
nyquistplot(2(s+0.5), s*(s+1)*(s+2)*(s+0.1))
飞机俯仰控制系统
分析飞行控制系统在俯仰轴上的稳定性
评估飞机在湍流中的姿态保持能力
确定自动驾驶模式下控制参数的合理范围
预测系统在极端飞行条件下的稳定性边界
nyquistplot(0.5(s+0.2)*(s+2), s^2*(s+1)*(s+5))
机械振动控制系统
评估主动减振系统的性能
分析系统在共振频率附近的稳定性
确定抑制机械振动的有效控制策略
预测系统在负载变化时的振动抑制效果
nyquistplot(0.8*(s^2+0.2s+1), s*(s^2+0.1s+4))
电源稳压系统
分析开关电源的稳定性
评估输出电压在负载突变时的恢复能力
确定防止振荡的最佳补偿网络参数
预测系统在输入电压波动时的稳定性
nyquistplot(10*(s+100), s*(s+50)*(s+200))
稳定性分析（条件稳定系统）
观察轨迹是否包围(-1,0)点。当K=5时轨迹包围(-1,0)点两次（闭环不稳定），K=1时轨迹不包围（闭环稳定）。用于设计航天器姿态控制器的增益范围。
nyquistplot([K,2K],[1,2,-3,0])
电机控制系统
确保工业机器人关节电机在负载变化时保持稳定（增益裕度>6dB）
nyquistplot([10],[1,5.5,2.5,0])
非最小相位系统（化工过程控制）
分析精馏塔温度控制系统的反常行为（右半平面零点导致反向响应）。
nyquistplot([-2,2],[1,3,2])
鲁棒性分析（飞行控制系统）
保证战斗机在气动参数摄动时（如湍流中）仍能稳定飞行（PM>40°, GM>6dB）。
nyquistplot([20,2],[1,12,20,0])
多曲线奈奎斯特图
将两个或多个传递函数的奈奎斯特图绘制在同一坐标系中（例如 nyquistplot([4s²+5s+1, 1], [3s²+2s+5, s+1/3])），这种可视化方法具有重要的工程意义和分析价值。
不同积分时间的 PI 控制器
nyquistplot([s+5,s^2+2s],[s+3,s^2+2s],[s+1,s^2+2s])
不同阻尼比的二阶系统
nyquistplot([1,s^2+0.4s+1],[1,s^2+s+1],[1,s^2+1.6s+1],[1,s^2+2s+1])
标称系统与鲁棒控制器对比
nyquistplot([s+0.5, s^3+3s^2+2s+0.5],[2s^2+s+0.2, s^3+3.5s^2+2.2s+0.6])