幻方矩阵
M = magic(n) 返回由1到n^2的整数构成并且总行数和总列数相等的n×n矩阵.n的阶数必须是大于或等于3的标量才能创建有效的幻方矩阵.
n — 矩阵的阶次, 整数标量
示例1. 数学教育：魔方性质验证
n = 5
matrix = magic(n)
magic_matrix = np.array(matrix, dtype=int)
验证魔方性质
row_sums = np.sum(magic_matrix, axis=1)
col_sums = np.sum(magic_matrix, axis=0)
main_diag_sum = np.trace(magic_matrix)
anti_diag_sum = np.trace(np.fliplr(magic_matrix))
print(f"{n}阶魔方阵:")
for row in magic_matrix:
print(row)
#5阶魔方阵:
#[17 24 1 8 15]
#[23 5 7 14 16]
#[ 4 6 13 20 22]
#[10 12 19 21 3]
#[11 18 25 2 9]
print(f"\n行和: {row_sums}")
print(f"列和: {col_sums}")
print(f"主对角线和: {main_diag_sum}")
print(f"副对角线和: {anti_diag_sum}")
print(f"魔方常数: {n * (n ** 2 + 1) // 2}")
#行和: [65 65 65 65 65]
#列和: [65 65 65 65 65]
#主对角线和: 65
#副对角线和: 65
#魔方常数: 65
验证所有和是否相等
all_sums = np.concatenate([row_sums, col_sums, [main_diag_sum, anti_diag_sum]])
is_magic = np.all(all_sums == all_sums[0])
print(f"是否为有效魔方阵: {is_magic}")
#是否为有效魔方阵: True
示例2. 艺术设计：对称图案生成
n = 4
matrix = magic(n)
magic_matrix = np.array(matrix, dtype=int)
将魔方阵转换为灰度图案
normalized = (magic_matrix - magic_matrix.min()) / (magic_matrix.max() - magic_matrix.min())
print("魔方矩阵:")
for row in magic_matrix:
print(row)
#魔方矩阵:
#[16 2 3 13]
#[ 5 11 10 8]
#[ 9 7 6 12]
#[ 4 14 15 1]
归一化值（可用于像素强度）:
for row in normalized:
print([f"{val:.2f}" for val in row])
#['1.00', '0.07', '0.13', '0.80']
#['0.27', '0.67', '0.60', '0.47']
#['0.53', '0.40', '0.33', '0.73']
#['0.20', '0.87', '0.93', '0.00']
#实际意义：在数字艺术中，魔方阵可用于生成对称纹理
示例3. 密码学：简单置换密码
n = 4
matrix = magic(n)
magic_matrix = np.array(matrix, dtype=int)
创建置换规则：按魔方阵数值排序的索引
flat_indices = np.argsort(magic_matrix.flatten())
permutation = flat_indices.reshape(magic_matrix.shape)
print("魔方矩阵:")
for row in magic_matrix:
print(row)
#魔方矩阵:
#[16 2 3 13]
#[ 5 11 10 8]
#[ 9 7 6 12]
#[ 4 14 15 1]
print("\n置换规则矩阵:")
for row in permutation:
print(row)
#置换规则矩阵:
#[15 1 2 12]
#[ 4 10 9 7]
#[ 8 6 5 11]
#[ 3 13 14 0]
加密消息
message = "HELLO WORLD123456"
# 将消息填充到矩阵大小
padded_msg = message[:n * n].ljust(n * n, 'X')
msg_matrix = np.array(list(padded_msg)).reshape(n, n)
print(f"\n原始消息矩阵:")
for row in msg_matrix:
print(row)
#原始消息矩阵:
#['H' 'E' 'L' 'L']
#['O' ' ' 'W' 'O']
#['R' 'L' 'D' '1']
#['2' '3' '4' '5']
应用置换
encrypted = np.empty_like(msg_matrix)
for i in range(n):
for j in range(n):
encrypted[i, j] = msg_matrix[permutation[i, j] // n, permutation[i, j] % n]
print(f"\n加密后矩阵:")
for row in encrypted:
print(row)
#加密后矩阵:
#['5' 'E' 'L' '2']
#['O' 'D' 'L' 'O']
#['R' 'W' ' ' '1']
#['L' '3' '4' 'H']
encrypted_msg = ''.join(encrypted.flatten())
print(f"加密消息: {encrypted_msg}")
#加密消息: 5EL2ODLORW 1L34H
示例4. 游戏开发：魔方数独变体
n = 3
matrix = magic(n)
magic_matrix = np.array(matrix, dtype=int)
print("基础魔方阵:")
for row in magic_matrix:
print(row)
#基础魔方阵:
#[8 1 6]
#[3 5 7]
#[4 9 2]
创建数独风格的谜题（隐藏部分数字）
puzzle = magic_matrix.copy()
hide_mask = np.random.choice([True, False], size=(n, n), p=[0.5, 0.5])
puzzle[hide_mask] = 0
print(f"\n魔方数独谜题 (0表示待填数字):")
for row in puzzle:
print(row)
#魔方数独谜题 (0表示待填数字):
#[8 1 6]
#[0 5 7]
#[0 0 0]
print(f"\n游戏规则: 每行、每列、对角线之和必须等于{np.sum(magic_matrix[0])}")
#游戏规则: 每行、每列、对角线之和必须等于15
#实际意义：基于魔方阵设计数学谜题游戏
示例5. 统计学：拉丁方实验设计
n = 4
matrix = magic(n)
magic_matrix = np.array(matrix, dtype=int)
使用魔方阵设计平衡的实验条件
treatments = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P']
print("实验设计矩阵 (基于魔方阵):")
design_matrix = np.empty_like(magic_matrix, dtype=object)
for i in range(n):
for j in range(n):
design_matrix[i, j] = treatments[magic_matrix[i, j] % len(treatments)]
for row in design_matrix:
print(list(row))
#实验设计矩阵 (基于魔方阵):
#['A', 'C', 'D', 'N']
#['F', 'L', 'K', 'I']
#['J', 'H', 'G', 'M']
#['E', 'O', 'P', 'B']
验证设计的平衡性
treatment_counts = {}
for treatment in design_matrix.flatten():
treatment_counts[treatment] = treatment_counts.get(treatment, 0) + 1
print(f"\n处理出现次数: {treatment_counts}")
print("设计是否平衡:", len(set(treatment_counts.values())) == 1)
#处理出现次数: {'A': 1, 'C': 1, 'D': 1, 'N': 1, 'F': 1, 'L': 1, 'K': 1, 'I': 1, 'J': 1, 'H': 1, 'G': 1, 'M': 1, 'E': 1, 'O': 1, 'P': 1, 'B': 1}
#设计是否平衡: True
#实际意义：在实验设计中确保各处理条件均匀分布
示例6. 计算机图形学：魔方纹理坐标
n = 4
matrix = magic(n)
magic_matrix = np.array(matrix.tolist(), dtype=int)
生成基于魔方阵的纹理坐标
texture_coords = magic_matrix / (n * n) # 归一化到[0,1]
print("魔方矩阵:")
for row in magic_matrix:
print(row)
#魔方矩阵:
#[16 2 3 13]
#[ 5 11 10 8]
#[ 9 7 6 12]
#[ 4 14 15 1]
print("\n对应的纹理坐标:")
for row in texture_coords:
print([f"({i // n / n:.2f}, {i % n / n:.2f})" for i in row])
#对应的纹理坐标:
#['(0.00, 0.25)', '(0.00, 0.03)', '(0.00, 0.05)', '(0.00, 0.20)']
#['(0.00, 0.08)', '(0.00, 0.17)', '(0.00, 0.16)', '(0.00, 0.12)']
#['(0.00, 0.14)', '(0.00, 0.11)', '(0.00, 0.09)', '(0.00, 0.19)']
#['(0.00, 0.06)', '(0.00, 0.22)', '(0.00, 0.23)', '(0.00, 0.02)']
# 实际意义：在计算机图形学中生成非重复但结构化的纹理
示例7. 音乐理论：基于魔方阵的音阶设计
n = 3
matrix = magic(n)
magic_matrix = np.array(matrix.tolist(), dtype=int)
音阶音符
notes = ['C', 'D', 'E', 'F', 'G', 'A', 'B', 'C5', 'D5']
print("魔方矩阵:")
for row in magic_matrix:
print(row)
#魔方矩阵:
#[8 1 6]
#[3 5 7]
#[4 9 2]
将魔方阵映射到音阶
music_matrix = np.empty_like(magic_matrix, dtype=object)
for i in range(n):
for j in range(n):
music_matrix[i, j] = notes[magic_matrix[i, j] % len(notes)]
print(f"\n音阶排列矩阵:")
for row in music_matrix:
print(list(row))
#音阶排列矩阵:
#['D5', 'D', 'B']
#['F', 'A', 'C5']
#['G', 'C', 'E']
生成旋律序列
melody = []
for i in range(n):
for j in range(n):
melody.append(music_matrix[i, j])
print(f"旋律序列: {'-'.join(melody)}")
#旋律序列: D5-D-B-F-A-C5-G-C-E
#实际意义：在算法作曲中使用魔方阵生成结构化音乐
示例8. 数据加密：魔方阵混淆
n = 4
matrix = magic(n)
magic_matrix = np.array(matrix, dtype=int)
要加密的数据（16个字节）
data = list(range(1, 17)) # [1, 2, 3, ..., 16]
print("原始数据:", data)
print("魔方矩阵:")
for row in magic_matrix:
print(row)
#原始数据: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]
#魔方矩阵:
#[16 2 3 13]
#[ 5 11 10 8]
#[ 9 7 6 12]
#[ 4 14 15 1]
基于魔方阵重新排列数据
data_matrix = np.array(data).reshape(n, n)
encrypted_data = np.empty_like(data_matrix)
使用魔方阵的值作为索引映射
for i in range(n):
for j in range(n):
# 将魔方阵值映射到有效索引
idx = (magic_matrix[i, j] - 1) % (n * n)
encrypted_data[i, j] = data[idx]
print(f"\n加密后数据矩阵:")
for row in encrypted_data:
print(row)
#加密后数据矩阵:
#[16 2 3 13]
#[ 5 11 10 8]
#[ 9 7 6 12]
#[ 4 14 15 1]
encrypted_flat = encrypted_data.flatten()
print(f"加密后数据: {encrypted_flat}")
#加密后数据: [16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1]
#实际意义：简单的数据混淆技术，可用于低安全要求的场景
示例9. 建筑学：基于魔方阵的对称布局
n = 5
matrix = magic(n)
magic_matrix = np.array(matrix, dtype=int)
建筑元素类型
elements = ['柱', '墙', '窗', '门', '装饰', '通道', '庭院', '楼梯', '屋顶']
print("魔方矩阵:")
for row in magic_matrix:
print(row)
#魔方矩阵:
#[17 24 1 8 15]
#[23 5 7 14 16]
#[ 4 6 13 20 22]
#[10 12 19 21 3]
#[11 18 25 2 9]
将魔方阵映射到建筑布局
layout_matrix = np.empty_like(magic_matrix, dtype=object)
for i in range(n):
for j in range(n):
layout_matrix[i, j] = elements[magic_matrix[i, j] % len(elements)]
print(f"\n建筑布局设计:")
for row in layout_matrix:
print(list(row))
#建筑布局设计:
#['屋顶', '庭院', '墙', '屋顶', '庭院']
#['通道', '通道', '楼梯', '通道', '楼梯']
#['装饰', '庭院', '装饰', '窗', '装饰']
#['墙', '门', '墙', '门', '门']
#['窗', '柱', '楼梯', '窗', '柱']
#实际意义：在建筑设计中创建对称且平衡的空间布局
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def magic_square(input_str):
try:
expr = sp.sympify(input_str)
n = int(expr)
error = False
if n < 3:
error = True
return f"输入错误: 阶数必须≥3，当前输入{n}"
def odd_magic_square(n):
"""生成奇数阶魔方阵，使用Siamese方法"""
magic = np.zeros((n, n), dtype=int)
i, j = 0, n // 2 # 初始位置在第一行中间
for num in range(1, n ** 2 + 1):
magic[i, j] = num
if num % n == 0: # 当前数是n的倍数，下移一行
i = (i + 1) % n
else: # 否则，右上移动并处理边界
i = (i - 1) % n
j = (j + 1) % n
return magic
def doubly_even_magic_square(n):
"""生成双偶数阶（4k）魔方阵，替换4x4子块的对角线元素"""
arr = np.arange(1, n * n + 1).reshape(n, n)
i, j = np.indices((n, n))
# 每个4x4子块中的对角线元素（主对角线和副对角线）
sub_i, sub_j = i % 4, j % 4
mask = ((sub_i == sub_j) | (sub_i + sub_j == 3))
arr[mask] = n * n + 1 - arr[mask]
return arr
def singly_even_magic_square(n):
"""生成单偶数阶（4k+2）魔方阵，采用LUX方法组合子块"""
half = n // 2
sub = odd_magic_square(half) # 生成奇数阶子块
# 组合四个子块
magic = np.block([
[sub, sub + 2 * half * half],
[sub + 3 * half * half, sub + half * half]
])
k = (n - 2) // 4
# 交换前k列的部分行
for i in range(half):
for j in range(k):
if i == k and j == 0:
continue # 跳过中间行第一列
magic[[i, i + half], j] = magic[[i + half, i], j]
# 交换后k列的部分行
for i in range(half):
for j in range(n - k, n):
magic[[i, i + half], j] = magic[[i + half, i], j]
# 调整中间部分
for i in range(k):
col = i + k
magic[[i + half, i], col] = magic[[i, i + half], col]
return magic
# 根据阶数选择生成方法
if n % 2 == 1:
result = odd_magic_square(n)
elif n % 4 == 0:
result = doubly_even_magic_square(n)
else:
result = singly_even_magic_square(n)
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
test_cases = [3, 5, 6]
for n in test_cases:
print(f"\n{n}阶魔方阵:")
matrix = magic_square(str(n))
if isinstance(matrix, sp.Matrix):
print(matrix)
else:
print(matrix)
#Matrix([[8, 1, 6],
# [3, 5, 7],
# [4, 9, 2]])
# Matrix([[17, 24, 1, 8, 15],
# [23, 5, 7, 14, 16],
# [4, 6, 13, 20, 22],
# [10, 12, 19, 21, 3],
# [11, 18, 25, 2, 9]])
# Matrix([[35, 28, 6, 26, 19, 15],
# [3, 5, 7, 21, 23, 16],
# [31, 9, 2, 22, 27, 11],
# [8, 1, 33, 17, 10, 24],
# [30, 32, 34, 12, 14, 25],
# [4, 36, 29, 13, 18, 20]])
幅频图
magnitudeplot(sys1) 关注增益随频率的变化. 直观地展示了系统对不同频率信号的“放大”或“衰减”特性.
sys1是控制系统的传递函数
低通滤波器（一阶系统）: RC电路噪声滤除。
在噪声滤除应用中，工程师只需关注幅频图即可确定截止频率和衰减程度。
幅频图可优化元件值（如R、C）以满足噪声抑制要求。
magnitudeplot(1,0.1s+1)
直流电机速度控制: 电机调速系统
幅频图可帮助确定系统的速度响应带宽和低频增益，这对于调速性能至关重要。
幅频图可帮助确定系统的速度响应带宽和低频增益，这对于调速性能至关重要。
magnitudeplot(5,0.5*s^2+s)
质量-弹簧-阻尼系统: 机械振动分析
在振动分析中，幅频图直接显示谐振频率和峰值增益，这对避免共振或设计减振器至关重要。
机械结构设计中，仅幅频图可用于评估共振风险和操作频率安全区。
magnitudeplot(1,s^2+0.6s+1)
超前补偿器（PD控制器）: 提高系统稳定裕度
幅频图清晰显示增益提升的频率范围（5-50 rad/s），这直接关联到超前补偿效果（如提高响应速度）。
用于调整补偿器参数以增加开环增益，从而减小稳态误差或扩展带宽。
magnitudeplot(2s+10,0.02s+1)
带通滤波器: 通信信号选频
在通信选频应用中，幅频图直接定义通带中心频率和带宽，这对信号提取（如调制解调）至关重要。
射频（RF）系统设计中，仅幅频图用于快速验证滤波器性能。
magnitudeplot(0.1s,(s+1)(0.01s+1))
低通滤波器（RC电路）
阻高频、通低频。幅频图显示高频段增益衰减（如音频系统滤除噪声）。
magnitudeplot([1],[0.001,1])
高通滤波器（CR电路）
阻低频、通高频（如去除传感器信号的直流偏移）。
magnitudeplot([0.002,0],[0.002,1])
带通滤波器（RLC谐振电路）
仅通过特定频段（如无线电接收目标频率信号）。
magnitudeplot([0.000637,0],[1,62.8,98596])
陷波滤波器（消除工频干扰）
在50Hz处深度衰减（如医疗设备中消除电源干扰）。
magnitudeplot([1, 0, 98596], [1, 31.4, 98596])
控制系统稳定性分析（二阶系统）
观察共振峰（$\zeta<1$时出现峰值），评估系统稳定性。
magnitudeplot([10000], [1, 60, 10000])
叠加幅频图
magnitudeplot(sys1,sys2...sysN) 能够高效地对比不同系统的频率特性或同一系统的不同设计参数，从而加速分析和决策。
同时显示多个系统的 带宽、截止频率、谐振峰值、衰减斜率 等参数，无需反复切换图表。
在控制器或滤波器设计中，叠加不同参数（如RC值、阻尼比）的曲线，快速找到最优解。
在教学中，叠加一阶/二阶系统的曲线，直观展示阶数对衰减斜率的影响
sys是控制系统的传递函数
滤波器设计迭代对比: 优化RC低通滤波器的电容值选择
对比不同电容值在关键噪声频段（如50Hz工频/10kHz开关噪声）的衰减效果
快速确定满足衰减要求的电容值（如1kHz处-20dB需选C=1μF）
magnitudeplot([1,0.1@e-6*s+1], # C=0.1μF (截止频率1.59kHz)
[1,@e-6*s+1], # C=1μF (截止频率159Hz)
[1,10@e-6*s+1]) # C=10μF (截止频率15.9Hz)
控制系统稳定性优化: 调整PD控制器参数提升电机系统稳定裕度
观察PD控制器对增益交越频率的影响（理想上移2-5倍）
验证高频段增益衰减是否加剧（需避免高频噪声放大）
magnitudeplot([5,0.5*s^2+s], # 原始系统
[s+5,0.5*s^2+s], # 添加PD: Kp=1, Td=0.2
[2.5s+5,0.5*s^2+s]) # 强化PD: Kp=1, Td=0.5
机械系统共振分析: 不同阻尼比对振动系统的影响
量化谐振峰值降低效果（ζ=0.6时峰值衰减＞10dB）
确定避免共振的工作频率范围（如>1.5ωₙ）
magnitudeplot([1,s^2+0.1*s+1], # ζ=0.05 (强共振)
[1,s^2+0.6*s+1], # ζ=0.3 (适度阻尼)
[1,s^2+1.2*s+1]) # ζ=0.6 (过阻尼)
通信滤波器组性能验证: 多通道带通滤波器频率隔离度分析
检查相邻通道隔离度（如15rad/s处需>30dB衰减）
验证通带平坦度和3dB带宽一致性
magnitudeplot([0.1*s,(s+1)*(0.01*s+1)], # 中心频率10rad/s
[0.2*s,(s+2)*(0.05*s+1)], # 中心频率20rad/s
[0.3*s,(s+5)*(0.02*s+1)]) # 中心频率50rad/s
电源设计EMI优化: 开关电源LC滤波器拓扑比较
对比衰减斜率（单级-40dB/dec vs 两级-80dB/dec）
评估开关频率（100kHz）处抑制能力（两级改善>20dB）
magnitudeplot([1,1@e-6*s^2+1@e-3*s+1], # 单级LC
[1,(1@e-6*s^2+1@e-3*s+1)^2]) # 两级级联
传感器信号调理对比: 不同抗混叠滤波器对测量精度影响
分析通带纹波（Bessel<0.1dB vs 巴特沃斯>3dB）
比较过渡带特性与采样频率关系（奈奎斯特频率50Hz）
magnitudeplot([1,0.01*s+1], # 一阶(采样率100Hz)
[1,(0.002*s)^2+0.55*0.002*s+1], # 二阶Bessel
[1,(0.001*s)^3+1]) # 三阶巴特沃斯
电机控制器频响验证: PWM载频谐波抑制效果
量化PWM载频（如20kHz）处增益衰减（需<-30dB）
观察增益提高对带宽影响（Kp=10时带宽扩展5倍）
magnitudeplot([5,0.5*s^2+s], # 开环
[5,0.5*s^2+s+5], # 闭环Kp=1
[5,0.5*s^2+s+50]) # 闭环Kp=10
修正Akima分段三次Hermite插值
yq = makima(x,y,xq) 使用采样点 x 处的值 y 执行修正 Akima 插值，以求出查询点 xq 处的插值 yq。
pp = makima(x,y) 返回一个分段多项式结构体以用于 ppval 和样条实用工具 unmkpp。
x — 样本点, 向量
y — 样本点处的函数值, 向量 | 矩阵
xq — 查询点, 标量 | 向量 | 矩阵
yq — 查询点位置的插值, 标量 | 向量 | 矩阵
pp — 分段多项式, 结构体
示例1：物理实验数据插值 - 弹簧伸长量与力的关系
弹簧实验：测量不同质量对应的伸长量
result_spring = makima([0.1, 0.2, 0.3, 0.4, 0.5], [2.1, 4.3, 6.2, 8.4, 10.1], [0.15, 0.25, 0.35, 0.45])
print("弹簧实验插值结果:", result_spring)
#弹簧实验插值结果: [3.2625 5.2484375 7.3 9.3203125]
#可以估算中间质量对应的伸长量
示例2：经济学 - 产品价格与销量关系
不同价格下的产品销量数据
result_sales = makima([10, 20, 30, 40, 50], [1000, 800, 600, 400, 200], [15, 25, 35, 45])
print("价格-销量关系插值:", result_sales)
#价格-销量关系插值: [900. 700. 500. 300.]
#帮助制定最优定价策略
示例3：医学 - 药物浓度随时间变化
药物在血液中的浓度监测数据（小时，mg/L）
result_medical = makima([0, 2, 4, 6, 8, 12], [0, 15.2, 12.8, 8.5, 5.1, 1.2], [1, 3, 5, 7, 10])
print("药物浓度插值:", result_medical)
#药物浓度插值: [11.88733154 14.33861441 10.59774808 6.68189108 2.74507979]
#估算未监测时间点的药物浓度
示例4：工程学 - 材料温度与膨胀系数
生成符号表达式，用于理论分析
expr_thermal = makima([0, 50, 100, 150, 200], [1.2, 1.8, 2.5, 3.3, 4.2])
print("温度-膨胀系数关系:")
print(expr_thermal)
#温度-膨胀系数关系:
#Piecewise((-2.66666666666666e-7*t**3 + 5.33333333333333e-5*t**2 + 0.01*t + 1.2, (t <= 50.0) & (t >= 0)),
(0.0133333333333333*t + 0.0166666666666666*(0.02*t - 1)**3 + 0.0166666666666667*(0.02*t - 1)**2 + 1.13333333333333, (t <= 100.0) & (t > 50.0)),
(0.015*t + 0.200000000000001*(0.01*t - 1)**2 + 1.0, (t <= 150.0) & (t > 100.0)),
(0.017*t - 9.36750677027476e-15*(0.00666666666666667*t - 1)**3 + 0.450000000000007*(0.00666666666666667*t - 1)**2 + 0.75, (t <= 200.0) & (t > 150.0)))
#获得连续的函数表达式，便于后续计算
示例5：气象学 - 多点温度监测
多个气象站的温度数据插值
result_weather = makima([0, 2, 4, 6, 8, 10, 12], [[15, 16, 18, 22, 25, 28, 30], [14, 15, 17, 21, 24, 26, 28]], [1, 3, 5, 7, 9, 11])
print("多站点温度插值:")
print(result_weather)
#多站点温度插值:
#[[15.3125,16.8125,20,23.5,26.5,29.1875]
[14.3125,15.8125,18.95833333,22.66666667,25,27]]
# 同时插值多个气象站的数据
示例6：金融学 - 股票价格预测（带外推）
股票历史价格，使用外推预测未来
result_finance = makima([1, 2, 3, 4, 5], [100, 105, 98, 110, 115], [6, 7], 0)
print("股票价格预测:", result_finance)
#股票价格预测: [109.76923077 84.61538462]
#注意：实际金融应用中需要更复杂的模型
示例7：信号处理 - 不规则采样数据重采样
不规则时间点采集的信号，重采样到规则时间点
result_signal = makima([0.1, 0.5, 1.2, 1.8, 2.3, 3.0], [0.5, 0.8, 1.2, 0.9, 0.7, 0.6], [0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0])
print("信号重采样结果:", result_signal)
#信号重采样结果: [0.3984933,0.8,1.16126374,1.07279777,0.81515045,0.64728863,0.6]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.interpolate import CubicHermiteSpline
def makima_interp(input_str):
"""
MATLAB makima插值法的Python实现，支持符号输入解析
参数格式：
"(x, y)" -> 返回插值函数
"(x, y, xi)" -> 返回插值结果
"(x, y, xi, extrap)" -> 指定外推值
输入要求：
- x 需为单调递增序列
- x, y 长度相同
- xi 在x范围内或提供外推值
返回：
符号表达式/数值结果 或 错误信息
"""
try:
# 符号解析输入
expr = sp.sympify(input_str)
x, y, xi, extrap = None, None, None, None
if isinstance(expr, tuple):
if len(expr) < 2:
return "错误：至少需要x,y两个参数"
x = np.array([float(i) for i in expr[0]], dtype=np.float64)
y = np.array([float(i) for i in expr[1]], dtype=np.float64)
if len(expr) >= 3:
xi = np.array(expr[2], dtype=np.float64)
if len(expr) >= 4:
extrap = float(expr[3])
else:
return "错误：输入应为元组格式"
if len(x) != len(y):
return f"错误：x(length={len(x)})与y(length={len(y)})长度不匹配"
if not np.all(np.diff(x) > 0):
return "错误：x必须单调递增"
# 修正后的Akima斜率计算函数
def akima_slopes(x, y):
n = len(x)
if n == 1:
return np.zeros(1)
m = np.zeros(n)
if n == 2:
m[:] = (y[1] - y[0]) / (x[1] - x[0])
return m
delta = (y[1:] - y[:-1]) / np.diff(x)
if len(delta) == 0:
return m
# 扩展delta数组
d_m2 = 2 * delta[0] - delta[1]
d_m1 = 2 * d_m2 - delta[0]
d_n = 2 * delta[-1] - delta[-2]
d_n1 = 2 * d_n - delta[-1]
delta_ext = np.concatenate([[d_m2, d_m1], delta, [d_n, d_n1]])
# 计算每个点的斜率
for i in range(n):
w1 = np.abs(delta_ext[i + 3] - delta_ext[i + 2])
w2 = np.abs(delta_ext[i + 1] - delta_ext[i])
sum_w = w1 + w2
if sum_w == 0:
m[i] = (delta_ext[i + 1] + delta_ext[i + 2]) / 2
else:
m[i] = (w1 * delta_ext[i + 1] + w2 * delta_ext[i + 2]) / sum_w
return m
m = akima_slopes(x, y)
interp = CubicHermiteSpline(x, y, m)
n_intervals = len(x) - 1
# 符号表达式构造
t = sp.symbols('t')
x_sym = list(map(sp.Float, x))
pieces = []
for i in range(n_intervals):
coeffs = interp.c[:, i]
poly = (
coeffs[0] * (t - x_sym[i]) ** 3 +
coeffs[1] * (t - x_sym[i]) ** 2 +
coeffs[2] * (t - x_sym[i]) +
coeffs[3]
)
# 区间条件处理
if i == 0:
cond = (t >= x_sym[i]) & (t <= x_sym[i + 1])
else:
cond = (t > x_sym[i]) & (t <= x_sym[i + 1])
pieces.append((poly, cond))
final_expr = sp.Piecewise(*pieces)
if extrap is not None:
final_expr = sp.Piecewise(
(extrap, t < x_sym[0]),
(final_expr, (t >= x_sym[0]) & (t <= x_sym[-1])),
(extrap, t > x_sym[-1])
)
# 根据输入类型返回结果
if xi is not None:
if xi.dtype == np.dtype('O'):
return final_expr.subs(t, xi)
else:
# 数值计算时处理外推
if extrap is not None:
result = np.full_like(xi, extrap)
in_bounds = (xi >= x[0]) & (xi <= x[-1])
result[in_bounds] = interp(xi[in_bounds])
return result
else:
return interp(xi)
else:
return final_expr
except Exception as e:
return f"错误：{str(e)}"
if __name__ == "__main__":
# 测试1：基础插值（数值计算）
input1 = "([1, 2, 3, 4, 5], [1, 2, 1, 3, 2], [1.5, 2.5, 4.5])"
print("测试1 插值结果:", makima_interp(input1))
# [1.89285714 1.45714286 2.875 ]
# 测试2：符号表达式生成
input2 = "([0, 1, 2], [0, 1, 4])"
expr = makima_interp(input2)
print("\n测试2 符号表达式:")
print(expr)
# Piecewise((-0.666666666666667*t**3 + 2.66666666666667*t**2 - 1.0*t, (t <= 1.0) & (t >= 0)),
# (2.33333333333333*t + 0.333333333333334*(t - 1.0)**3 + 0.333333333333333*(t - 1.0)**2 - 1.33333333333333, (t <= 2.0) & (t > 1.0)))
# 测试3：外推测试
input3 = "([1, 3, 5], [1, 2, 1], [0, 6], 0)"
print("\n测试3 外推结果:", makima_interp(input3))
# [0. 0.]
求解线性分配问题
M = matchpairs(Cost) 求解矩阵 Cost 的行和列的线性分配问题. 该函数按总成本最小的方法将每行分配给一列.
M = matchpairs(Cost,costUnmatched) 求解矩阵 Cost 的行和列的线性分配问题. 该函数按总成本最小的方法将每行分配给一列. costUnmatched指定存在未分配行时的每行成本, 以及存在未分配行的列时的每列成本.
Cost — 成本矩阵,矩阵
costUnmatched — 未匹配成本,标量
M — 匹配项,矩阵
匹配，以矩阵形式返回. M 是 p×2 矩阵, 其中 M(i,1) 和 M(i,2) 是成本矩阵中匹配对组的行和列索引. M的行按第二列升序排序.
每行和每列只能匹配一次,因此每个 M(i,1) 值和每个 M(i,2) 值均唯一.
M 包含 p 个匹配, p 小于或等于最大匹配数 min(size(Cost)).
M 中的匹配成本是 sum([Cost(M(1,1),M(1,2)), Cost(M(2,1),M(2,2)), ..., Cost(M(p,1),M(p,2))]).
示例1：任务分配 - 最小化总完成时间
"""
任务分配问题：将任务分配给工人，最小化总完成时间
行：工人，列：任务，值：完成时间
"""
4个工人，4个任务
cost_matrix = [
[5, 8, 6, 7], # 工人1完成各任务的时间
[9, 6, 4, 5], # 工人2
[7, 5, 6, 8], # 工人3
[8, 7, 9, 6] # 工人4
]
result = matchpairs(cost_matrix)
print("任务分配结果 (工人, 任务):")
print(result)
#任务分配结果 (工人, 任务):
#[[1, 1],
[2, 3],
[3, 2],
[4, 4]]
示例2：车辆调度 - 出租车接客问题
"""
出租车调度：将乘客分配给最近的出租车
行：出租车位置，列：乘客位置，值：距离
"""
距离矩阵 (公里)
distance_matrix = [
[2.5, 1.8, 3.2, 4.1], # 出租车1到各乘客的距离
[3.1, 2.9, 1.5, 2.8], # 出租车2
[1.9, 3.5, 2.7, 3.9], # 出租车3
[4.2, 2.1, 3.8, 2.3] # 出租车4
]
设置最大可接受距离为3.5公里
max_distance = 3.5
result = matchpairs(distance_matrix, max_distance)
print("出租车调度结果 (出租车, 乘客):")
print(result)
#出租车调度结果 (出租车, 乘客):
#[[1, 2],
[2, 3],
[3, 1],
[4, 4]]
示例3：项目团队组建 - 技能匹配
"""
项目团队组建：将人员分配到项目角色，最大化技能匹配度
这里使用成本矩阵，值越小表示匹配度越高
"""
技能不匹配度矩阵 (1-10，1表示完美匹配，10表示完全不匹配)
skill_mismatch = [
[2, 5, 8, 3], # 人员1对各角色的不匹配度
[7, 3, 2, 6], # 人员2
[4, 6, 5, 4], # 人员3
[3, 4, 7, 2] # 人员4
]
最大可接受的不匹配度为5
max_mismatch = 5
result = matchpairs(skill_mismatch, max_mismatch)
print("团队组建结果 (人员, 角色):")
print(result)
#团队组建结果 (人员, 角色):
#[[1, 1],
[2, 3],
[3, 4],
[4, 2]]
示例4：生产调度 - 机器任务分配
"""
生产调度：将生产任务分配给机器，最小化总加工时间
"""
加工时间矩阵 (小时)
processing_time = [
[3, 6, 4, 5, 7], # 机器1加工各任务的时间
[5, 4, 3, 6, 5], # 机器2
[6, 5, 7, 4, 6], # 机器3
[4, 7, 5, 5, 4], # 机器4
[5, 6, 6, 7, 5] # 机器5
]
result = matchpairs(processing_time)
print("生产调度结果 (机器, 任务):")
print(result)
#生产调度结果 (机器, 任务):
#[[1, 1],
[2, 3],
[3, 4],
[4, 5],
[5, 2]]
示例5：学生选课 - 偏好最大化
"""
学生选课分配：基于学生偏好分配课程名额
成本值表示不满意度 (越低越好)
"""
学生选课偏好矩阵 (1-非常喜欢, 5-非常不喜欢)
preference_cost = [
[1, 3, 2, 4, 5], # 学生1对各课程的偏好成本
[2, 1, 4, 3, 5], # 学生2
[3, 2, 1, 5, 4], # 学生3
[4, 5, 3, 1, 2], # 学生4
[5, 4, 5, 2, 1] # 学生5
]
只接受偏好成本≤3的分配
max_cost = 3
result = matchpairs(preference_cost, max_cost)
print("选课分配结果 (学生, 课程):")
print(result)
#选课分配结果 (学生, 课程):
#[[1, 1],
[2, 2],
[3, 3],
[4, 4],
[5, 5]]
示例6：传感器-目标配对 - 军事应用
"""
传感器-目标配对：在军事应用中分配传感器监视目标
成本表示检测不确定性
"""
检测不确定性矩阵 (0-1，0表示完美检测，1表示完全不确定)
uncertainty_matrix = [
[0.1, 0.4, 0.2, 0.3], # 传感器1检测各目标的不确定性
[0.3, 0.1, 0.5, 0.2], # 传感器2
[0.2, 0.3, 0.1, 0.4], # 传感器3
[0.4, 0.2, 0.3, 0.1] # 传感器4
]
最大可接受的不确定性
max_uncertainty = 0.3
result = matchpairs(uncertainty_matrix, max_uncertainty)
print("传感器-目标配对结果 (传感器, 目标):")
print(result)
#传感器-目标配对结果 (传感器, 目标):
#[[1, 1],
[2, 2],
[3, 3],
[4, 4]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.optimize import linear_sum_assignment
def match_cost_assignment(C, max_cost=None):
"""
解决线性分配问题，对标MATLAB的matchpairs函数。
参数:
C (list, sympy.Matrix, numpy.ndarray): 成本矩阵，元素为数值。
max_cost (float, optional): 允许的最大成本，超过此值的边不会被选中。默认为None。
返回:
sympy.Matrix: 一个M×2的矩阵，每行表示一对匹配的行和列索引（1-based索引）。
异常:
当输入无法处理时返回错误信息字符串。
"""
try:
# 将输入转换为NumPy数组
if isinstance(C, sp.Matrix):
C_np = np.array(C.tolist(), dtype=float)
elif isinstance(C, (list, np.ndarray)):
C_np = np.array(C, dtype=float)
else:
raise ValueError("输入类型不支持，必须是列表、SymPy矩阵或NumPy数组。")
# 处理max_cost：将超过阈值的元素设为无穷大
if max_cost is not None:
C_np = np.where(C_np > max_cost, np.inf, C_np)
# 检查整个矩阵是否全为无穷大（无可行解）
if np.all(C_np == np.inf):
return sp.Matrix(0, 2, []) # 返回0行2列的空矩阵
# 使用匈牙利算法寻找最小权重匹配
row_ind, col_ind = linear_sum_assignment(C_np)
# 过滤掉成本为无穷大的匹配（无效边）
valid_mask = C_np[row_ind, col_ind] != np.inf
row_ind = row_ind[valid_mask]
col_ind = col_ind[valid_mask]
# 转换为1-based索引并生成结果矩阵
if row_ind.size > 0:
matches = np.column_stack((row_ind + 1, col_ind + 1))
else:
matches = np.empty((0, 2), dtype=int) # 无匹配时返回空矩阵
return sp.Matrix(matches.tolist())
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：基础用法
C1 = [[1, 3], [4, 2]]
print("示例1 - 无max_cost:")
print("匹配结果:\n", match_cost_assignment(C1))
# Matrix([[1, 1],
# [2, 2]])
# 示例2：使用max_cost过滤
C2 = [[5.7, 6.3, 3.1, 4.8, 3.5], [5.8, 6.4, 3.3, 4.7, 3.2], [5.7, 6.3, 3.2, 4.9, 3.4]]
max_cost2 = 0
print("\n示例2 - 使用max_cost=0:")
print("匹配结果:\n", match_cost_assignment(C2, max_cost2))
# Matrix(0, 2, [])
数组的均值
M = mean(A)返回A的所有元素的均值.
M = mean(A, 0)返回包含每列均值的行向量.
M = mean(A, 1)返回包含每行均值的行向量.
A - 输入数组, 向量, 矩阵.
示例1：分析学生各科成绩的平均分
学生成绩矩阵：4个学生，3门课程
行：学生，列：课程 [数学, 英语, 物理]
grades = [
[85, 92, 78],
[76, 88, 90],
[92, 85, 87],
[81, 79, 83]
]
学生成绩矩阵:
for i, student in enumerate(grades):
print(f"学生{i + 1}: 数学{student[0]}, 英语{student[1]}, 物理{student[2]}")
#学生1: 数学85, 英语92, 物理78
#学生2: 数学76, 英语88, 物理90
#学生3: 数学92, 英语85, 物理87
#学生4: 数学81, 英语79, 物理83
按列计算 - 各科目平均分
subject_avg = mean(str(grades))
print(f"\n各科目平均分: 数学{subject_avg[0]:.1f}, 英语{subject_avg[1]:.1f}, 物理{subject_avg[2]:.1f}")
#各科目平均分: 数学83.5, 英语86.0, 物理84.5
按行计算 - 每个学生的平均分
student_avg = mean(grades, 2)
print("\n学生个人平均分:")
for i, avg in enumerate(student_avg):
print(f"学生{i + 1}: {avg:.1f}分")
#学生个人平均分:
#学生1: 85.0分
#学生2: 84.7分
#学生3: 88.0分
#学生4: 81.0分
整体平均分
overall_avg = mean(np.array(grades).flatten())
print(f"\n班级整体平均分: {overall_avg:.1f}分")
#班级整体平均分: 84.7分
示例2：分析多个气象站的温度数据
5个气象站，7天的每日最高温度数据
行：气象站，列：天数
temperature_data = [
[25, 26, 24, 27, 28, 26, 25], # 站点1
[23, 24, 22, 25, 26, 24, 23], # 站点2
[27, 28, 26, 29, 30, 28, 27], # 站点3
[22, 23, 21, 24, 25, 23, 22], # 站点4
[26, 27, 25, 28, 29, 27, 26] # 站点5
]
各气象站周温度数据
days = ["周一", "周二", "周三", "周四", "周五", "周六", "周日"]
for i, station in enumerate(temperature_data):
print(f"站点{i + 1}: {station}")
#站点1: [25, 26, 24, 27, 28, 26, 25]
#站点2: [23, 24, 22, 25, 26, 24, 23]
#站点3: [27, 28, 26, 29, 30, 28, 27]
#站点4: [22, 23, 21, 24, 25, 23, 22]
#站点5: [26, 27, 25, 28, 29, 27, 26]
按列计算 - 每天的平均温度
daily_avg = mean(temperature_data)
每日平均温度:
for day, temp in zip(days, daily_avg):
print(f"{day}: {temp:.1f}°C")
#周一: 24.6°C
#周二: 25.6°C
#周三: 23.6°C
#周四: 26.6°C
#周五: 27.6°C
#周六: 25.6°C
#周日: 24.6°C
按行计算 - 各站点的周平均温度
station_avg = mean(temperature_data, 2)
各站点周平均温度:
for i, avg in enumerate(station_avg):
print(f"站点{i + 1}: {avg:.1f}°C")
#站点1: 25.9°C
#站点2: 23.9°C
#站点3: 27.9°C
#站点4: 22.9°C
#站点5: 26.9°C
区域整体平均温度
overall_avg = mean(np.array(temperature_data).flatten())
print(f"\n区域整体平均温度: {overall_avg:.1f}°C")
#区域整体平均温度: 25.5°C
示例3：分析多只股票的日收益率
4只股票，5个交易日的收益率 (%)
行：股票，列：交易日
returns = [
[1.2, -0.5, 2.3, -1.1, 0.8], # 股票A
[0.8, 1.5, -0.2, 0.9, 1.1], # 股票B
[-0.3, 2.1, 1.8, -0.7, 1.5], # 股票C
[1.5, -1.2, 0.7, 1.3, -0.4] # 股票D
]
stocks = ["股票A", "股票B", "股票C", "股票D"]
days = ["Day1", "Day2", "Day3", "Day4", "Day5"]
股票日收益率 (%):
for i, stock_returns in enumerate(returns):
print(f"{stocks[i]}: {stock_returns}")
#股票A: [1.2, -0.5, 2.3, -1.1, 0.8]
#股票B: [0.8, 1.5, -0.2, 0.9, 1.1]
#股票C: [-0.3, 2.1, 1.8, -0.7, 1.5]
#股票D: [1.5, -1.2, 0.7, 1.3, -0.4]
按列计算 - 每日平均收益率
daily_avg_returns = mean(returns)
每日平均收益率:
for day, ret in zip(days, daily_avg_returns):
print(f"{day}: {ret:.2f}%")
#Day1: 0.80%
#Day2: 0.48%
#Day3: 1.15%
#Day4: 0.10%
#Day5: 0.75%
按行计算 - 各股票的平均收益率
stock_avg_returns = mean(returns, 2)
各股票平均收益率:
for stock, avg_ret in zip(stocks, stock_avg_returns):
print(f"{stock}: {avg_ret:.2f}%")
#股票A: 0.54%
#股票B: 0.82%
#股票C: 0.88%
#股票D: 0.38%
投资组合整体平均收益率
overall_avg = mean(np.array(returns).flatten())
print(f"\n投资组合整体平均收益率: {overall_avg:.2f}%")
#投资组合整体平均收益率: 0.66%
示例4：分析生产批次的产品质量指标
4个生产批次，每个批次测量3个质量指标
[尺寸精度, 表面光洁度, 硬度]
quality_metrics = [
[98.5, 95.2, 87.3], # 批次1
[97.8, 96.1, 88.5], # 批次2
[99.2, 94.8, 86.9], # 批次3
[98.1, 95.9, 87.8] # 批次4
]
metrics = ["尺寸精度", "表面光洁度", "硬度"]
各批次产品质量指标:
for i, batch in enumerate(quality_metrics):
print(f"批次{i + 1}: {batch}")
#批次1: [98.5, 95.2, 87.3]
#批次2: [97.8, 96.1, 88.5]
#批次3: [99.2, 94.8, 86.9]
#批次4: [98.1, 95.9, 87.8]
按列计算 - 各质量指标的平均值
metric_avg = mean(quality_metrics)
各质量指标平均值:
for metric, avg in zip(metrics, metric_avg):
print(f"{metric}: {avg:.1f}%")
#尺寸精度: 98.4%
#表面光洁度: 95.5%
#硬度: 87.6%
按行计算 - 各批次的综合质量评分
batch_avg = mean(quality_metrics, 2)
各批次综合质量评分:
for i, avg in enumerate(batch_avg):
print(f"批次{i + 1}: {avg:.1f}%")
#批次1: 93.7%
#批次2: 94.1%
#批次3: 93.6%
#批次4: 93.9%
整体质量水平
overall_avg = mean(np.array(quality_metrics).flatten())
print(f"\n整体质量水平: {overall_avg:.1f}%")
#整体质量水平: 93.8%
示例5：分析多个门店的销售数据
3个门店，4个季度的销售额（万元）
sales_data = [
[120, 150, 180, 140], # 门店A
[90, 110, 130, 100], # 门店B
[200, 180, 220, 190] # 门店C
]
stores = ["门店A", "门店B", "门店C"]
quarters = ["Q1", "Q2", "Q3", "Q4"]
各门店季度销售额（万元）:
for i, sales in enumerate(sales_data):
print(f"{stores[i]}: {sales}")
#门店A: [120, 150, 180, 140]
#门店B: [90, 110, 130, 100]
#门店C: [200, 180, 220, 190]
按列计算 - 各季度平均销售额
quarterly_avg = mean(sales_data)
各季度平均销售额:
for q, avg in zip(quarters, quarterly_avg):
print(f"{q}: {avg:.1f}万元")
#Q1: 136.7万元
#Q2: 146.7万元
#Q3: 176.7万元
#Q4: 143.3万元
按行计算 - 各门店年平均销售额
store_avg = mean(sales_data, 2)
各门店年平均销售额:
for store, avg in zip(stores, store_avg):
print(f"{store}: {avg:.1f}万元")
#门店A: 147.5万元
#门店B: 107.5万元
#门店C: 197.5万元
公司整体平均季度销售额
overall_avg = mean(np.array(sales_data).flatten())
print(f"\n公司整体平均季度销售额: {overall_avg:.1f}万元")
#公司整体平均季度销售额: 150.8万元
示例6：分析科学实验的重复测量数据
3个实验条件，每个条件5次重复测量的结果
experimental_data = [
[12.3, 11.8, 12.1, 12.5, 11.9], # 条件A
[15.6, 16.1, 15.8, 15.3, 16.2], # 条件B
[9.8, 10.2, 9.9, 10.1, 10.0] # 条件C
]
conditions = ["条件A", "条件B", "条件C"]
实验重复测量数据:
for i, measurements in enumerate(experimental_data):
print(f"{conditions[i]}: {measurements}")
#条件A: [12.3, 11.8, 12.1, 12.5, 11.9]
#条件B: [15.6, 16.1, 15.8, 15.3, 16.2]
#条件C: [9.8, 10.2, 9.9, 10.1, 10.0]
按行计算 - 各条件的平均值
condition_means = mean(experimental_data, 2)
各实验条件的平均值:
for cond, mean_val in zip(conditions, condition_means):
print(f"{cond}: {mean_val:.2f}")
#条件A: 12.12
#条件B: 15.80
#条件C: 10.00
计算标准差（使用numpy）
condition_stds = [np.std(measurements) for measurements in experimental_data]
各实验条件的标准差:
for cond, std_val in zip(conditions, condition_stds):
print(f"{cond}: {std_val:.3f}")
#条件A: 0.256
#条件B: 0.329
#条件C: 0.141
所有测量值的总平均
overall_mean = mean(np.array(experimental_data).flatten())
print(f"\n所有测量值的总平均: {overall_mean:.2f}")
#所有测量值的总平均: 12.64
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def mean_cal_array(input_str):
"""
对标 MATLAB 的 mean 函数，计算数组的平均值
参数:
input_str (str): 输入的矩阵表达式，可以是：
- 普通矩阵："[[1,2],[3,4]]"
- 带维度参数的元组："([[1,2],[3,4]], 2)"
返回:
float/list/str: 计算结果（标量/列表）或错误信息
MATLAB 对应行为:
- mean(A) 默认按列计算
- mean(A,1) 按列计算
- mean(A,2) 按行计算
- 对向量返回整体平均值
"""
try:
expr = sp.sympify(input_str)
result = None
error = False
# 情况1：带维度参数的输入 (矩阵, 维度)
if isinstance(expr, tuple):
# 验证元组格式 (必须包含2个元素，第二个元素为整数)
if len(expr) != 2 or not expr[1].is_Integer:
error = True
else:
M = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
dim = int(expr[1])
# 检查矩阵合法性和维度值有效性
if M is not None and dim in [1, 2]:
data = np.array(M, dtype=np.float64)
axis = dim - 1 # 将 MATLAB 维度转换为 numpy 轴
result = np.mean(data, axis=axis)
else:
error = True
# 情况2：单个矩阵输入
elif isinstance(expr, list):
data = np.array(expr, dtype=np.float64)
# 判断是否为向量（行向量或列向量）
if data.ndim == 1 or 1 in data.shape:
result = np.mean(data) # 整体平均值
else:
# 默认按列计算（MATLAB 的默认行为）
result = np.mean(data, axis=0)
else:
error = True
# 错误处理
if error:
return f"输入错误: {input_str}"
else:
# 将 numpy 类型转换为 Python 原生类型
return result.item() if np.isscalar(result) else result.tolist()
except Exception as e:
return f"Error: {e}"
# 示例测试
if __name__ == "__main__":
# 示例1：矩阵默认按列计算
print(mean_cal_array("[[1,2], [3,4]]"))
# [2.0, 3.0]
# 示例2：按行计算平均值
print(mean_cal_array("([[1,2], [3,4]], 2)"))
# [1.5, 3.5]
# 示例3：行向量计算整体平均值
print(mean_cal_array("[[1,2,3,4]]"))
# 2.5
# 示例4：列向量计算整体平均值
print(mean_cal_array("[[1], [2], [3], [4]]"))
# 2.5
# 示例6：三维矩阵测试
print(mean_cal_array("[[[1,2],[3,4]]]"))
# 2.5
数组的中位值
M = median(A)返回A的所有元素的中位值.
M = median(A, 0)返回包含每列中位值的行向量.
M = median(A, 1)返回包含每行中位值的行向量.
A - 输入数组, 向量, 矩阵.
示例1：收入分析 - 避免极端值影响
5个地区，每个地区抽样7个家庭的年收入（万元）
注意：包含一些极端高收入，平均数会被拉高
income_data = [
[15, 18, 20, 22, 25, 28, 200], # 地区A（有一个超高收入）
[12, 15, 16, 18, 20, 22, 24], # 地区B
[20, 22, 25, 28, 30, 35, 40], # 地区C
[8, 10, 12, 15, 18, 20, 150], # 地区D（有一个超高收入）
[25, 28, 30, 32, 35, 38, 40] # 地区E
]
regions = ["地区A", "地区B", "地区C", "地区D", "地区E"]
各地区家庭年收入数据（万元）:
for i, incomes in enumerate(income_data):
print(f"{regions[i]}: {incomes}")
#地区A: [15, 18, 20, 22, 25, 28, 200]
#地区B: [12, 15, 16, 18, 20, 22, 24]
#地区C: [20, 22, 25, 28, 30, 35, 40]
#地区D: [8, 10, 12, 15, 18, 20, 150]
#地区E: [25, 28, 30, 32, 35, 38, 40]
按行计算中位数 - 各地区收入中位数
region_medians = median(income_data, 2)
各地区收入中位数:
for region, median_income in zip(regions, region_medians):
print(f"{region}: {median_income}万元")
#地区A: 22.0万元
#地区B: 18.0万元
#地区C: 28.0万元
#地区D: 15.0万元
#地区E: 32.0万元
对比平均数（受极端值影响）
region_means = [np.mean(incomes) for incomes in income_data]
对比平均数（受极端值影响）:
for region, mean_income, median_income in zip(regions, region_means, region_medians):
difference = mean_income - median_income
print(f"{region}: 平均数{mean_income:.1f}万, 中位数{median_income}万, 差异{difference:.1f}万")
#地区A: 平均数46.9万, 中位数22.0万, 差异24.9万
#地区B: 平均数18.1万, 中位数18.0万, 差异0.1万
#地区C: 平均数28.6万, 中位数28.0万, 差异0.6万
#地区D: 平均数33.3万, 中位数15.0万, 差异18.3万
#地区E: 平均数32.6万, 中位数32.0万, 差异0.6万
示例2：房价分析
4个区域，每个区域的房价样本（万元/平米）
housing_prices = [
[3.5, 4.2, 5.1, 6.8, 8.9, 12.5], # 市中心（价格差异大）
[2.8, 3.2, 3.5, 3.8, 4.1, 4.3], # 郊区（价格相对集中）
[1.5, 1.8, 2.2, 2.5, 15.0, 18.0], # 新区（有高端楼盘拉高）
[4.0, 4.5, 5.0, 5.5, 6.0, 6.5] # 学区房（价格稳定）
]
areas = ["市中心", "郊区", "新区", "学区房"]
各区域房价样本（万元/平米）:
for i, prices in enumerate(housing_prices):
print(f"{areas[i]}: {prices}")
#市中心: [3.5, 4.2, 5.1, 6.8, 8.9, 12.5]
#郊区: [2.8, 3.2, 3.5, 3.8, 4.1, 4.3]
#新区: [1.5, 1.8, 2.2, 2.5, 15.0, 18.0]
#学区房: [4.0, 4.5, 5.0, 5.5, 6.0, 6.5]
计算各区域房价中位数
area_medians = median(housing_prices, 2)
各区域房价中位数:
for area, median_price in zip(areas, area_medians):
print(f"{area}: {median_price}万元/平米")
#市中心: 5.95万元/平米
#郊区: 3.65万元/平米
#新区: 2.35万元/平米
#学区房: 5.25万元/平米
对比分析
area_means = [np.mean(prices) for prices in housing_prices]
中位数 vs 平均数分析:
for area, mean_price, median_price in zip(areas, area_means, area_medians):
if mean_price > median_price:
skew = "右偏（有高价房拉高）"
elif mean_price < median_price:
skew = "左偏（有低价房拉低）"
else:
skew = "对称分布"
print(f"{area}: 中位数{median_price}, 平均数{mean_price:.1f} -> {skew}")
#市中心: 中位数5.949999999999999, 平均数6.8 -> 右偏（有高价房拉高）
#郊区: 中位数3.65, 平均数3.6 -> 左偏（有低价房拉低）
#新区: 中位数2.35, 平均数6.8 -> 右偏（有高价房拉高）
#学区房: 中位数5.25, 平均数5.2 -> 对称分布
示例3：响应时间分析
3个系统，每个系统10次请求的响应时间（毫秒）
偶尔会有异常高的响应时间（网络波动等）
response_times = [
[45, 48, 52, 55, 58, 60, 62, 65, 68, 250], # 系统A（有一次超时）
[38, 40, 42, 45, 47, 50, 52, 55, 58, 60], # 系统B
[55, 58, 60, 62, 65, 68, 70, 75, 80, 300] # 系统C（有一次超时）
]
systems = ["系统A", "系统B", "系统C"]
各系统响应时间（毫秒）:
for i, times in enumerate(response_times):
print(f"{systems[i]}: {times}")
#系统A: [45, 48, 52, 55, 58, 60, 62, 65, 68, 250]
#系统B: [38, 40, 42, 45, 47, 50, 52, 55, 58, 60]
#系统C: [55, 58, 60, 62, 65, 68, 70, 75, 80, 300]
计算各系统中位数响应时间
system_medians = median(response_times, 2)
各系统中位数响应时间:
for system, median_time in zip(systems, system_medians):
print(f"{system}: {median_time}毫秒")
#系统A: 59.0毫秒
#系统B: 48.5毫秒
#系统C: 66.5毫秒
对比分析
system_means = [np.mean(times) for times in response_times]
性能分析（中位数更能反映典型性能）:
for system, mean_time, median_time in zip(systems, system_means, system_medians):
print(f"{system}: 中位数{median_time}ms, 平均数{mean_time:.1f}ms")
if mean_time > median_time + 10:
print(f" → 有异常高响应时间影响平均数")
#系统A: 中位数59.0ms, 平均数76.3ms
# → 有异常高响应时间影响平均数
#系统B: 中位数48.5ms, 平均数48.7ms
#系统C: 中位数66.5ms, 平均数89.3ms
# → 有异常高响应时间影响平均数
示例4：考试成绩分析
4个班级，每个班级学生的考试成绩
exam_scores = [
[45, 58, 62, 65, 68, 72, 75, 78, 82, 85, 88, 92, 95], # 班级A
[35, 42, 55, 58, 62, 65, 68, 72, 75, 78, 82, 85, 95], # 班级B
[75, 78, 80, 82, 85, 88, 90, 92, 95, 98, 55, 40, 35], # 班级C（有几个低分）
[60, 62, 65, 68, 70, 72, 75, 78, 80, 82, 85, 88, 90] # 班级D
]
classes = ["班级A", "班级B", "班级C", "班级D"]
各班级考试成绩:
for i, scores in enumerate(exam_scores):
print(f"{classes[i]}: {scores}")
#班级A: [45, 58, 62, 65, 68, 72, 75, 78, 82, 85, 88, 92, 95]
#班级B: [35, 42, 55, 58, 62, 65, 68, 72, 75, 78, 82, 85, 95]
#班级C: [75, 78, 80, 82, 85, 88, 90, 92, 95, 98, 55, 40, 35]
#班级D: [60, 62, 65, 68, 70, 72, 75, 78, 80, 82, 85, 88, 90]
计算各班级成绩中位数
class_medians = median(exam_scores, 2)
各班级成绩中位数:
for cls, median_score in zip(classes, class_medians):
print(f"{cls}: {median_score}分")
#班级A: 75.0分
#班级B: 68.0分
#班级C: 82.0分
#班级D: 75.0分
对比分析
class_means = [np.mean(scores) for scores in exam_scores]
教学效果分析:
for cls, mean_score, median_score in zip(classes, class_means, class_medians):
if mean_score < median_score:
analysis = "低分学生较多"
elif mean_score > median_score:
analysis = "高分学生较多"
else:
analysis = "分布相对对称"
print(f"{cls}: 中位数{median_score}, 平均数{mean_score:.1f} -> {analysis}")
#班级A: 中位数75.0, 平均数74.2 -> 低分学生较多
#班级B: 中位数68.0, 平均数67.1 -> 低分学生较多
#班级C: 中位数82.0, 平均数76.4 -> 低分学生较多
#班级D: 中位数75.0, 平均数75.0 -> 分布相对对称
示例5：客户消费分析
4个产品类别，每个类别的客户消费金额（元）
spending_data = [
[50, 80, 120, 150, 180, 200, 250, 300, 500, 1000], # 电子产品（消费差异大）
[30, 35, 40, 45, 50, 55, 60, 65, 70, 75], # 日用品（消费集中）
[100, 120, 150, 180, 200, 250, 300, 400, 600, 50], # 服装（有低价清仓）
[200, 220, 250, 280, 300, 320, 350, 380, 400, 450] # 奢侈品（消费稳定）
]
categories = ["电子产品", "日用品", "服装", "奢侈品"]
各产品类别客户消费金额（元）:
for i, spendings in enumerate(spending_data):
print(f"{categories[i]}: {spendings}")
#电子产品: [50, 80, 120, 150, 180, 200, 250, 300, 500, 1000]
#日用品: [30, 35, 40, 45, 50, 55, 60, 65, 70, 75]
#服装: [100, 120, 150, 180, 200, 250, 300, 400, 600, 50]
#奢侈品: [200, 220, 250, 280, 300, 320, 350, 380, 400, 450]
计算各类别消费金额中位数
category_medians = median(spending_data, 2)
各产品类别消费金额中位数:
for category, median_spending in zip(categories, category_medians):
print(f"{category}: {median_spending}元")
#电子产品: 190.0元
#日用品: 52.5元
#服装: 190.0元
#奢侈品: 310.0元
对比分析
category_means = [np.mean(spendings) for spendings in spending_data]
消费行为分析:
for category, mean_spending, median_spending in zip(categories, category_means, category_medians):
ratio = mean_spending / median_spending
if ratio > 1.5:
pattern = "高端消费主导"
elif ratio < 0.7:
pattern = "低价消费主导"
else:
pattern = "消费分布均匀"
print(f"{category}: 中位数{median_spending}元, 平均数{mean_spending:.1f}元 -> {pattern}")
#电子产品: 中位数190.0元, 平均数283.0元 -> 消费分布均匀
#日用品: 中位数52.5元, 平均数52.5元 -> 消费分布均匀
#服装: 中位数190.0元, 平均数235.0元 -> 消费分布均匀
#奢侈品: 中位数310.0元, 平均数315.0元 -> 消费分布均匀
示例6：医疗数据分析
3种治疗方案，患者的恢复时间（天）
recovery_times = [
[5, 6, 7, 8, 9, 10, 12, 15, 20, 45], # 方案A（有一个异常恢复）
[8, 9, 10, 11, 12, 13, 14, 15, 16, 18], # 方案B
[6, 7, 8, 9, 10, 11, 12, 30, 35, 40] # 方案C（有几个异常恢复）
]
treatments = ["治疗方案A", "治疗方案B", "治疗方案C"]
各治疗方案患者恢复时间（天）:
for i, times in enumerate(recovery_times):
print(f"{treatments[i]}: {times}")
#治疗方案A: [5, 6, 7, 8, 9, 10, 12, 15, 20, 45]
#治疗方案B: [8, 9, 10, 11, 12, 13, 14, 15, 16, 18]
#治疗方案C: [6, 7, 8, 9, 10, 11, 12, 30, 35, 40]
计算各治疗方案恢复时间中位数
treatment_medians = median(recovery_times, 2)
各治疗方案恢复时间中位数:
for treatment, median_time in zip(treatments, treatment_medians):
print(f"{treatment}: {median_time}天")
#治疗方案A: 9.5天
#治疗方案B: 12.5天
#治疗方案C: 10.5天
对比分析
treatment_means = [np.mean(times) for times in recovery_times]
治疗效果评估（中位数更可靠）:
for treatment, mean_time, median_time in zip(treatments, treatment_means, treatment_medians):
print(f"{treatment}: 中位数{median_time}天, 平均数{mean_time:.1f}天")
if mean_time > median_time + 5:
print(f" → 有异常恢复时间病例")
#治疗方案A: 中位数9.5天, 平均数13.7天
#治疗方案B: 中位数12.5天, 平均数12.6天
#治疗方案C: 中位数10.5天, 平均数16.8天
# → 有异常恢复时间病例
异常值识别:
for i, times in enumerate(recovery_times):
q75, q25 = np.percentile(times, [75, 25])
iqr = q75 - q25
upper_bound = q75 + 1.5 * iqr
outliers = [t for t in times if t > upper_bound]
if outliers:
print(f"{treatments[i]}: 异常恢复时间 {outliers}")
#治疗方案A: 异常恢复时间 [45]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def median_cal_array(input_str):
"""
对标 MATLAB 的 median 函数，计算矩阵的中位数。
参数:
input_str (str): 输入的矩阵表达式，可以是普通矩阵或带维度参数的元组。
示例: "[[1,2],[3,4]]" 或 "([[1,2],[3,4]], 2)"
返回:
如果输入合法，返回中位数计算结果（标量或列表）；否则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
# 情况1：输入是元组 (矩阵, 维度)
if isinstance(expr, tuple):
# 检查元组长度和元素类型
if len(expr) != 2 or not expr[1].is_Integer:
error = True
else:
M = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
dim = int(expr[1])
# 检查矩阵合法性和维度值 (MATLAB 风格，1=列，2=行)
if M is not None and dim in [1, 2]:
data = np.array(M, dtype=np.float64)
axis = dim - 1 # 转换为 numpy 的轴 (0=列，1=行)
result = np.median(data, axis=axis)
else:
error = True
# 情况2：输入是单个矩阵
elif isinstance(expr, list):
data = np.array(expr, dtype=np.float64)
# 判断是否为向量（行向量或列向量）
if data.ndim == 1 or 1 in data.shape:
result = np.median(data) # 整个向量的中位数
else:
result = np.median(data, axis=0) # 默认按列计算
else:
error = True
# 错误处理
if error:
return f"输入错误: {input_str}"
else:
# 将 numpy 结果转换为 Python 原生类型
return result.item() if np.isscalar(result) else result.tolist()
except Exception as e:
return f"Error: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1：矩阵按列计算（默认）
print(median_cal_array("[[1,2,3], [4,5,6]]"))
# [2.5, 3.5, 4.5]
# 示例2：按行计算
print(median_cal_array("([[1,2,3], [4,5,6]], 2)"))
# [2.0, 5.0]
# 示例3：行向量默认计算
print(median_cal_array("[[1,2,3]]"))
# 2.0
# 示例4：列向量默认计算
print(median_cal_array("[[1], [2], [3]]"))
# 2.0
MedianFit([x], y)：使用中位数-中位数直线法将数据拟合为形式 y = a*x + b 的线性方程。
[x]（可选）-- 自变量，向量。
y -- 因变量，向量，与x维度一致。
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
import re
import ast
def parse_list_string(input_string):
"""解析输入字符串，提取所有列表格式的子串"""
# 正则匹配方括号内的内容（非贪婪模式）
pattern = re.compile(r'\[.*?\]')
matches = pattern.findall(input_string)
return matches if matches else None
def evaluate_list(data):
"""将字符串列表转换为Python对象（支持数值、符号表达式）"""
try:
# 尝试直接解析为Python列表
return ast.literal_eval(data)
except (ValueError, SyntaxError):
try:
# 尝试作为表达式解析（如包含SymPy符号）
return eval(data)
except NameError:
# 作为符号表达式解析
return sp.sympify(data)
def median_fit(data_lists):
"""中位拟合：y = a*x + b，基于中位数的鲁棒回归"""
# 将输入数据转换为合法列表（支持单列表或双列表输入）
data_lists = [evaluate_list(data) for data in data_lists]
# 处理输入数据
if len(data_lists) == 1:
# 单列表输入：y值，x自动生成索引 [0,1,2,...]
y = np.array(data_lists[0], dtype=np.float64)
x = np.arange(len(y))
elif len(data_lists) >= 2:
# 双列表输入：x和y值
x = np.array(data_lists[0], dtype=np.float64)
y = np.array(data_lists[1], dtype=np.float64)
# 检查长度匹配
if len(x) != len(y):
raise ValueError("x和y的长度必须相同")
else:
raise ValueError("至少需要提供一组数据")
# 计算中位数坐标点 (x_med, y_med)
x_med, y_med = np.median(x), np.median(y)
# 计算斜率a = median[(x_i - x_med)(y_i - y_med)] / median[(x_i - x_med)^2]
numerator = np.median((x - x_med) * (y - y_med))
denominator = np.median((x - x_med) ** 2)
if denominator == 0:
raise ZeroDivisionError("分母为零，无法计算斜率（可能所有x值相同）")
a = numerator / denominator
# 计算截距b = y_med - a*x_med
b = y_med - a * x_med
# 构建SymPy表达式
x_sym = sp.symbols('x')
equation = a * x_sym + b
return equation
def median_fit_expression(input_str):
"""对外接口：解析输入字符串并返回拟合方程"""
lists = parse_list_string(input_str)
if not lists:
raise ValueError("输入格式错误，应包含至少一个列表")
return median_fit(lists)
# 示例验证
if __name__ == "__main__":
# 示例1：单列表输入（自动生成x索引）
y_data = [3, 4, 5, 7, 8, 10, 4, 7, 6]
eq1 = median_fit_expression("[3, 4, 5, 7, 8, 10, 4, 7, 6]")
print(f"示例1拟合方程: {eq1}")
# 0.5*x + 4.0
# 示例2：双列表输入（指定x和y）
x_data = [1, 4, 9, 5, 7, 5, 4, 2, 9]
y_data = [3, 4, 5, 7, 8, 10, 4, 7, 6]
eq2 = median_fit_expression(f"{x_data}, {y_data}")
print(f"示例2拟合方程: {eq2}")
# 6.00000000000000
# 对比普通最小二乘法（numpy.polyfit）
a_ls, b_ls = np.polyfit(x_data, y_data, 1)
print(f"最小二乘法结果: y = {a_ls:.3f}x + {b_ls:.3f}")
# y = 0.207x + 4.943
a, b = MedianFitModel([x], y)：使用中位数-中位数直线法将数据拟合为线性模型 y = a \cdot x + by=a⋅x+b，并返回斜率 a 和截距 b。
[x] -- 可选, 自变量，向量。
y -- 因变量（目标）的观测值，长度必须与 x 一致。
a -- 斜率，表示x每增加1单位时y的中位数变化量
b -- 截距，表示x=0时y的中位数基准值
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
import re
import ast
def parse_list_string(input_string):
"""解析输入字符串，提取所有列表格式的子串"""
# 正则匹配方括号内的内容（非贪婪模式）
pattern = re.compile(r'\[.*?\]')
matches = pattern.findall(input_string)
return matches if matches else None
def evaluate_list(data):
"""将字符串列表转换为Python对象（支持数值、符号表达式）"""
try:
# 尝试直接解析为Python列表
return ast.literal_eval(data)
except (ValueError, SyntaxError):
try:
# 尝试作为表达式解析（如包含SymPy符号）
return eval(data)
except NameError:
# 作为符号表达式解析
return sp.sympify(data)
# --- 求a, b, 中位拟合 start
def median_fit_with_correlation(data_lists):
data_lists = [evaluate_list(data) for data in data_lists]
if len(data_lists) == 1:
y_data = data_lists[0]
x_data = np.arange(len(y_data))
elif len(data_lists) > 1:
x_data = data_lists[0]
y_data = data_lists[1]
x_med = np.median(x_data)
y_med = np.median(y_data)
a = np.median((x_data - x_med) * (y_data - y_med)) / np.median((x_data - x_med) ** 2)
b = y_med - a * x_med
# Returning the result as a formatted string
result_string = f'a: {a}, b: {b}'
return result_string
def median_fit_model_expression(input_str):
# 判断元素是列表还是随机函数
res_expression = parse_list_string(input_str)
if res_expression:
return median_fit_with_correlation(res_expression)
else:
print(f"{input_str} 既不是列表字符串")
# --- 求a, b 中位拟合 end
# 示例验证
if __name__ == "__main__":
# 示例1：单列表输入（自动生成x索引）
y_data = [3, 4, 5, 7, 8, 10, 4, 7, 6]
eq1 = median_fit_model_expression("[3, 4, 5, 7, 8, 10, 4, 7, 6]")
print(f"示例1拟合方程: {eq1}")
# a: 0.5, b: 4.0
# 示例2：双列表输入（指定x和y）
x_data = [1, 4, 9, 5, 7, 5, 4, 2, 9]
y_data = [3, 4, 5, 7, 8, 10, 4, 7, 6]
eq2 = median_fit_model_expression(f"{x_data}, {y_data}")
print(f"示例2拟合方程: {eq2}")
# a: 0.0, b: 6.0
求解关于x的线性方程组
x = A\B 对线性方程组A*x=B求解.矩阵A和B必须具有相同的行数.如果A未正确缩放或接近奇异值,还是会执行计算.
A, B — 操作数,向量,满矩阵,稀疏矩阵
x — 解,向量,满矩阵,稀疏矩阵
示例 1: 电路分析 - 求解节点电压
电路方程: 2V1 - V2 = 5, -V1 + 3V2 = 2
A_circuit = [[2, -1], [-1, 3]]
B_circuit = [[5], [2]]
print("电路方程解:", mldivide(A_circuit, B_circuit))
#电路方程解:
#[[3.4],
[1.8]]
示例 2: 力学系统 - 弹簧质点系统
弹簧系统: 3x + 2y = 10, x + 4y = 8
A_spring = [[3, 2], [1, 4]]
B_spring = [[10], [8]]
print("弹簧系统解:", mldivide(A_spring,B_spring))
#弹簧系统解:
#[[2.4],
[1.4]]
示例 3: 经济学 - 市场均衡
供需方程: 2P + Q = 20, P - Q = 4
A_econ = [[2, 1], [1, -1]]
B_econ = [[20], [4]]
print("市场均衡解:", mldivide(A_econ,B_econ))
#市场均衡解:
#[[8],
[4]]
示例 4: 化学计量 - 化学反应平衡
反应: aCH4 + bO2 → cCO2 + dH2O
原子平衡: C: a = c, H: 4a = 2d, O: 2b = 2c + d
A_chem = [[1, 0, -1, 0], [4, 0, 0, -2], [0, 2, -2, -1]]
B_chem = [[0], [0], [0]]
print("化学反应系数:", mldivide(A_chem,B_chem))
#化学反应系数:
#[[0],
[0],
[0],
[0]]
示例 5: 数据拟合 - 线性回归
拟合 y = ax + b
数据点: (1,2), (2,3), (3,5), (4,6)
A_fit = [[1, 1], [1, 2], [1, 3], [1, 4]]
B_fit = [[2], [3], [5], [6]]
print("线性回归系数 [b, a]:", mldivide(A_fit,B_fit))
#线性回归系数 [b, a]:
#[[0.500000000000001],
[1.4]]
示例 6: 结构分析 - 桁架受力
节点平衡方程
A_struct = [[1, 0, 1, 0], [0, 1, 0, 1], [-1, 0, 0, 0], [0, -1, 0, 0]]
B_struct = [[0], [10], [5], [0]]
print("桁架受力解:", mldivide(A_struct,B_struct))
#桁架受力解:
#[[-5],
[0],
[5],
[10]]
示例 7: 网络流 - 交通流量
交通网络流量平衡
A_flow = [[1, -1, 0, 0], [0, 1, -1, 0], [1, 0, 0, -1]]
B_flow = [[100], [50], [50]]
print("交通流量解:", mldivide(A_flow, B_flow))
#交通流量解:
#[[75],
[-25],
[-75],
[25]]
示例 8: 热传导 - 温度分布
一维热传导离散方程
A_heat = [[2, -1, 0], [-1, 2, -1], [0, -1, 2]]
B_heat = [[100], [0], [50]]
print("温度分布解:", mldivide(A_heat, B_heat))
#温度分布解:
#[[87.5],
[75],
[62.5]]
示例 9: 机器人运动学 - 逆运动学
简化的机器人关节角度求解
A_robot = [[1, 0.5], [0.5, 1]]
B_robot = [[1.2], [0.8]]
print("关节角度解:", min_linsolvedivide(A_robot, B_robot))
#关节角度解:
#[[1.06666666666667],
[0.266666666666667]]
示例 10: 金融 - 资产组合权重
最小方差组合
A_finance = [[0.1, 0.03, 1], [0.03, 0.15, 1], [1, 1, 0]]
B_finance = [[0], [0], [1]]
print("资产权重解:", mldivide(A_finance, B_finance))
#资产权重解:
#[[0.631578947368421],
[0.368421052631579],
[-0.0742105263157895]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
from scipy import linalg as dense_linalg
import sympy as sp
import numpy as np
def min_linsolvedivide(input_str):
"""
对标 MATLAB 的 mldivide 函数，求解线性方程组 Ax = B。
参数:
input_str: 字符串形式的输入，应为包含两个矩阵的元组，例如 "(Matrix([[1,2],[3,4]]), Matrix([[5],[6]]))"
返回:
SymPy 矩阵形式的结果，或错误信息。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
A = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
B = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
if A is not None and B is not None:
# 将 SymPy 矩阵转换为 NumPy 二维数组，保持结构
N_A = np.array(A.tolist(), dtype=float)
N_B = np.array(B.tolist(), dtype=float)
try:
# 尝试直接求解 Ax = B（适用于方阵且可逆的情况）
result = np.linalg.solve(N_A, N_B)
except np.linalg.LinAlgError:
# 若无法直接求解，使用最小二乘法（处理超定、欠定等情况）
result, _, _, _ = np.linalg.lstsq(N_A, N_B, rcond=None)
else:
error = True
else:
error = True
if error:
return f"输入错误: 输入应为包含两个矩阵的元组。"
else:
return sp.Matrix(result)
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例 1: 方阵可逆
A = [[1, 2], [3, 4]]
B = [[5], [6]]
input_str = f"({A}, {B})"
print("示例1 解:", min_linsolvedivide(input_str))
# Matrix([[-4.00000000000000],
# [4.50000000000000]])
# 示例 2: 超定系统（最小二乘解）
A = [[1], [2], [3]]
B = [[4], [5], [6]]
input_str = f"({A}, {B})"
print("示例2 解:", min_linsolvedivide(input_str))
# Matrix([[2.28571428571429]])
# 示例 3: 欠定系统（最小范数解）
A = [[1, 1, 1], [2, 2, 2]]
B = [[3], [6]]
input_str = f"({A}, {B})"
print("示例3 解:", min_linsolvedivide(input_str))
# Matrix([[1.00000000000000],
# [1.00000000000000],
# [1.00000000000000]])
生成分段多项式
pp = mkpp(breaks,coefs) 根据其间断数和系数生成分段多项式 pp。使用 ppval 计算特定点处的分段多项式，或使用 unmkpp 提取有关分段多项式的详细信息。
breaks — 断点,向量
coefs — 多项式系数,矩阵
示例 1: 三次样条插值（机械工程）
描述机械臂轨迹的平滑运动
spline_pp = mkpp(
[0, 1, 2, 3],
[[1, 0, -2, 1], # 在 [0,1]: x^3 - 2x + 1
[0.5, -1, 1, 0], # 在 [1,2]: 0.5x^3 - x^2 + x
[-1, 3, -2, 1]]) # 在 [2,3]: -x^3 + 3x^2 - 2x + 1
)
print("三次样条分段多项式:")
print(spline_pp)
#三次样条分段多项式:
#Piecewise((x**3 - 2*x + 1, (x >= 0) & (x < 1)), (0.5*x**3 - x**2 + x, (x >= 1) & (x < 2)),
(-x**3 + 3*x**2 - 2*x + 1, (x >= 2) & (x <= 3)), (0, True))
示例 2: 温度控制分段策略（过程控制）
不同温度区间的加热功率控制
temp_control = mkpp(
[0, 50, 100, 150, 200],
[[0.02, 0], # 0-50°C: 线性加热 0.02x
[0.1, -2], # 50-100°C: 0.1x - 2
[0.05, 3], # 100-150°C: 0.05x + 3
[0, 10]] # 150-200°C: 恒功率 10
)
print("温度控制策略:")
print(temp_control)
#温度控制策略:
#Piecewise((0.02*x, (x >= 0) & (x < 50)), (0.1*x - 2, (x >= 50) & (x < 100)),
(0.05*x + 3, (x >= 100) & (x < 150)), (10, (x >= 150) & (x <= 200)), (0, True))
示例 3: 税收计算分段函数（经济学）
tax_brackets = mkpp(
[0, 50000, 100000, 200000, 500000],
[[0.1, 0], # 0-50k: 10%
[0.2, -5000], # 50-100k: 20% - 5000
[0.3, -15000], # 100-200k: 30% - 15000
[0.4, -35000]] # 200-500k: 40% - 35000
)
print("累进税计算函数:")
print(tax_brackets)
#累进税计算函数:
#Piecewise((0.1*x, (x >= 0) & (x < 50000)), (0.2*x - 5000, (x >= 50000) & (x < 100000)),
(0.3*x - 15000, (x >= 100000) & (x < 200000)), (0.4*x - 35000, (x >= 200000) & (x <= 500000)), (0, True))
示例 4: 材料应力-应变关系（材料科学）
stress_strain = mkpp(
[0, 0.002, 0.01, 0.02],
[[200000, 0], # 弹性阶段: 200000ε
[50000, 300], # 屈服阶段: 50000ε + 300
[0, 800]] # 塑性阶段: 恒定应力800
)
print("材料应力-应变关系:")
print(stress_strain)
#材料应力-应变关系:
#Piecewise((200000*x, (x < 0.002) & (x >= 0)), (50000*x + 300, (x >= 0.002) & (x < 0.01)),
(800, (x >= 0.01) & (x <= 0.02)), (0, True))
示例 5: 信号处理分段滤波器（电子工程）
filter_response = mkpp(
[0, 1000, 5000, 10000, 20000],
[[0, 1], # 0-1kHz: 恒定增益1
[-0.0002, 1.2], # 1-5kHz: 线性衰减
[0, 0.2], # 5-10kHz: 恒定增益0.2
[-0.00001, 0.3]] # 10-20kHz: 缓慢衰减
)
print("滤波器频率响应:")
print(filter_response)
#滤波器频率响应:
#Piecewise((1, (x >= 0) & (x < 1000)),
(1.2 - 0.0002*x, (x >= 1000) & (x < 5000)),
(0.2, (x >= 5000) & (x < 10000)),
(0.3 - 1.0e-5*x, (x >= 10000) & (x <= 20000)), (0, True))
示例 6: 药物剂量-效应关系（药理学）
dose_response = mkpp(
[0, 10, 50, 100, 200],
[[0.05, 0], # 0-10mg: 线性增加
[0.01, 0.4], # 10-50mg: 缓慢增加
[0, 0.9], # 50-100mg: 平台期
[-0.002, 1.1]] # 100-200mg: 毒性下降
)
print("药物剂量-效应关系:")
print(dose_response)
#药物剂量-效应关系:
#Piecewise((0.05*x, (x >= 0) & (x < 10)),
(0.01*x + 0.4, (x >= 10) & (x < 50)),
(0.9, (x >= 50) & (x < 100)),
(1.1 - 0.002*x, (x >= 100) & (x <= 200)), (0, True))
示例 7: 机器人速度规划（自动化）
velocity_profile = mkpp(
([0, 2, 5, 8, 10],
[[1, 0], # 0-2s: 匀加速 a=1
[0, 2], # 2-5s: 匀速 v=2
[-1, 11], # 5-8s: 匀减速
[0, 0]]) # 8-10s: 停止
)
print("机器人速度规划:")
print(velocity_profile)
#机器人速度规划:
#Piecewise((x, (x >= 0) & (x < 2)),
(2, (x >= 2) & (x < 5)),
(11 - x, (x >= 5) & (x < 8)), (0, True))
示例 8: 经济周期分析（宏观经济学）
business_cycle = mkpp(
[0, 2, 4, 6, 8],
[[0.5, 2], # 复苏期: 0.5t + 2
[1, 1], # 扩张期: t + 1
[-0.5, 7], # 衰退期: -0.5t + 7
[0, 4]]) # 萧条期: 恒定4
)
print("经济周期GDP增长模型:")
print(business_cycle)
#经济周期GDP增长模型:
#Piecewise((0.5*x + 2, (x >= 0) & (x < 2)),
(x + 1, (x >= 2) & (x < 4)),
(7 - 0.5*x, (x >= 4) & (x < 6)),
(4, (x >= 6) & (x <= 8)), (0, True))
示例 9: 光学透镜像差校正（光学工程）
lens_correction = mkpp(
[0, 0.3, 0.7, 1.0],
[[-2, 0], # 中心区域: -2r^2
[1, -1.18], # 中间区域: r^2 - 1.18
[3, -3.18]] # 边缘区域: 3r^2 - 3.18
)
print("透镜像差校正函数:")
print(lens_correction)
#透镜像差校正函数:
#Piecewise((-2*x, (x < 0.3) & (x >= 0)),
(x - 1.18, (x >= 0.3) & (x < 0.7)),
(3*x - 3.18, (x >= 0.7) & (x <= 1.0)), (0, True))
示例 10: 神经网络激活函数（人工智能）
custom_activation = mkpp(
[-5, -1, 0, 1, 5],
[[0, -0.1], # x<-1: 恒定-0.1
[0.3, 0.2], # -1~0: 0.3x + 0.2
[0.7, 0.2], # 0~1: 0.7x + 0.2
[0, 0.9]]) # x>1: 恒定0.9
)
print("自定义神经网络激活函数:")
print(custom_activation)
#自定义神经网络激活函数:
#Piecewise((-0.1, (x >= -5) & (x < -1)),
(0.3*x + 0.2, (x >= -1) & (x < 0)),
(0.7*x + 0.2, (x >= 0) & (x < 1)),
(0.9, (x >= 1) & (x <= 5)), (0, True))
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def make_piece_poly(input_str):
"""
生成分段多项式（SymPy 符号版），修复节点外返回 NaN 的问题，添加默认值
参数格式示例:
"([-4, 0], [[0.25, -1, 0]])" # 二维结构需要明确括号层次
改进点:
1. 支持 SymPy 的 Tuple 结构解析
2. 添加默认区间外返回 0 的逻辑
3. 增强类型转换鲁棒性
"""
try:
x = sp.Symbol('x')
expr = sp.sympify(input_str, evaluate=False)
error_msg = None
# 解析输入结构
if not isinstance(expr, (tuple, sp.Tuple)) or len(expr) != 2:
raise ValueError("输入必须为 (x_nodes, coefs) 的元组形式")
# 转换节点为浮点数列表
'''
x_nodes = [float(node) for node in expr[0]]
# 转换系数为二维浮点数列表
coefs = []
for row in expr[1]:
if not isinstance(row, (list, tuple, sp.Tuple)):
raise ValueError("系数矩阵每行必须为元组")
coefs.append([float(c) for c in row])
'''
x_nodes, coefs = expr
# 验证节点和系数维度
n_nodes = len(x_nodes)
if n_nodes < 2:
raise ValueError("至少需要2个节点")
if len(coefs) != n_nodes - 1:
raise ValueError("系数矩阵行数必须等于节点数减一")
# 构造分段多项式
pieces = []
for i in range(n_nodes - 1):
a, b = x_nodes[i], x_nodes[i + 1]
degree = len(coefs[i]) - 1
# 构建多项式: c0*x^d + c1*x^(d-1) + ... + cd
poly = sum(coefs[i][k] * x ** (degree - k) for k in range(degree + 1))
# 设置区间条件（最后一个区间闭右端）
if i < n_nodes - 2:
cond = sp.And(x >= a, x < b)
else:
cond = sp.And(x >= a, x <= b)
pieces.append((poly, cond))
# 添加默认区间外返回 0
pieces.append((0, True))
return sp.Piecewise(*pieces)
except Exception as e:
return f"错误：{str(e)}"
# 测试示例
if __name__ == "__main__":
# 示例 1: 正常情况
pp = make_piece_poly("([-4, 0], [[1/4, -1, 0]])")
if isinstance(pp, sp.Piecewise):
print("分段多项式:")
print(pp)
# Piecewise((x**2/4 - x, (x >= -4) & (x <= 0)), (0, True))
print("x=-2:", pp.subs(sp.Symbol('x'), -2))
# 3
print("x=1:", pp.subs(sp.Symbol('x'), 1))
# 0
else:
print(pp)
最大似然估计参数
phat=mle(data)使用样本数据返回正态分布参数的最大似然估计（mle）
[phat，pci]=mle（___）还使用前面语法中的任何输入参数组合返回参数的置信区间。
data —— 样本数据和审查信息,向量,二维矩阵
示例 1: 产品质量检验（制造业）
检测100个产品中的缺陷品数量
quality_data = [0,1,0,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0]
result = mle(quality_data, binomial)
print(f"产品缺陷率估计: {result}")
#产品缺陷率估计: (0.14, [0.0438, 0.2362])
示例 2: 不同生产批次的合格率（质量管理）
每行表示一个批次：[合格数量, 总产量]
batch_data = [[95,100], [98,100], [92,100], [96,100], [97,100], [94,100], [99,100], [93,100], [96,100], [95,100]]
result = mle(batch_data, binomial)
print(f"整体合格率估计: {result}")
#整体合格率估计: (0.955, [0.9422, 0.9678])
示例 3: 临床试验药物有效性（医药研发）
不同剂量组的疗效数据：[有效病例数, 总病例数]
clinical_trial = [[8,20], [12,20], [15,20], [18,20]]
result = mle(clinical_trial, binomial)
print(f"综合有效率估计: {result}")
#综合有效率估计: (0.6625, [0.5589, 0.7661])
示例 4: 客户转化率分析（市场营销）
不同广告渠道的转化数据：[转化用户数, 曝光用户数]
conversion_data = [[150,10000], [80,5000], [200,15000], [120,8000]]
result = mle(conversion_data, binomial)
print(f"平均转化率估计: {result}")
#平均转化率估计: (0.0145, [0.0133, 0.0157])
示例 5: 网络设备故障率（IT运维）
每月故障设备数量：[故障设备数, 总设备数]
failure_data = [[2,100], [3,100], [1,100], [4,100], [2,100], [3,100], [1,100], [2,100], [3,100], [2,100], [1,100], [4,100]]
result = mle(failure_data, binomial)
print(f"月故障率估计: {result}")
#月故障率估计: (0.0233, [0.0148, 0.0319])
示例 6: 泊松分布 - 呼叫中心来电数量（服务业）
每小时来电数量
call_data = [12,15,8,10,14,9,11,13,10,12,16,7,14,11,9,13,10,15,8,12,14,11,10,13]
result = mle(call_data, poisson)
print(f"平均每小时来电数: {result}")
#平均每小时来电数: (11.5417, [10.1825, 12.9008])
示例 7: 泊松分布 - 网站访问量（互联网）
每分钟页面访问量
web_traffic = [45,38,52,41,47,39,50,43,48,40,51,42,46,37,53,44,49,41,47,38,52,43,46,39]
result = mle(web_traffic, poisson)
print(f"平均每分钟访问量: {result}")
#平均每分钟访问量: (44.625, [41.9524, 47.2976])
示例 8: A/B测试成功率（数据科学）
A组和B组的转化数据
group_a = [1,0,1,1,0,1,0,1,1,1,0,1,1,0,1,1,1,0,1,1,0,1,1,1,0,1,1,1,1,0] # 30次试验
group_b = [1,0,0,1,0,1,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0] # 30次试验
result_a = mle(group_a, binomial)
result_b = mle(group_b, binomial)
print(f"A组转化率: {result_a}")
print(f"B组转化率: {result_b}")
#A组转化率: (0.7, [0.536, 0.864])
#B组转化率: (0.4, [0.2247, 0.5753])
示例 9: 流行病发病率估计（公共卫生）
不同地区发病率：[发病人数, 人口数]
disease_incidence = [[25,10000], [18,8000], [32,12000], [15,7500], [28,11000]]
result = mle(disease_incidence, binomial)
print(f"疾病发病率估计: {result}")
#疾病发病率估计: (0.0024, [0.002, 0.0029])
示例 10: 信用违约概率（金融风控）
不同信用等级的违约数据：[违约客户数, 总客户数]
default_data = [[5,500], [12,300], [25,200], [45,150], [80,100]]
result = mle(default_data, binomial)
print(f"平均违约概率: {result}")
#平均违约概率: (0.1336, [0.1147, 0.1525])
示例 11: 指数分布 - 设备寿命测试（可靠性工程）
设备故障时间（小时）
lifetime_data = [1200,800,1500,600,2000,900,1100,1700,1300,750,1800,950,1400,850,1600,700,1900,1000,1250,650]
result = mle(lifetime_data, exponential)
print(f"设备寿命参数估计: {result}")
#设备寿命参数估计: (0.0008, [0.0005, 0.0012])
示例 12: 指数分布 - 客户服务时间（运营管理）
每次客户服务所需时间（分钟）
service_time = [8.5,12.3,6.7,15.2,9.8,11.5,7.9,13.6,10.4,14.1,8.2,12.8,7.5,16.3,9.1,11.9,8.8,13.2,7.1,15.7]
result = mle(service_time, exponential)
print(f"服务时间参数估计: {result}")
#服务时间参数估计: (0.0907, [0.0554, 0.1345])
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.stats import binom, norm, chi2
from functools import partial
from scipy.optimize import minimize
def maximum_likelihood_estimate(input_str):
"""
改进的最大似然估计函数，支持多种分布
参数:
input_str: 字符串形式的输入数据
distribution: 分布类型 ('binomial', 'poisson', 'exponential')
confidence_level: 置信水平 (默认0.95)
返回:
估计参数和置信区间
"""
try:
# SymPy 表达式解析
expr = sp.sympify(input_str)
data = None
def _binomial_mle(data, confidence_level):
"""二项分布的最大似然估计"""
# 列数校验与处理
if data.shape[1] == 1: # 向量输入
n_col = np.ones((data.shape[0], 1))
data = np.hstack([data, n_col])
elif data.shape[1] != 2:
return "输入矩阵列数必须为1（向量）或2（k和n）"
k = data[:, 0]
n_values = data[:, 1]
# 数据校验
if np.any(k < 0) or np.any(k > n_values) or np.any(n_values <= 0):
return "数据无效：成功次数应在[0, n]范围内且n>0"
# 最大似然估计（解析解）
total_successes = np.sum(k)
total_trials = np.sum(n_values)
p_hat = total_successes / total_trials
# 计算置信区间
z = norm.ppf(1 - (1 - confidence_level) / 2)
se = np.sqrt(p_hat * (1 - p_hat) / total_trials)
ci = (max(0, p_hat - z * se), min(1, p_hat + z * se))
return {
'distribution': 'binomial',
'parameter': 'p',
'estimate': round(p_hat, 4),
'confidence_interval': [round(ci[0], 4), round(ci[1], 4)],
'total_successes': int(total_successes),
'total_trials': int(total_trials)
}
def _poisson_mle(data, confidence_level):
"""泊松分布的最大似然估计"""
if data.shape[1] != 1:
return "泊松分布需要一维数据向量"
observations = data.flatten()
if np.any(observations < 0):
return "泊松分布数据不能为负"
# 最大似然估计（解析解）
lambda_hat = np.mean(observations)
# 计算置信区间
z = norm.ppf(1 - (1 - confidence_level) / 2)
se = np.sqrt(lambda_hat / len(observations))
ci = (max(0, lambda_hat - z * se), lambda_hat + z * se)
return {
'distribution': 'poisson',
'parameter': 'lambda',
'estimate': round(lambda_hat, 4),
'confidence_interval': [round(ci[0], 4), round(ci[1], 4)],
'sample_size': len(observations)
}
def _exponential_mle(data, confidence_level):
"""指数分布的最大似然估计"""
if data.shape[1] != 1:
return "指数分布需要一维数据向量"
observations = data.flatten()
if np.any(observations <= 0):
return "指数分布数据必须为正数"
# 最大似然估计（解析解）
lambda_hat = 1 / np.mean(observations)
# 计算置信区间（使用卡方分布）
n = len(observations)
chi2_lower = chi2.ppf((1 - confidence_level) / 2, 2 * n)
chi2_upper = chi2.ppf(1 - (1 - confidence_level) / 2, 2 * n)
ci = (chi2_lower / (2 * n * np.mean(observations)),
chi2_upper / (2 * n * np.mean(observations)))
return {
'distribution': 'exponential',
'parameter': 'lambda',
'estimate': round(lambda_hat, 4),
'confidence_interval': [round(ci[0], 4), round(ci[1], 4)],
'sample_size': n,
'mean_lifetime': round(1 / lambda_hat, 4)
}
distribution = 'binomial',
confidence_level = 0.95
if isinstance(expr, tuple):
if len(expr) == 3:
data_expr, distribution, confidence_level = expr
elif len(expr) == 2:
data_expr, distribution = expr
else:
raise ValueError(f"处理过程中发生错误: {input_str}")
else:
data_expr = expr
if isinstance(data_expr,list):
data = np.array(data_expr, dtype=float)
else:
raise ValueError(f"处理过程中发生错误: {input_str}")
distribution = str(distribution)
# 数据预处理
if data.ndim == 1:
data = data.reshape(-1, 1)
if distribution == 'binomial':
return _binomial_mle(data, confidence_level)
elif distribution == 'poisson':
return _poisson_mle(data, confidence_level)
elif distribution == 'exponential':
return _exponential_mle(data, confidence_level)
else:
return f"不支持的分布类型: {distribution}"
except Exception as e:
return f"处理过程中发生错误: {str(e)}"
# 测试用例
if __name__ == "__main__":
# 示例 1: 产品质量检验（制造业）
print("=== 制造业：产品质量检验 ===")
# 检测100个产品中的缺陷品数量
quality_data = "[0,1,0,0,1,0,0,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0]"
result = maximum_likelihood_estimate(f"{quality_data}, binomial")
print(f"产品缺陷率估计: {result}")
# 示例 2: 不同生产批次的合格率（质量管理）
print("\n=== 质量管理：多批次合格率 ===")
# 每行表示一个批次：[合格数量, 总产量]
batch_data = "[[95,100], [98,100], [92,100], [96,100], [97,100], [94,100], [99,100], [93,100], [96,100], [95,100]]"
result = maximum_likelihood_estimate(f"{batch_data}, binomial")
print(f"整体合格率估计: {result}")
# 示例 3: 临床试验药物有效性（医药研发）
print("\n=== 医药研发：药物有效性试验 ===")
# 不同剂量组的疗效数据：[有效病例数, 总病例数]
clinical_trial = "[[8,20], [12,20], [15,20], [18,20]]"
result = maximum_likelihood_estimate(f"{clinical_trial}, binomial")
print(f"综合有效率估计: {result}")
多项式概率密度函数
Y = mnpdf(X,PROB)返回具有概率PROB的多项式分布的pdf,在X的每一行进行评估.X和PROB是m-by-k矩阵或1-by-k向量,其中k是多项式箱或类别的数量.
PROB的每一行必须加起来为1,每个观测值（X行）的样本量由行和总和（X,2）给出. Y是一个m乘1的向量,mnpdf使用输入的相应行计算Y的每一行,或者在需要时复制它们.
示例1：市场调研 - 消费者偏好调查
在100名受访者中，35人喜欢产品A，40人喜欢产品B，25人喜欢产品C
result1 = mnpdf([35, 40, 25], [0.3, 0.4, 0.3])
print(f"输出：{result1}")
print(f"解释：观察到这种偏好分布的概率约为 {float(result1):.6f}\n")
#输出：0.00365725460712191
#解释：观察到这种偏好分布的概率约为 0.003657
示例2：质量控制 - 产品缺陷分类
在50个样本中，2个有外观缺陷，3个有功能缺陷，45个合格
result2 = mnpdf([2, 3, 45], [0.05, 0.08, 0.87])
print(f"输出：{result2}")
print(f"解释：观察到这种缺陷分布的概率约为 {float(result2):.6f}\n")
#输出：0.0514823191784355
#解释：观察到这种缺陷分布的概率约为 0.051482
示例3：医学研究 - 疾病症状分布
在80名患者中，30人有发热，25人有咳嗽，25人有其他症状
result3 = mnpdf([30, 25, 25], [0.4, 0.3, 0.3])
print(f"输出：{result3}")
print(f"解释：观察到这种症状分布的概率约为 {float(result3):.6f}\n")
#输出：0.00928194019203596
#解释：观察到这种症状分布的概率约为 0.009282
示例4：符号计算 - 理论概率分析
分析三种结果的理论概率分布
result4 = mnpdf([x, y, z], [0.2, 0.5, 0.3])
print(f"输出：{result4}")
print("解释：得到了包含符号变量的概率表达式\n")
#输出：Piecewise((9.33262154439442e+157*0.2**x*0.3**z*0.5**y/(factorial(x)*factorial(y)*factorial(z)),
Eq(x + y + z, 100)), (0, True))
#解释：得到了包含符号变量的概率表达式
示例5：选举预测 - 选民投票倾向
在200名选民中，85人支持候选人A，75人支持候选人B，40人未决定
result5 = mnpdf([85, 75, 40], [0.45, 0.4, 0.15])
print(f"输出：{result5}")
print(f"解释：观察到这种投票分布的概率约为 {float(result5):.8f}\n")
#输出：0.000727850747946614
#解释：观察到这种投票分布的概率约为 0.00072785
示例6：教育评估 - 学生成绩分布
在60名学生中，15人得A，20人得B，15人得C，10人得D
result7 = mnpdf([15, 20, 15, 10], [0.2, 0.35, 0.25, 0.2])
print(f"输出：{result7}")
print(f"解释：观察到这种成绩分布的概率约为 {float(result7):.8f}\n")
#输出：0.00131068668057065
#解释：观察到这种成绩分布的概率约为 0.00131069
示例7：符号计算 - 参数化分析
使用多个符号变量进行参数化概率分析
result8 = mnpdf([a, b, c], [0.25, 0.35, 0.4])
print(f"输出：{result8}")
print("解释：得到了多符号变量的通用概率表达式\n")
#输出：Piecewise((9.33262154439442e+157*0.25**a*0.35**b*0.4**c/(factorial(a)*factorial(b)*factorial(c)),
Eq(a + b + c, 100)), (0, True))
#解释：得到了多符号变量的通用概率表达式
示例8：生物学 - 基因型分布
在植物杂交实验中，观察到的基因型分布
result9 = mnpdf([22, 47, 31], [0.25, 0.5, 0.25])
print(f"输出：{result9}")
print(f"解释：观察到这种基因型分布的概率约为 {float(result9):.10f}\n")
#输出：0.00341946691595411
#解释：观察到这种基因型分布的概率约为 0.0034194669
示例9：风险管理 - 事故类型统计
在工厂安全记录中，不同类型事故的分布
result10 = mnpdf([5, 8, 12, 25], [0.1, 0.15, 0.25, 0.5])
print(f"输出：{result10}")
print(f"解释：观察到这种事故分布的概率约为 {float(result10):.12f}\n")
#输出：0.00385168484670193
#解释：观察到这种事故分布的概率约为 0.003851684847
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import density, Multinomial, sample, MultivariateNormal
def multinomial_probability_density(input_str):
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
# 解析次数列表和概率列表
counts, probs = expr
# 检查输入有效性
if (len(counts) == len(probs)) and (sum(probs) <= 1):
# 处理符号情况
if any(sp.Symbol in type(e).__mro__ for e in counts):
n = 100 # 符号情况默认总次数
M = Multinomial('M', n, probs)
result = density(M)(*counts)
result = result.evalf()
else:
n = sum(counts) # 数值情况计算总次数
M = Multinomial('M', n, probs)
result = density(M)(*counts).evalf()
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
# 示例1：数值计算
print("示例1：数值计算")
input_str1 = "([2, 3, 5], [0.2, 0.3, 0.5])"
print(f"输入：{input_str1}")
print(f"输出：{multinomial_probability_density(input_str1)}\n")
# 输出：0.0850500000000000
# 示例2：符号计算
print("示例2：符号计算（需要预先定义符号）")
x, y, z = sp.symbols('x y z', integer=True, nonnegative=True)
input_str2 = "([x, y, z], [0.3, 0.3, 0.4])"
print(f"输入：{input_str2}")
print(f"输出：{multinomial_probability_density(input_str2)}\n")
# Piecewise((9.33262154439442e+157*0.3**x*0.3**y*0.4**z/(factorial(x)*factorial(y)*factorial(z)), Eq(x + y + z, 100)), (0, True))
多项式随机数
r=mnrnd(n,p)从具有参数n和p的多项式分布中返回随机值r. n是一个正整数,指定每个多项式结果的试验次数(样本量). p是多项式概率的1-byk向量, 其中k是多项式区间或类别的数量.
p的总和必须为1.(如果p不等于1,则r完全由NaN值组成.)r是一个1-by-k向量,包含k个多项式区间中每个区间的计数.
示例1：市场调研 - 消费者偏好模拟
模拟1000名消费者在三种产品中的选择分布
result1 = mnrnd(1000, [0.3, 0.4, 0.3])
print(f"输出：\n{result1}")
print(f"解释：模拟显示约{result1[0]}人选择产品A，{result1[1]}人选择产品B，{result1[2]}人选择产品C\n")
#[[276],
[404],
[320]]
#解释：模拟显示约276人选择产品A，404人选择产品B，320人选择产品C
示例2：选举预测 - 选民投票模拟
模拟5000名选民在三个候选人之间的投票分布
result2 = mnrnd(5000, [0.45, 0.35, 0.2])
print(f"输出：\n{result2}")
print(f"解释：候选人A得票约{result2[0]}，候选人B得票约{result2[1]}，候选人C得票约{result2[2]}\n")
#[[2263],
[1725],
[1012]]
#解释：候选人A得票约2263，候选人B得票约1725，候选人C得票约1012
示例3：质量控制 - 产品缺陷类型模拟
模拟10000个产品中不同类型缺陷的分布
result3 = mnrnd(10000, [0.02, 0.03, 0.95])
print(f"输出：\n{result3}")
print(f"解释：预计有{result3[0]}个外观缺陷，{result3[1]}个功能缺陷，{result3[2]}个合格产品\n")
#[[205],
[284],
[9511]]
#解释：预计有205个外观缺陷，284个功能缺陷，9511个合格产品
示例4：医学研究 - 疾病症状分布模拟
模拟200名患者中不同症状组合的分布
result4 = mnrnd(200, [0.25, 0.35, 0.2, 0.2])
print(f"输出：\n{result4}")
print(f"解释：症状A约{result4[0]}人，症状B约{result4[1]}人，症状C约{result4[2]}人，症状D约{result4[3]}人\n")
#[[49],
[62],
[46],
[43]]
#解释：症状A约49人，症状B约62人，症状C约46人，症状D约43人
示例5：生物学 - 基因型分布模拟
模拟孟德尔遗传实验中的基因型分布
result8 = mnrnd(1000, [0.25, 0.5, 0.25])
print(f"输出：\n{result8}")
print(f"解释：显性纯合子约{result8[0]}，杂合子约{result8[1]}，隐性纯合子约{result8[2]}\n")
#[[263],
[497],
[240]]
#解释：显性纯合子约263，杂合子约497，隐性纯合子约240
示例6：风险管理 - 事故类型统计模拟
模拟工厂不同类型安全事故的年度分布
result9 = mnrnd(365, [0.05, 0.1, 0.15, 0.7])
print(f"输出：\n{result9}")
print(f"解释：机械事故{result9[0]}天，电气事故{result9[1]}天，人为失误{result9[2]}天，安全无事故{result9[3]}天\n")
#[[20],
[30],
[62],
[253]]
#解释：机械事故20天，电气事故30天，人为失误62天，安全无事故253天
示例7：产品开发 - 功能需求优先级模拟
模拟用户对产品不同功能的需求优先级分布
result10 = mnrnd(500, [0.3, 0.25, 0.2, 0.15, 0.1])
print(f"输出：\n{result10}")
print(f"解释：核心功能{result10[0]}票，用户体验{result10[1]}票，性能优化{result10[2]}票，安全性{result10[3]}票，其他{result10[4]}票\n")
#[[167],
[106],
[109],
[80],
[38]]
#解释：核心功能167票，用户体验106票，性能优化109票，安全性80票，其他38票
示例8：大规模样本模拟
模拟大型电商平台的用户行为分布
result11 = mnrnd(100000, [0.6, 0.25, 0.1, 0.05])
print(f"输出样本数：{result11.rows} × {result11.cols}")
print(f"总样本量：{sum(result11)}")
print("解释：可用于大数据分析和机器学习模型训练\n")
#输出样本数：4 × 1
#总样本量：100000
#解释：可用于大数据分析和机器学习模型训练
示例9：小样本精确模拟
模拟小型实验的精确结果分布
result12 = mnrnd(10, [0.5, 0.3, 0.2])
print(f"输出：\n{result12}")
print(f"解释：小样本实验的精确随机结果，总和为{sum(result12)}\n")
#[[6],
[2],
[2]]
#解释：小样本实验的精确随机结果，总和为10
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import Multinomial, sample
def multinomial_random_numbers(input_str):
"""
生成多项分布随机数（对标 MATLAB 的 mnrnd 函数）
参数:
input_str (str): 输入表达式，格式为 (n, p)
- n: 试验次数（正整数）
- p: 概率向量（和为1）或概率矩阵（每行和为1）
返回:
list/sp.Matrix/str: 随机数结果或错误信息
示例:
>>> multinomial_random_numbers("(5, [0.3,0.7])")
Matrix([[1, 4]]) # 表示5次试验中类别1出现1次，类别2出现4次
>>> multinomial_random_numbers("(10, [[0.2,0.8], [0.5,0.5]])")
Matrix([[3,7], [5,5]]) # 生成两行样本
"""
try:
expr = sp.sympify(input_str)
if not isinstance(expr, tuple) or len(expr) != 2:
return "输入格式错误：应为 (n, p)"
n_expr, p_expr = expr[0], expr[1]
# ----------------------------
# 步骤 1: 参数解析和校验
# ----------------------------
# 解析试验次数 n
if not n_expr.is_integer or n_expr <= 0:
return "试验次数 n 必须是正整数"
n = int(n_expr)
# 解析概率 p (支持向量和矩阵)
p = sp.Matrix(p_expr) if isinstance(p_expr, list) else None
if p is None:
return "概率 p 必须是向量或矩阵"
# 检查概率有效性
def check_probability(row):
prob_sum = sum(row)
if not sp.utilities.lambdify((), prob_sum)() == 1.0:
return False
return all(x >= 0 for x in row)
# ----------------------------
# 步骤 2: 分情况处理
# ----------------------------
# 情况1: p是单行向量 (shape=[1, k])
if p.rows == 1 or p.cols == 1:
# 转换为概率向量
prob_vec = [float(x) for x in p] if p.is_symbolic() is False else p
if not check_probability(prob_vec):
return "概率和必须为1且元素非负"
# 生成随机样本
X = Multinomial('X', n, prob_vec)
sample_vec = sample(X)
return sp.Matrix([sample_vec]) if p.rows == 1 else sp.Matrix(sample_vec)
# 情况2: p是多行矩阵 (shape=[m, k])
else:
# 检查每行概率有效性
for row in p.tolist():
if not check_probability(row):
return f"概率行 {row} 无效"
# 为每行生成独立样本
samples = []
for i, row in enumerate(p.tolist()):
X = Multinomial(f'X_{i}', n, row)
samples.append(sample(X))
return sp.Matrix(samples)
except Exception as e:
return f"错误: {str(e)}"
# 单行概率向量
print(multinomial_random_numbers("(5, [0.3,0.7])"))
# Matrix([[1],
# [4]])
# 多行概率矩阵
print(multinomial_random_numbers("(10, [[0.2,0.8], [0.5,0.5]])"))
# Matrix([[3, 7],
# [8, 2]])
print(multinomial_random_numbers("(1000,[0.2,0.3,0.5])"))
# Matrix([[218],
# [259],
# [523]])
除后的余数(取模运算)
b = mod(a,m) 返回a除以m后的余数,其中a是被除数,m是除数.此函数通常称为取模运算,表达式为b=a-m*floor(a/m).
mod函数遵从mod(a,0)返回a的约定
a — 被除数,标量,向量,矩阵
m — 除数,标量,向量,矩阵
示例1：时间计算 - 时钟运算
计算当前时间经过若干小时后的时间
result1 = mod(14 + 20, 24) # 下午2点再过20小时
print(f"输出：{result1}")
print(f"解释：14点 + 20小时 = 第二天10点 (14+20 mod 24 = {result1})\n")
#输出：10.0
#解释：14点 + 20小时 = 第二天10点 (14+20 mod 24 = 10.0)
示例2：日历计算 - 星期计算
计算100天后是星期几（假设今天是星期1）
result2 = mod(1 + 100, 7)
print(f"输出：{result2}")
print(f"解释：100天后是星期{result2} (1+100 mod 7 = {result2})\n")
#输出：3.0
#解释：100天后是星期3.0 (1+100 mod 7 = 3.0)
示例3：密码学 - 哈希函数简化
对大数进行取模运算，用于哈希表索引
result3 = mod(123456789, 1000)
print(f"输出：{result3}")
print(f"解释：大数123456789在哈希表中的位置索引为{result3}\n")
#输出：789.0
#解释：大数123456789在哈希表中的位置索引为789.0
示例4：计算机科学 - 内存地址映射
计算虚拟地址到物理页框的映射
addresses = [0x1234, 0x5678, 0x9ABC, 0xDEF0]
page_size = 4096 # 4KB页大小
result4 = mod(addresses, page_size)
print(f"输出：{result4}")
print(f"解释：地址在页内的偏移量分别为{result4}\n")
#输出：[ 564. 1656. 2748. 3824.]
#解释：地址在页内的偏移量分别为[ 564. 1656. 2748. 3824.]
示例5：数学问题 - 奇偶性判断
判断一组数字的奇偶性
numbers = [15, 22, 37, 48, 59]
result5 = mod(numbers, 2)
print(f"输出：{result5}")
print(f"解释：余数为1是奇数，0是偶数 → 奇偶性: {['偶' if x==0 else '奇' for x in result5]}\n")
#输出：[1. 0. 1. 0. 1.]
#解释：余数为1是奇数，0是偶数 → 奇偶性: ['奇', '偶', '奇', '偶', '奇']
示例6：工程应用 - 角度归一化
将角度归一化到0-360度范围内
angles = [450, -90, 720, 180]
result6 = mod(angles, 360)
print(f"输出：{result6}")
print(f"解释：角度归一化结果: 450°→{result6[0]}°, -90°→{result6[1]}°, 720°→{result6[2]}°, 180°→{result6[3]}°\n")
#输出：[ 90. 270. 0. 180.]
#解释：角度归一化结果: 450°→90.0°, -90°→270.0°, 720°→0.0°, 180°→180.0°
示例7：矩阵运算 - 图像处理中的像素值归一化
将像素值限制在0-255范围内
result7 = mod([[300, -50, 128], [255, 500, 0], [100, 200, 400]], 256)
print(f"输出：\n{result7}")
print("解释：超出范围的像素值通过取模回到有效范围内\n")
#输出：
#[[44,206,128]
[255,244,0]
[100,200,144]]
#解释：超出范围的像素值通过取模回到有效范围内
示例8：金融计算 - 利息周期计算
计算投资在不同计息周期下的剩余天数
investment_days = 237
periods = [30, 90, 180, 365]
result8 = mod(investment_days, periods)
print(f"输出：{result8}")
print(f"解释：237天投资在30/90/180/365天周期下的剩余天数分别为{result8}\n")
#输出：[27,57,57,237]
#解释：237天投资在30/90/180/365天周期下的剩余天数分别为[ 27. 57. 57. 237.]
示例9：符号计算 - 数学公式推导
使用符号变量进行模运算公式推导
result9 = mod(x + 2*y, z)
print(f"输出：{result9}")
print(f"解释：得到了符号表达式 {result9}\n")
#输出：Mod(x + 2*y, z)
#解释：得到了符号表达式 Mod(x + 2*y, z)
示例10：数据分片 - 分布式计算
根据键值哈希确定数据分片位置
keys = [12345, 67890, 23456, 78901, 34567]
num_shards = 4
result12 = mod(keys, num_shards)
print(f"输出：{result12}")
print(f"解释：数据键分配到分片 {result12}\n")
#输出：[1,2,0,1,3]
#解释：数据键分配到分片 [1. 2. 0. 1. 3.]
示例11：颜色计算 - RGB值循环
实现颜色值的循环渐变
color_values = [200, 250, 300, 50, 150]
result14 = mod(color_values, 256)
print(f"输出：{result14}")
print(f"解释：颜色值循环结果: {result14}\n")
#输出：[200. 250. 44. 50. 150.]
#解释：颜色值循环结果: [200. 250. 44. 50. 150.]
示例12：复杂矩阵运算
矩阵与标量的模运算
result15 = mod([[15, 22, 37], [48, 59, 64], [71, 83, 99]], 7)
print(f"输出：\n{result15}")
print("解释：矩阵每个元素对7取模的结果\n")
#输出：
#[[1,1,2]
[6,3,1]
[1,6,1]]
#解释：矩阵每个元素对7取模的结果
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def modulo_dividend_operation(input_str):
"""
执行取模运算 (对标 MATLAB 的 mod 函数)
参数:
input_str (str): 输入表达式，格式为 (a, m)
- a 和 m 可以是标量或矩阵
- 若均为矩阵，则需同形；若一个为标量，则广播到另一矩阵
返回:
sp.Matrix/float/str: 取模结果或错误信息
示例:
>>> modulo_dividend_operation("(7, 3)") # 7 mod 3
1
>>> modulo_dividend_operation("([5,6;7,8], 3)") # 矩阵每个元素 mod 3
Matrix([[2, 0], [1, 2]])
>>> modulo_dividend_operation("(10, [2,3;4,5])") # 10 mod 每个元素
Matrix([[0, 1], [2, 0]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
if all(e.is_number for e in expr):
params = tuple(float(e.evalf()) for e in expr)
result = np.mod(*params)
elif any(e.free_symbols for e in expr):
result = sp.Mod(*expr)
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}" # 对标matlab mod函数(除后的余数（取模运算）), 请检查代码，增加注释和示范代码。
# 标量运算
print(modulo_dividend_operation("(7, 3)"))
# 输出: 1.0
print(modulo_dividend_operation("(10.5, x)"))
# Mod(10.5, x)
# 符号运算
print(modulo_dividend_operation("(x, 3)"))
# Mod(x, 3)
M = movmean(A,k) 返回由局部 k 个数据点的均值组成的数组，其中每个均值是基于 A 的相邻元素的长度为 k 的滑动窗计算得出。当 k 为奇数时，窗以当前位置的元素为中心。当 k 为偶数时，窗以当前元素及其前一个元素为中心。当没有足够的元素填满窗时，窗将自动在端点处截断。当窗口被截断时，只根据窗口内的元素计算平均值。M 与 A 的大小相同。
如果 A 是向量，movmean 将沿向量 A 的长度运算。
如果 A 是矩阵，则 movmean 沿 A 的大小不等于 1 的第一个维度进行运算。
M = movmean(A,[kb kf]) 通过长度为 kb+kf+1 的窗口计算均值，其中包括当前位置的元素、前面的 kb 个元素和后面的 kf 个元素。
M = movmean(___,dim) 为上述任一语法指定 A 的运算维度。例如，如果 A 是矩阵，则 movmean(A,k,2) 沿 A 的列运算，计算每行的 k 元素移动均值。
M = movmean(___,nanflag) 指定包含还是省略 A 中的 NaN 值。例如，movmean(A,k,"omitnan") 在计算每个均值时会忽略 NaN 值。默认情况下，movmean 包括 NaN 值。
A — 输入数组, 向量 | 矩阵
k — 窗长度, 数值或持续时间标量
dim — 沿其运算的维度, 正整数标量
nanflag — 缺失值条件, "omitmissing" (默认) | "omitnan" | "includemissing" | "includenan"
示例1. 简单向量移动平均
result1 = movmean([1,2,3,4,5,6,7,8,9,10], 3)
print(f"结果: {result1}\n")
#结果: [1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 9.5]
示例2. 非对称窗口 [2,1] 表示向前看2个，向后看1个
result3 = movmean([1,2,3,4,5,6,7,8,9,10], [2,1])
print(f"结果: {result3}\n")
#结果: [1.5, 2.0, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.0]
示例3. 2D矩阵按行计算移动平均
result4 = movmean([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], 2, 1)
print(f"结果: {result4}\n")
#结果:
#[[2.5, 3.5, 4.5],
[5.5, 6.5, 7.5],
[8.5, 9.5, 10.5],
[10.0, 11.0, 12.0]]
示例4. 2D矩阵按列计算移动平均
result5 = movmean([[1,2,3],[4,5,6],[7,8,9],[10,11,12]], 2, 2)
print(f"结果: {result5}\n")
#结果: [[1.5, 2.5, 3.0],
[4.5, 5.5, 6.0],
[7.5, 8.5, 9.0],
[10.5, 11.5, 12.0]]
示例5. 股票价格平滑处理 (5日移动平均)
stock_prices = "[100, 102, 101, 105, 108, 107, 110, 112, 115, 113]"
result = movmean(stock_prices, 5)
print(f"5日移动平均: {result}\n")
#5日移动平均:
#[101.0, 102.0, 103.2, 104.6, 106.2, 108.4, 110.4, 111.4, 112.5, 113.333333333333]
示例6. 传感器数据滤波
sensor_data = [1.0, 1.2, 0.9, 1.1, 3.0, 1.0, 0.8, 1.3, 1.1, 1.0]
result = movmean(sensor_data, [1,1])
print(f"滤波后数据: {result}\n")
#滤波后数据: [1.1, 1.03333333333333, 1.06666666666667, 1.66666666666667, 1.7, 1.6, 1.03333333333333, 1.06666666666667, 1.13333333333333, 1.05]
示例7. 销售数据趋势分析 (7日移动平均)
sales = [120, 125, 118, 130, 140, 135, 128, 132, 145, 150, 148, 155, 160, 158]
result = movmean(sales, 7)
print(f"7日移动平均: {result}\n")
#7日移动平均: [123.25, 126.6, 128.0, 128.0, 129.714285714286, 132.571428571429, 137.142857142857, 139.714285714286,
141.857142857143, 145.428571428571, 149.714285714286, 152.666666666667, 154.2, 155.25]
示例8. 边界情况测试
窗口大小等于数据长度
result1 = movmean([1,2,3,4,5], 5)
print(f"结果: {result1}\n")
#结果: [2.0, 2.5, 3.0, 3.5, 4.0]
窗口大小为1
result2 = movmean([1,2,3,4,5], 1)
print(f"结果: {result2}\n")
#结果: [1.0, 2.0, 3.0, 4.0, 5.0]
单元素数组
result4 = movmean([5], 1)
print(f"结果: {result4}\n")
#结果: [5.0]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
def move_mean_symbolic(input_str):
"""
处理数值移动均值的入口函数，解析输入参数并调用核心计算函数。
参数:
input_expr: 输入的表达式，可以是矩阵、列表或包含矩阵及参数的元组
返回:
移动均值的结果或错误信息
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def movmean(A, k=None, dim=None, nanflag="includenan"):
# 如果 k 是列表形式 [kb kf]
if isinstance(k, list):
kb, kf = k
window_size = kb + kf + 1
elif k is not None:
window_size = k
else:
raise ValueError("k must be specified.")
A = np.asarray(A)
# 如果没有指定维度，根据 A 的类型确定维度
if dim is None:
if A.ndim == 1:
dim = 0
else:
dim = np.argmax(A.shape)
# 处理 NaN 值
if nanflag == "omitnan":
def window_mean(window):
valid_values = window[~np.isnan(window)]
if valid_values.size == 0:
return np.nan
return np.mean(valid_values)
else:
window_mean = np.mean
# 处理向量情况
if A.ndim == 1:
M = np.zeros_like(A, dtype=float)
for i in range(len(A)):
if isinstance(k, list):
start = max(0, i - kb)
end = min(len(A), i + kf + 1)
else:
if k % 2 == 1:
half_k = k // 2
start = max(0, i - half_k)
end = min(len(A), i + half_k + 1)
else:
half_k = k // 2
start = max(0, i - half_k + 1)
end = min(len(A), i + half_k + 1)
M[i] = window_mean(A[start:end])
return M
# 处理矩阵情况
elif A.ndim == 2:
M = np.zeros_like(A, dtype=float)
if dim == 0:
for j in range(A.shape[1]):
for i in range(A.shape[0]):
if isinstance(k, list):
start = max(0, i - kb)
end = min(A.shape[0], i + kf + 1)
else:
if k % 2 == 1:
half_k = k // 2
start = max(0, i - half_k)
end = min(A.shape[0], i + half_k + 1)
else:
half_k = k // 2
start = max(0, i - half_k + 1)
end = min(A.shape[0], i + half_k + 1)
M[i, j] = window_mean(A[start:end, j])
elif dim == 1:
for i in range(A.shape[0]):
for j in range(A.shape[1]):
if isinstance(k, list):
start = max(0, j - kb)
end = min(A.shape[1], j + kf + 1)
else:
if k % 2 == 1:
half_k = k // 2
start = max(0, j - half_k)
end = min(A.shape[1], j + half_k + 1)
else:
half_k = k // 2
start = max(0, j - half_k + 1)
end = min(A.shape[1], j + half_k + 1)
M[i, j] = window_mean(A[i, start:end])
return M
else:
raise ValueError("Input array A must be 1D or 2D.")
if isinstance(expr, tuple) and isinstance(expr[0], list):
A = np.array(expr[0], dtype=float)
if len(expr) > 3:
k = expr[1]
dim = int(expr[2])
nanflag = str(expr[3]).lower()
elif len(expr) == 3:
k = expr[1]
if expr[2].is_integer:
dim = int(expr[2])
nanflag = 'includemissing'
else:
dim = None
nanflag = str(expr[2]).lower()
elif len(expr) == 2:
k = expr[1]
dim = None
nanflag = 'includemissing'
if nanflag not in ['includenan', 'includemissing', 'omitmissing', 'omitnan']:
error = True
result = movmean(A=A, k=k, dim=dim, nanflag=nanflag)
elif isinstance(expr, list):
A = sp.Matrix(expr)
result = movmean(A=A)
else:
error = True
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
def main():
"""
主程序，用于测试 cumulative_max_symbolic 函数。
"""
# 定义测试用的输入字符串
test_inputs = [
"[4,8,nan,-1,-2,-3,nan,3,4,5],3,omitnan",
# Matrix([[6.00000000000000],
# [6.00000000000000],
# [3.50000000000000],
# [-1.50000000000000],
# [-2.00000000000000],
# [-2.50000000000000],
# [0],
# [3.50000000000000],
# [4.00000000000000],
# [4.50000000000000]])
"[4,8,6,-1,-2,-3,-1,3,4,5],3",
# Matrix([[6.00000000000000],
# [6.00000000000000],
# [4.33333333333333],
# [1.00000000000000],
# [-2.00000000000000],
# [-2.00000000000000],
# [-0.333333333333333],
# [2.00000000000000],
# [4.00000000000000],
# [4.50000000000000]])
"[4,8,6,-1,-2,-3,-1,3,4,5],[2,0]",
# Matrix([[4.00000000000000],
# [6.00000000000000],
# [6.00000000000000],
# [4.33333333333333],
# [1.00000000000000],
# [-2.00000000000000],
# [-2.00000000000000],
# [-0.333333333333333],
# [2.00000000000000],
# [4.00000000000000]])
]
# 遍历测试输入，调用 movsum 函数并打印结果
for input_str in test_inputs:
result = move_mean_symbolic(input_str)
print(f"输入: {input_str}")
print(f"结果: {result}")
print()
if __name__ == "__main__":
# 调用主程序
main()
M = movsum(A,k) 返回由局部 k 个数据点总和组成的数组，其中每个总和基于 A 的相邻元素的长度为 k 的滑动窗计算得出。当 k 为奇数时，窗以当前位置的元素为中心。当 k 为偶数时，窗以当前元素及其前一个元素为中心。当没有足够的元素填满窗时，窗将自动在端点处截断。当窗口被截断时，只根据窗口内的元素计算总和。M 与 A 的大小相同。
如果 A 是向量，movsum 将沿向量 A 的长度运算。
M = movsum(A,[kb kf]) 通过长度为 kb+kf+1 的窗口计算和，其中包括当前位置的元素、前面的 kb 个元素和后面的 kf 个元素。
M = movsum(___,dim) 为上述任一语法指定 A 的运算维度。例如，如果 A 是矩阵，则 movsum(A,k,2) 沿 A 的列运算，计算每行的 k 个元素的移动总和。
M = movsum(___,nanflag) 指定包含还是省略 A 中的 NaN 值。例如，movsum(A,k,"omitnan") 在计算每个和时忽略所有 NaN 值。默认情况下，movsum 包括 NaN 值。
A — 输入数组, 向量 | 矩阵
k — 窗长度, 数值或持续时间标量
dim — 沿其运算的维度, 正整数标量
nanflag — 缺失值条件, "omitmissing" (默认) | "omitnan" | "includemissing" | "includenan"
示例1. 简单向量移动总和
result1 = movsum([1,2,3,4,5,6,7,8,9,10], 3)
print(f"结果: {result1}\n")
#结果: [3.0, 6.0, 9.0, 12.0, 15.0, 18.0, 21.0, 24.0, 27.0, 19.0]
示例2. 非对称窗口 [2,1] 表示向前看2个，向后看1个
result3 = movsum([1,2,3,4,5,6,7,8,9,10], [2,1])
print(f"结果: {result3}\n")
#结果: [3.0, 6.0, 10.0, 14.0, 18.0, 22.0, 26.0, 30.0, 34.0, 27.0]
示例3. 销售额滚动总和(7日滚动总和)
daily_sales = [1000, 1200, 800, 1500, 2000, 1800, 1600, 1400, 1700, 1900, 2100, 2200]
result = movsum(daily_sales, 7)
print(f"7日滚动总和: {result}\n")
#7日滚动总和: [4500.0, 6500.0, 8300.0, 9900.0, 10300.0, 10800.0, 11900.0, 12500.0, 12700.0, 10900.0, 9300.0, 7900.0]
示例4. 库存累积计算
inventory_changes = [50, -20, 30, -15, 40, -25, 35, -10, 45, -30]
result = movsum(inventory_changes, [0,2])
print(f"未来3期累积变化: {result}\n")
#未来3期累积变化: [60.0, -5.0, 55.0, 0.0, 50.0, 0.0, 70.0, 5.0, 15.0, -30.0]
示例5. 电力消耗滚动计算 (24小时总和)
hourly_power = [100, 95, 90, 85, 80, 75, 120, 150, 180, 200, 190, 180, 170, 160, 150, 140, 130, 200, 220, 210, 190, 170, 150, 130]
result = movsum(hourly_power, 24)
print(f"24小时电力消耗滚动总和: {result}\n")
#24小时电力消耗滚动总和:
#[1715.0, 1875.0, 2025.0, 2165.0, 2295.0, 2495.0, 2715.0, 2925.0, 3115.0, 3285.0, 3435.0, 3565.0, 3465.0, 3370.0,
3280.0, 3195.0, 3115.0, 3040.0, 2920.0, 2770.0, 2590.0, 2390.0, 2200.0, 2020.0]
示例6. 网站访问量滚动统计
daily_visits = [1000, 1200, 800, 1500, 2000, 1800, 1600, 1400, 1700, 1900]
result = movsum(daily_visits, [6,0])
print(f"过去7天访问量总和: {result}\n")
#过去7天访问量总和: [1000.0, 2200.0, 3000.0, 4500.0, 6500.0, 8300.0, 9900.0, 10300.0, 10800.0, 11900.0]
示例7. 生产缺陷数量统计
defects = [0, 1, 0, 2, 0, 0, 1, 3, 0, 1, 0, 0, 2, 1, 0]
result = movsum(defects, 5)
print(f"5日缺陷滚动总和: {result}\n")
#5日缺陷滚动总和: [1.0, 3.0, 3.0, 3.0, 3.0, 6.0, 4.0, 5.0, 5.0, 4.0, 3.0, 4.0, 3.0, 3.0, 3.0]
示例8. 交易量滚动总和
volume = [10000, 15000, 8000, 12000, 20000, 18000, 22000, 19000, 16000, 14000]
result = movsum(volume, 5)
print(f"5日交易量滚动总和: {result}\n")
#5日交易量滚动总和: [33000.0, 45000.0, 65000.0, 73000.0, 80000.0, 91000.0, 95000.0, 89000.0, 71000.0, 49000.0]
示例9. 资金流入流出净额
cash_flow = [100, -50, 200, -100, 150, -80, 300, -120, 250, -90]
result = movsum(cash_flow, [0,4])
print(f"未来5日资金流预期总和: {result}\n")
#未来5日资金流预期总和: [300.0, 120.0, 470.0, 150.0, 500.0, 260.0, 340.0, 40.0, 160.0, -90.0]
示例10. 风险管理 - 损失累积
losses = [0, 0, 500, 0, 0, 1000, 0, 200, 0, 0, 300, 0, 1500]
result = movsum(losses, [0,11])
print(f"未来12个月损失预期总和: {result}\n")
#未来12个月损失预期总和:
#[2000.0, 3500.0, 3500.0, 3000.0, 3000.0, 3000.0, 2000.0, 2000.0, 1800.0, 1800.0, 1800.0, 1500.0, 1500.0]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import numpy as np
import sympy as sp
from itertools import product
def move_sum_symbolic(input_str):
"""
处理符号移动总和的入口函数，解析输入参数并调用核心计算函数。
参数:
input_expr: 输入的表达式，可以是矩阵、列表或包含矩阵及参数的元组
返回:
移动总和的结果或错误信息
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def movsum(A, k=None, dim=None, nanflag='includemissing'):
# 处理k参数
if isinstance(k, list) and len(k) == 2:
kb, kf = k
window_size = kb + kf + 1
elif k.is_integer:
if k <= 0:
raise ValueError("k must be a positive integer.")
window_size = k
kb = (k - 1) // 2
kf = (k - 1) - kb
else:
raise ValueError("k must be an integer or a list of two integers.")
A = np.asarray(A)
# 处理dim参数
if dim is None:
dim = next((i for i, size in enumerate(A.shape) if size > 1), 0)
else:
if dim < 0 or dim >= A.ndim:
raise ValueError(f"dim must be between 0 and {A.ndim - 1}")
# 处理NaN
if nanflag == 'omitnan' or 'omitmissing':
sum_func = np.nansum
elif nanflag == 'includenan' or 'includemissing':
sum_func = np.sum
else:
raise ValueError("nanflag must be 'includenan' or 'omitnan'.")
# 准备输出数组
result_shape = A.shape
M = np.empty_like(A)
# 生成索引组合
other_dims = [d for d in range(A.ndim) if d != dim]
other_indices = [range(A.shape[d]) for d in other_dims]
for indices in product(*other_indices):
# 构建当前向量的切片
slice_list = []
idx_iter = iter(indices)
for d in range(A.ndim):
if d == dim:
slice_list.append(slice(None))
else:
slice_list.append(next(idx_iter))
vector = A[tuple(slice_list)]
# 确保是一维数组
vector = np.atleast_1d(vector)
n = len(vector)
summed = np.zeros(n)
for i in range(n):
start = max(0, i - kb)
end = min(n, i + kf + 1)
summed[i] = sum_func(vector[start:end])
# 将结果放回M中
M[tuple(slice_list)] = summed
return M
if isinstance(expr, tuple) and isinstance(expr[0], list):
A = expr[0]
if len(expr) > 3:
k = expr[1]
dim = int(expr[2])
nanflag = str(expr[3]).lower()
elif len(expr) == 3:
k = expr[1]
if expr[2].is_integer:
dim = int(expr[2])
nanflag = 'includemissing'
else:
dim = None
nanflag = str(expr[2]).lower()
elif len(expr) == 2:
k = expr[1]
dim = None
nanflag = 'includemissing'
if nanflag not in ['includenan', 'includemissing', 'omitmissing', 'omitnan']:
error = True
result = movsum(A=A, k=k, dim=dim, nanflag=nanflag)
elif isinstance(expr, list):
A = sp.Matrix(expr)
result = movsum(A=A)
else:
error = True
return sp.Matrix(result) if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误:{e}"
def main():
"""
主程序，用于测试 cumulative_max_symbolic 函数。
"""
# 定义测试用的输入字符串
test_inputs = [
"[4,8,nan,-1,-2,-3,nan,3,4,5],3,omitnan",
# Matrix([[12.0000000000000],
# [nan],
# [nan],
# [nan],
# [-6.00000000000000],
# [nan],
# [nan],
# [nan],
# [12.0000000000000],
# [9.00000000000000]])
"[4,8,6,-1,-2,-3,-1,3,4,5],3",
# Matrix([[12.0000000000000],
# [18.0000000000000],
# [13.0000000000000],
# [3.00000000000000],
# [-6.00000000000000],
# [-6.00000000000000],
# [-1.00000000000000],
# [6.00000000000000],
# [12.0000000000000],
# [9.00000000000000]])
"[4,8,6,-1,-2,-3,-1,3,4,5],[2,0]",
# Matrix([[4.00000000000000],
# [12.0000000000000],
# [18.0000000000000],
# [13.0000000000000],
# [3.00000000000000],
# [-6.00000000000000],
# [-6.00000000000000],
# [-1.00000000000000],
# [6.00000000000000],
# [12.0000000000000]])
]
# 遍历测试输入，调用 movsum 函数并打印结果
for input_str in test_inputs:
result = move_sum_symbolic(input_str)
print(f"输入: {input_str}")
print(f"结果: {result}")
print()
if __name__ == "__main__":
# 调用主程序
main()
矩阵幂
C = A^B计算A的B次幂并将结果返回给C.
C = mpower(A,B)是执行A^B的替代方法.
A,B -- 操作数,标量,矩阵
基数A和指数B均为标量,在这种情况下A^B等效于A.^B。
基数A是方阵,指数B是标量.如果B为正整数,则按重复平方计算幂.对于B的其他值,计算使用特征值分解(对于大多数矩阵)或舒尔分解(对于亏损矩阵).
基数A是标量,指数B是方阵,计算使用特征值分解.
示例1. 基本矩阵幂运算
result1 = mpower([[1,2],[3,4]], 2)
print(f"结果: {result1}\n")
#结果:
#[[7, 10],
[15, 22]]
示例2. 单位矩阵的幂
result2 = mpower([[1,0],[0,1]], 5)
print(f"结果: {result2}\n")
#结果:
#[[1, 0],
[0, 1]]
示例3. 零矩阵的幂
result3 = mpower([[0,0],[0,0]], 3)
print(f"结果: {result3}\n")
#结果:
#[[0, 0],
[0, 0]]
示例4. 马尔可夫链 - 状态转移矩阵
天气转移矩阵: [晴, 雨]
今天晴->明天晴: 0.7, 晴->雨: 0.3
今天雨->明天晴: 0.4, 雨->雨: 0.6
P = [[0.7,0.3],[0.4,0.6]]
result = mpower(P, 3)
print(f"3步转移矩阵: {result}")
#3步转移矩阵:
#[[0.583, 0.417],
[0.556, 0.444]]
#说明: 矩阵的3次幂表示3天后的天气概率分布
示例5. 图论 - 邻接矩阵幂
4个节点的图邻接矩阵
A = [[0,1,1,0],[1,0,1,1],[1,1,0,0],[0,1,0,0]]
result = mpower(A, 2)
print(f"邻接矩阵的2次幂: {result}")
#邻接矩阵的2次幂:
#[[2, 1, 1, 1],
[1, 3, 1, 0],
[1, 1, 2, 1],
[1, 0, 1, 1]]
#说明: A^2[i,j] 表示从节点i到节点j长度为2的路径数量
示例6. 线性变换的复合
旋转矩阵
theta = 45 # 45度
cos_t = math.cos(math.radians(theta))
sin_t = math.sin(math.radians(theta))
R = [[cos_t,-sin_t],[sin_t,cos_t]]
result = mpower(R, 3)
print(f"旋转矩阵的3次幂 (135度旋转): {result}")
print("说明: 旋转矩阵的3次幂相当于连续旋转3次45度，即135度旋转\n")
#旋转矩阵的3次幂 (135度旋转):
#[[-0.706786486, -0.706786486],
[0.706786486, -0.706786486]]
#说明: 旋转矩阵的3次幂相当于连续旋转3次45度，即135度旋转
示例7. 标量的矩阵幂 (指数矩阵)
A = [[1,0],[0,2]]
result = mpower(2, A)
print(f"2^A (逐元素指数): {result}")
#2^A (逐元素指数):
#[[2, 0],
[0, 4]]
#说明: 对于对角矩阵，2^A相当于对每个对角线元素计算2的幂
示例8. 量子力学 - 时间演化算子
简单的2能级系统哈密顿量
H = [[1,-1],[-1,1]] # 哈密顿矩阵
result = mpower(H, 2)
print(f"哈密顿矩阵的2次幂: {result}")
#哈密顿矩阵的2次幂:
#[[2, -2],
[-2, 2]]
#说明: 在量子力学中，矩阵幂用于计算时间演化算子
示例9. 控制系统 - 状态转移矩阵
离散时间系统的状态矩阵
A = [[0.8,0.2],[0.1,0.9]]
result = mpower(A, 5)
print(f"5步状态转移矩阵: {result}")
#5步状态转移矩阵:
#[[0.44538, 0.55462],
[0.27731, 0.72269]]
#说明: 在控制系统中，A^k表示k步后的状态转移
示例10. 网络分析 - 连接性矩阵
社交网络连接矩阵
network = [[0,1,0,1],[1,0,1,0],[0,1,0,1],[1,0,1,0]]
result = mpower(network, 3)
print(f"3度连接矩阵: {result}")
#3度连接矩阵:
#[[0, 4, 0, 4],
[4, 0, 4, 0],
[0, 4, 0, 4],
[4, 0, 4, 0]]
#说明: 在网络分析中，A^k[i,j]表示通过k步连接的关系强度
示例11. 投资组合相关性矩阵
资产相关性矩阵
corr_matrix = [[1.0,0.6,0.3],[0.6,1.0,0.4],[0.3,0.4,1.0]]
result = mpower(corr_matrix, 2)
print(f"相关性矩阵的2次幂: {result}")
#相关性矩阵的2次幂:
#[[1.45, 1.32, 0.84],
[1.32, 1.52, 0.98],
[0.84, 0.98, 1.25]]
#说明: 在风险管理中，相关性矩阵的幂用于分析高阶相关性
示例12. 经济增长模型
经济部门间的投入产出系数矩阵
io_matrix = [[0.2,0.1,0.3],[0.1,0.3,0.2],[0.2,0.1,0.4]]
result = mpower(io_matrix, 3)
print(f"3轮效应矩阵: {result}")
#3轮效应矩阵:
#[[0.07, 0.055, 0.129],
[0.064, 0.062, 0.119],
[0.083, 0.064, 0.153]]
#说明: 在投入产出分析中，A^k表示k轮经济波及效应
示例13. 斐波那契数列的矩阵表示
fib_matrix = [[1,1],[1,0]]
result = mpower(fib_matrix, 5)
print(f"矩阵的5次幂: {result}")
print(f"斐波那契数列 F6 = {result[0, 1]}, F5 = {result[1, 1]}")
#矩阵的5次幂:
#[[8, 5],
[5, 3]]
#斐波那契数列 F6 = 5, F5 = 3
#说明: 通过矩阵幂可以快速计算斐波那契数列
示例14. 图像处理 - 变换组合
缩放矩阵
scale = [[2,0],[0,2]]
result = mpower(scale, 3)
print(f"缩放矩阵的3次幂: {result}")
#缩放矩阵的3次幂:
#[[8, 0],
[0, 8]]
#说明: 在图像处理中，变换矩阵的幂表示重复应用变换
示例15. 负整数幂 (需要可逆矩阵)
A = [[1,2],[3,4]]
result = mpower(A, -1)
print(f"结果: {result}")
#结果:
#[[-2, 1],
[3/2, -1/2]]
示例16. 零次幂
A = [[1,2],[3,4]]
result = mpower(A, 0)
print(f"结果: {result}")
#结果:
#[[1, 0],
[0, 1]]
#说明: 任何矩阵的0次幂都是单位矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import math
def matrix_power(input_str):
"""
执行矩阵幂运算（对标 MATLAB 的 mpower 函数）
参数:
input_str (str): 输入表达式，格式为 (A, k) 或 (k, A)
- A 必须为方阵（SymPy 矩阵）
- k 为标量（整数或浮点数）
返回:
sp.Matrix/str: 结果矩阵或错误信息
规则:
1. 若输入为 (A, k): 计算矩阵幂 A^k（A 必须为方阵，k 为整数）
2. 若输入为 (k, A): 计算标量幂 k^A = exp(ln(k) * A)（A 必须为方阵）
示例:
>>> matrix_power("([[1,2],[3,4]], 2)") # A^2
Matrix([[7, 10], [15, 22]])
>>> matrix_power("(2, [[1,0],[0,1]])") # 2^I = [[2,0],[0,2]]
Matrix([[2.0, 0], [0, 2.0]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
# Case 1: (k, A) => 计算 k^A = exp(ln(k) * A)
if isinstance(expr[1], list) and expr[0].is_number:
A = sp.Matrix(expr[1])
if A.is_square:
k = float(expr[0])
ln_k = sp.log(k)
result = (ln_k * A).exp()
else:
error = True
# Case 2: (A, k) => 计算 A^k
elif expr[1].is_integer and isinstance(expr[0], list):
A = sp.Matrix(expr[0])
if A.is_square:
k = expr[1]
result = A ** int(k)
else:
error = True
else:
error = True
else:
error = True
return result if not error else f"输入错误: {input_str}"
except Exception as e:
return f"错误: {str(e)}"
# 矩阵整数幂
print(matrix_power("[[1,2],[3,4]], 2"))
# Matrix([[7, 10],
# [15, 22]])
# 标量幂（k^A）
print(matrix_power("2,[[0,1],[1,0]]"))
# Matrix([[1.25000000000000, 0.750000000000000],
# [0.750000000000000, 1.25000000000000]])
求解关于x的线性方程组xA=B
x = B/A 对线性方程组x*A = B求解x. 矩阵A和B必须具有相同的列数.
如果A是标量,那么B/A等于B/A
如果A是n×n方阵,B是n列矩阵,那么x = B/A是方程x*A = B的解（如果存在解的话）
如果A是矩形m×n矩阵,且m ~= n,B是n列矩阵,那么x=B/A返回方程组x*A = B的最小二乘解
A, B — 操作数,向量,满矩阵,稀疏矩阵
x — 解,向量,满矩阵,稀疏矩阵
示例1. 精确解 - 方阵系统
result1 = mrdivide([[1,1,3],[2,0,4],[-1,6,-1]],[2,19,8])
print(f"方程: x * [[1,1,3],[2,0,4],[-1,6,-1]] = [2,19,8]")
print(f"解: {result1}\n")
#方程: x * [[1,1,3],[2,0,4],[-1,6,-1]] = [2,19,8]
#解: [[0.999999999999999],
[2.0],
[3.0]]
示例2. 最小二乘解 - 超定系统
result2 = mrdivide([[1,0],[2,0],[1,0]],[1,2])
print(f"方程: x * [[1,0],[2,0],[1,0]] = [1,2]")
print(f"最小二乘解: {result2}\n")
#方程: x * [[1,0],[2,0],[1,0]] = [1,2]
#最小二乘解: [[0.166666666666667],
[0.333333333333333],
[0.166666666666667]]
示例3. 电路分析 - 节点电压法
电路导纳矩阵 * 节点电压 = 电流源
方程: V * Y = I
Y_matrix = [[0.1,-0.05,0],[-0.05,0.15,-0.1],[0,-0.1,0.2]]
I_sources = [1,0,0.5]
result = mrdivide(Y_matrix,I_sources)
print(f"导纳矩阵 Y: {Y_matrix}")
print(f"电流源 I: {I_sources}")
print(f"节点电压 V: {result}")
#导纳矩阵 Y: [[0.1,-0.05,0],[-0.05,0.15,-0.1],[0,-0.1,0.2]]
#电流源 I: [1,0,0.5]
#节点电压 V: [[15],
[10],
[7.5]]
#说明: 求解电路节点电压
示例4. 结构力学 - 力平衡方程
刚度矩阵 * 位移 = 力
方程: u * K = F
K_matrix = [[100,-50,0],[-50,100,-50],[0,-50,100]]
F_forces = [10,0,5]
result = mrdivide(K_matrix,F_forces)
print(f"刚度矩阵 K: {K_matrix}")
print(f"外力 F: {F_forces}")
print(f"位移 u: {result}")
#刚度矩阵 K: [[100,-50,0],[-50,100,-50],[0,-50,100]]
#外力 F: [10,0,5]
#位移 u: [[0.175],
[0.15],
[0.125]]
#说明: 求解结构位移
示例5. 控制系统 - 状态反馈
希望的特征值配置
方程: K * A = desired_poles
A_system = [[0,1],[-2,-3]]
desired_poles = [-1,-2]
result = mrdivide(A_system,desired_poles)
print(f"系统矩阵 A: {A_system}")
print(f"期望极点: {desired_poles}")
print(f"反馈增益 K: {result}")
#系统矩阵 A: [[0,1],[-2,-3]]
#期望极点: [-1,-2]
#反馈增益 K: [[-0.5],
[0.5]]
#说明: 设计状态反馈控制器
示例6. 投入产出分析
技术系数矩阵 * 总产出 = 最终需求
方程: X * (I - A) = Y
A_tech = [[0.2,0.1],[0.3,0.2]]
Y_demand = [100,150]
构造 (I - A) 矩阵
result = mrdivide([[0.8,-0.1],[-0.3,0.8]],Y_demand)
print(f"技术系数矩阵 A: {A_tech}")
print(f"最终需求 Y: {Y_demand}")
print(f"总产出 X: {result}")
#技术系数矩阵 A: [[0.2,0.1],[0.3,0.2]]
#最终需求 Y: [100,150]
#总产出 X: [[204.918032786885],
[213.114754098361]]
#说明: 列昂惕夫投入产出模型
示例7. 投资组合优化
协方差矩阵 * 权重 = 期望收益
方程: w * Σ = μ
cov_matrix = [[0.04,0.02,0.01],[0.02,0.09,0.03],[0.01,0.03,0.16]]
expected_returns = [0.08,0.12,0.15]
result = mrdivide(cov_matrix,expected_returns)
print(f"协方差矩阵 Σ: {cov_matrix}")
print(f"期望收益 μ: {expected_returns}")
print(f"最优权重 w: {result}")
#协方差矩阵 Σ: [[0.04,0.02,0.01],
[0.02,0.09,0.03],
[0.01,0.03,0.16]]
#期望收益 μ: [0.08,0.12,0.15]
#最优权重 w: Matrix([[1.43423799582463], [0.780793319415449], [0.701461377870564]])
#说明: 马科维茨投资组合理论
示例8. 主成分分析
协方差矩阵 * 特征向量 = 特征值 * 特征向量
简化形式求解主要方向
cov_data = [[4,2],[2,3]]
target_direction = [1,0]
result = mrdivide(cov_data,target_direction)
print(f"数据协方差矩阵: {cov_data}")
print(f"目标方向: {target_direction}")
print(f"调整后的方向: {result}")
#数据协方差矩阵: [[4,2],[2,3]]
#目标方向: [1,0]
#调整后的方向: [[0.375],
[-0.25]]
#说明: PCA中求解特征方向
示例9. 滤波器设计
频率响应方程
方程: h * A = desired_response
A_freq = [[1,1,1],[1,0.5,0.25],[1,0,-1]]
desired_response = [1,0.7,0]
result = mrdivide(A_freq,desired_response)
print(f"频率点矩阵 A: {A_freq}")
print(f"期望响应: {desired_response}")
print(f"滤波器系数 h: {result}")
#频率点矩阵 A: [[1,1,1],[1,0.5,0.25],[1,0,-1]]
#期望响应: [1,0.7,0]
#滤波器系数 h: [[1.5],
[-1.6],
[1.1]]
#说明: FIR滤波器系数设计
示例10. 量子力学 - 期望值
算符矩阵 * 态向量 = 期望值
方程: ψ * H = E
H_operator = [[1,0.5],[0.5,1]]
E_energy = [1,1]
result = mrdivide(H_operator,E_energy)
print(f"哈密顿算符 H: {H_operator}")
print(f"能量期望 E: {E_energy}")
print(f"波函数 ψ: {result}")
#哈密顿算符 H: [[1,0.5],[0.5,1]]
#能量期望 E: [1,1]
#波函数 ψ: [[0.666666666666667],
[0.666666666666667]]
#说明: 求解量子态的期望值
示例11. 热传导 - 温度分布
热传导矩阵 * 温度 = 热源
方程: T * K = Q
K_thermal = [[2,-1,0],[-1,2,-1],[0,-1,2]]
Q_source = [10,0,5]
result = mrdivide(K_thermal,Q_source)
print(f"热传导矩阵 K: {K_thermal}")
print(f"热源 Q: {Q_source}")
print(f"温度分布 T: {result}")
#热传导矩阵 K: [[2,-1,0],[-1,2,-1],[0,-1,2]]
#热源 Q: [10,0,5]
#温度分布 T: [[8.75],
[7.5],
[6.25]]
#说明: 稳态热传导问题求解
示例12. 奇异矩阵 - 最小二乘解
result1 = mrdivide([[1,2,3],[2,4,6],[3,6,9]],[1,2,3])
print(f"奇异矩阵方程求解: {result1}\n")
#奇异矩阵方程求解:
#[[0.0714285714285714],
[0.142857142857143],
[0.214285714285714]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
from scipy import linalg as dense_linalg
import sympy as sp
import numpy as np
def mrdivide(input_str):
"""
对标Matlab的mrdivide函数，求解矩阵方程 xA = B 的解x。
参数:
B: 右侧矩阵，可以是SymPy矩阵、列表的列表或NumPy数组。
A: 系数矩阵，可以是SymPy矩阵、列表的列表或NumPy数组。
返回:
x: SymPy矩阵形式的最小二乘解或精确解。
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
A = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
B = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
if A is not None and B is not None:
# 将 SymPy 矩阵转换为 NumPy 二维数组，保持结构
# 转换为NumPy数组进行数值计算
if A.shape[1] == 1:
N_A = np.array(np.ravel(A), dtype=float)
else:
N_A = np.array(A, dtype=float)
if B.shape[1] == 1:
N_B = np.array(np.ravel(B), dtype=float)
else:
N_B = np.array(B, dtype=float)
# 转置 A 和 B
A_T = N_A.T
B_T = N_B.T
try:
# 尝试直接求解 Ax = B（适用于方阵且可逆的情况）
x_T = np.linalg.solve(A_T, B_T)
except np.linalg.LinAlgError:
# 若无法直接求解，使用最小二乘法（处理超定、欠定等情况）
x_T, _, _, _ = np.linalg.lstsq(A_T, B_T, rcond=None)
# 转置回原来的 x
result = x_T.T
else:
error = True
else:
error = True
return sp.Matrix(result) if not error else f"错误"
except Exception as e:
return f"错误: {e}"
# 示例代码
if __name__ == "__main__":
# 示例1: 精确解
x = mrdivide("[[1,1,3],[2,0,4],[-1,6,-1]],[2,19,8]")
print("示例1的解 (精确解):")
print(x)
# Matrix([[0.999999999999999],
# [2.00000000000000],
# [3.00000000000000]])
# 示例2: 最小二乘解
x_ls = mrdivide("[[1,0],[2,0],[1,0]],[1,2]")
print("\n示例2的解 (最小二乘解):")
print(x_ls)
# Matrix([[0.166666666666667],
# [0.333333333333333],
# [0.166666666666667]])
矩阵乘法
C=A*B是A和B的矩阵乘积. 如果A是一个m-by-p矩阵，B是一个p-by-n矩阵.
这个定义说C(i，j)是A的第i行与B的第j列的内积
A,B -- 操作对象,标量,向量,矩阵
C--产品,标量,向量,矩阵
示例1. 矩阵乘法函数测试
基本矩阵乘法:
result1 = mtimes([[1,2], [3,4]], [[5,6], [7,8]])
print(f"结果: {result1}\n")
#结果:
#[[19, 22],
[43, 50]]
符号矩阵乘法
result2 = mtimes([[x,2*x**2,3*x**3,4*x**4]], [[1/x],[2/x**2],[3/x**3],[4/x**4]])
print(f"结果: {result2}\n")
#结果: [[30]]
示例2. 线性变换应用
二维旋转
theta = 45 # 45度
cos_t = math.cos(math.radians(theta))
sin_t = math.sin(math.radians(theta))
rotation_matrix = [[cos_t,-sin_t],[sin_t,cos_t]]
point = [[2],[0]] # 点 (2,0)
result = mtimes(rotation_matrix, point)
print(f"旋转矩阵 R(45°): {rotation_matrix}")
print(f"原始点: {point}")
print(f"旋转后点: {result}")
#旋转矩阵 R(45°): [[0.707,-0.707],[0.707,0.707]]
#原始点: [[2],[0]]
#旋转后点: [[1.414],
[1.414]]
#说明: 矩阵乘法实现坐标旋转
三维变换组合
先缩放，再旋转，最后平移
scale_matrix = [[2,0,0],[0,3,0],[0,0,1]]
rotation_matrix_3d = [[1,0,0],[0,0.866,-0.5],[0,0.5,0.866]] # 绕x轴旋转30度
translation_matrix = [[1,0,1],[0,1,2],[0,0,1]]
point_3d = [[1],[1],[1]]
组合变换: 平移 * 旋转 * 缩放
step1 = mtimes(scale_matrix, point_3d)
step2 = mtimes(rotation_matrix_3d, step1)
result = mtimes(translation_matrix, step2)
print(f"缩放矩阵: {scale_matrix}")
print(f"旋转矩阵: {rotation_matrix_3d}")
print(f"平移矩阵: {translation_matrix}")
print(f"原始点: {point_3d}")
print(f"变换后点: {result}")
#缩放矩阵: [[2,0,0],[0,3,0],[0,0,1]]
#旋转矩阵: [[1,0,0],[0,0.866,-0.5],[0,0.5,0.866]]
#平移矩阵: [[1,0,1],[0,1,2],[0,0,1]]
#原始点: [[1],[1],[1]]
#变换后点: [[4.366],
[6.83],
[2.366]]
#说明: 多个变换矩阵的复合
投影变换
projection_matrix = [[1,0,0],[0,1,0],[0,0,0]] # 投影到xy平面
point_3d = [[3],[4],[5]]
result = mtimes(projection_matrix, point_3d)
print(f"投影矩阵: {projection_matrix}")
print(f"三维点: {point_3d}")
print(f"投影点: {result}")
#投影矩阵: [[1,0,0],[0,1,0],[0,0,0]]
#三维点: [[3],[4],[5]]
#投影点: [[3],
[4],
[0]]
#说明: 将3D点投影到2D平面
示例3. 图论应用
邻接矩阵幂 - 路径计数
4个节点的有向图
adjacency_matrix = [[0,1,1,0],[1,0,0,1],[0,1,0,1],[1,0,0,0]]
result = mtimes(adjacency_matrix, adjacency_matrix)
print(f"邻接矩阵 A: {adjacency_matrix}")
print(f"A²: {result}")
#邻接矩阵 A: [[0,1,1,0],[1,0,0,1],[0,1,0,1],[1,0,0,0]]
#A²: [[1, 1, 0, 2],
[1, 1, 1, 0],
[2, 0, 0, 1],
[0, 1, 1, 0]]
#说明: A²[i,j] 表示从节点i到节点j长度为2的路径数量
可达性矩阵
计算 A + A² + A³
A_squared = mtimes(adjacency_matrix, adjacency_matrix)
A_cubed = mtimes(adjacency_matrix, A_squared)
print(f"A³: {A_cubed}")
#A³: [[3, 1, 1, 1],
[1, 2, 1, 2],
[1, 2, 2, 0],
[1, 1, 0, 2]]
#说明: A³[i,j] 表示从节点i到节点j长度为3的路径数量
示例4. 经济学应用
投入产出分析
技术系数矩阵
technology_matrix = [[0.2,0.1,0.3],[0.3,0.2,0.1],[0.1,0.3,0.2]]
最终需求
final_demand = [[100],[200],[150]]
总产出 = (I - A)⁻¹ * 最终需求
identity_minus_A = [[0.8,-0.1,-0.3],[-0.3,0.8,-0.1],[-0.1,-0.3,0.8]]
result = mtimes(identity_minus_A, final_demand)
print(f"技术系数矩阵 A: {technology_matrix}")
print(f"最终需求: {final_demand}")
print(f"总产出: {result}")
#技术系数矩阵 A: [[0.2,0.1,0.3],[0.3,0.2,0.1],[0.1,0.3,0.2]]
#最终需求: [[100],[200],[150]]
#总产出: [[15],
[115],
[50]]
#说明: 列昂惕夫投入产出模型
马尔可夫链
状态转移矩阵
transition_matrix = [[0.7,0.3],[0.4,0.6]]
初始状态分布
initial_state = [[0.5],[0.5]]
2步后的状态分布
result = mtimes(transition_matrix, initial_state)
result_2step = mtimes(transition_matrix, result)
print(f"转移矩阵 P: {transition_matrix}")
print(f"初始状态: {initial_state}")
print(f"1步后状态: {result}")
print(f"2步后状态: {result_2step}")
#转移矩阵 P: [[0.7,0.3],[0.4,0.6]]
#初始状态: [[0.5],[0.5]]
#1步后状态: [[0.5],
[0.5]]
#2步后状态: [[0.5],
[0.5]]
#说明: 马尔可夫链状态预测
示例5. 物理学应用
量子力学 - 算符作用
泡利矩阵
pauli_x = [[0,1],[1,0]]
量子态
quantum_state = [[1/sqrt(2)],[1/sqrt(2)]]
result = mtimes(pauli_x, quantum_state)
print(f"泡利X矩阵: {pauli_x}")
print(f"量子态: {quantum_state}")
print(f"作用后状态: {result}")
#泡利X矩阵: [[0,1],[1,0]]
#量子态: [[1/sqrt(2)],[1/sqrt(2)]]
#作用后状态: [[0.707106781186548],
[0.707106781186548]]
#说明: 量子算符对量子态的作用
刚体转动
惯性张量
inertia_tensor = [[10,0,0],[0,15,0],[0,0,12]]
角速度
angular_velocity = [[2],[3],[1]]
角动量 = 惯性张量 * 角速度
result = mtimes(inertia_tensor, angular_velocity)
print(f"惯性张量 I: {inertia_tensor}")
print(f"角速度 ω: {angular_velocity}")
print(f"角动量 L = Iω: {result}")
#惯性张量 I: [[10,0,0],[0,15,0],[0,0,12]]
#角速度 ω: [[2],[3],[1]]
#角动量 L = Iω: [[20],
[45],
[12]]
#说明: 刚体动力学中的角动量计算
示例6. 计算机图形学应用
模型视图变换
模型矩阵
model_matrix = [[1,0,0,2],[0,1,0,3],[0,0,1,1],[0,0,0,1]]
视图矩阵
view_matrix = [[1,0,0,-1],[0,1,0,-2],[0,0,1,-1],[0,0,0,1]]
顶点坐标
vertex = [[1],[2],[3],[1]]
组合变换: 视图矩阵 * 模型矩阵 * 顶点
model_view = mtimes(view_matrix, model_matrix)
result = mtimes(model_view, vertex)
print(f"模型矩阵: {model_matrix}")
print(f"视图矩阵: {view_matrix}")
print(f"顶点: {vertex}")
print(f"变换后顶点: {result}")
#模型矩阵: [[1,0,0,2],[0,1,0,3],[0,0,1,1],[0,0,0,1]]
#视图矩阵: [[1,0,0,-1],[0,1,0,-2],[0,0,1,-1],[0,0,0,1]]
#顶点: [[1],[2],[3],[1]]
#变换后顶点: [[2],
[3],
[3],
[1]]
#说明: 3D图形中的模型视图变换
颜色混合
RGB颜色混合矩阵
color_matrix = [[0.3,0.6,0.1],[0.2,0.7,0.1],[0.1,0.1,0.8]]
原始颜色
original_color = [[255],[128],[64]] # RGB值
result = mtimes(color_matrix, original_color)
print(f"颜色混合矩阵: {color_matrix}")
print(f"原始颜色: {original_color}")
print(f"混合后颜色: {result}")
#颜色混合矩阵: [[0.3,0.6,0.1],[0.2,0.7,0.1],[0.1,0.1,0.8]]
#原始颜色: [[255],[128],[64]]
#混合后颜色: [[159.7],
[147],
[89.5]]
#说明: 图像处理中的颜色变换
示例7. 机器学习应用
神经网络前向传播
权重矩阵
weights = [[0.1,0.2,0.3],[0.4,0.5,0.6]]
输入向量
inputs = [[0.5],[0.3],[0.8]]
result = mtimes(weights, inputs)
print(f"权重矩阵 W: {weights}")
print(f"输入向量 x: {inputs}")
print(f"输出 z = Wx: {result}")
#权重矩阵 W: [[0.1,0.2,0.3],[0.4,0.5,0.6]]
#输入向量 x: [[0.5],[0.3],[0.8]]
#输出 z = Wx: [[0.35],
[0.83]]
#说明: 神经网络单层前向传播
主成分分析
数据协方差矩阵
covariance = [[4,2,1],[2,3,1],[1,1,2]]
特征向量矩阵
eigenvectors = [[0.7,0.5,0.3],[0.5,-0.7,0.4],[0.3,0.4,-0.9]]
数据投影
data_point = [[2],[1],[3]]
result = mtimes(eigenvectors, data_point)
print(f"特征向量矩阵: {eigenvectors}")
print(f"数据点: {data_point}")
print(f"主成分坐标: {result}")
#特征向量矩阵: [[0.7,0.5,0.3],[0.5,-0.7,0.4],[0.3,0.4,-0.9]]
#数据点: [[2],[1],[3]]
#主成分坐标: [[2.8],
[1.5],
[-1.7]]
#说明: PCA数据降维
示例8. 符号计算应用
多项式矩阵
A = [[a,b],[c,d]]
B = [[x],[y]]
result = mtimes(A, B)
print(f"矩阵 A: {A}")
print(f"向量 B: {B}")
print(f"乘积 AB: {result}\n")
#矩阵 A: [[a,b],[c,d]]
#向量 B: [[x],[y]]
#乘积 AB: [[a*x + b*y],
[c*x + d*y]]
参数化变换
rotation_param = [[cos(theta),-sin(theta)],[sin(theta),cos(theta)]]
point = [[r*cos(phi)],[r*sin(phi)]]
result = mtimes(rotation_param, point)
print(f"旋转矩阵: {rotation_param}")
print(f"极坐标点: {point}")
print(f"旋转后坐标: {result}\n")
#旋转矩阵: [[cos(theta),-sin(theta)],[sin(theta),cos(theta)]]
#极坐标点: [[r*cos(phi)],[r*sin(phi)]]
#旋转后坐标: [[-r*sin(phi)*sin(theta) + r*cos(phi)*cos(theta)],
[r*sin(phi)*cos(theta) + r*sin(theta)*cos(phi)]]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
import math
def matrix_multiplication(input_str):
"""
执行矩阵乘法（对标 MATLAB 的 mtimes 函数）
参数:
input_str (str): 输入表达式，格式为 (Matrix_A, Matrix_B)
返回:
sp.Matrix/str: 乘积矩阵或错误信息
示例:
>>> matrix_multiplication("([[1,2], [3,4]], [[5,6], [7,8]])")
Matrix([
[19, 22],
[43, 50]])
>>> matrix_multiplication("([a,b], [c;d])") # 符号计算
Matrix([[a*c + b*d]])
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
if isinstance(expr, tuple) and len(expr) == 2:
A = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
B = sp.Matrix(expr[1]) if isinstance(expr[1], list) else None
# 检查矩阵是否包含符号或复数
contains_symbols = any(
element.has(sp.Symbol) or element.has(sp.I)
for matrix in [A, B]
for element in matrix
)
# 符号矩阵使用 SymPy 乘法
if contains_symbols:
result = A * B # SymPy 自动检查维度
# 数值矩阵使用 NumPy 加速
else:
np_A = np.array(A.tolist(), dtype=np.float64)
np_B = np.array(B.tolist(), dtype=np.float64)
np_C = np_A @ np_B # NumPy 矩阵乘法
result = sp.Matrix(np_C.tolist())
else:
error = True
return result if not error else f"输入错误: {input_str}"
except sp.matrices.common.ShapeError as e:
return f"错误:{e}"
except Exception as e:
return f"错误:{e}"
# 测试用例
print(matrix_multiplication("([[x,2*x**2,3*x**3,4*x**4]], [[1/x],[2/x**2],[3/x**3],[4/x**4]])"))
# Matrix([[30]])
print(matrix_multiplication("([[1,2], [3,4]], [[5,6], [7,8]])"))
# Matrix([[19.0000000000000, 22.0000000000000],
# [43.0000000000000, 50.0000000000000]])
print(matrix_multiplication("([[2]], [[1,2,3]])"))
# Matrix([[2.00000000000000, 4.00000000000000, 6.00000000000000]])
多项式展开式的系数
multinomial(n,k1,k2,k3)返回多项式展开式的系数.
n, k1, k2, k3 -- 标量, 符号变量, 复数.
示例1. 多项式系数函数测试
基本多项式系数
result1 = multinomial(2,1,1)
print(f"结果: {result1}\n")
#结果: 12.0
二项式系数特例
result2 = multinomial(3,2)
print(f"结果: {result2}\n")
#结果: 10.0
符号计算
result3 = multinomial(a,b)
print(f"结果: {result3}\n")
#结果: factorial(a + b)
示例2. 组合数学应用
1. 排列组合问题
将5个不同的球放入3个盒子，第一个盒子放2个，第二个盒子放2个，第三个盒子放1个
result = multinomial(2,2,1)
print(f"分配方案数 (2,2,1): {result}")
print("说明: 5个球分配到3个盒子的不同方法数\n")
#分配方案数 (2,2,1): 30.0
#说明: 5个球分配到3个盒子的不同方法数
2. 字符串排列
单词 "MISSISSIPPI" 中字母的排列数
M:1, I:4, S:4, P:2
result = multinomial(1,4,4,2)
print(f"'MISSISSIPPI' 排列数: {result}")
#'MISSISSIPPI' 排列数: 34650.0
#说明: 计算有重复字母的单词的不同排列数
3. 多项展开式系数
(x+y+z)^4 展开式中 x²yz 的系数
result = multinomial(2,1,1)
print(f"(x+y+z)⁴ 中 x²yz 的系数: {result}")
#(x+y+z)⁴ 中 x²yz 的系数: 12.0
#说明: 多项式展开中特定项的系数
示例3. 概率论应用
1. 多项分布概率
掷骰子10次，得到1点3次，2点2次，3点2次，4点1次，5点1次，6点1次的概率系数
result = multinomial(3,2,2,1,1,1)
probability_coeff = result * ((1 / 6) ** 10)
print(f"骰子实验系数: {result}")
print(f"概率系数部分: {probability_coeff:.8f}")
#骰子实验系数: 151200.0
#概率系数部分: 0.00250057
#说明: 多项分布概率计算中的组合系数
2. 遗传学概率
孟德尔遗传实验：在8个后代中，3个显性纯合，2个杂合，3个隐性纯合
result = multinomial(3,2,3)
print(f"遗传分布系数: {result}")
#遗传分布系数: 560.0
#说明: 遗传学中基因型分布的组合系数
3. 抽样调查
从3个群体中抽取样本：群体A抽4人，群体B抽3人，群体C抽3人
result = multinomial(4,3,3)
print(f"抽样方案数: {result}")
print("说明: 分层抽样中的组合方案数\n")
#抽样方案数: 4200.0
#说明: 分层抽样中的组合方案数
示例4. 统计学应用
1. 列联表分析
3×3列联表的边际分布系数
row_totals = (5,4,6) # 各行总和
col_totals = (3,4,8) # 各列总和
result = multinomial(5,4,6)
print(f"行分布系数: {result}")
#行分布系数: 630630.0
#说明: 列联表分析中的多项系数
2. 聚类分析
将10个数据点分配到3个聚类：聚类1有4个点，聚类2有3个点，聚类3有3个点
result = multinomial(4,3,3)
print(f"聚类分配方案数: {result}")
#聚类分配方案数: 4200.0
#说明: 聚类分析中数据点分配的组合数
示例5. 物理学应用
1. 统计力学 - 粒子分布
10个不可区分粒子在3个能级上的分布：能级1有4个粒子，能级2有3个粒子，能级3有3个粒子
result = multinomial(4,3,3)
print(f"玻色子分布微观状态数: {result}")
#玻色子分布微观状态数: 4200.0
#说明: 玻色-爱因斯坦统计中的分布系数
2. 量子力学 - 态计数
多个量子态中的粒子分布
result = multinomial(2,3,1,4)
print(f"量子态分布系数: {result}")
#量子态分布系数: 12600.0
#说明: 多能级系统中粒子分布的状态数
示例6. 计算机科学应用
1. 数据分区
将数据分配到3个服务器：服务器1接收40%数据，服务器2接收35%数据，服务器3接收25%数据
在100个数据块中的分配方案数
result = multinomial(40,35,25)
print(f"数据分区方案数: {result}")
#数据分区方案数: 7.13641766177020E+44
#说明: 分布式系统中数据分布的组合数
2. 负载均衡
10个任务分配到4个处理器：处理器1处理3个，处理器2处理2个，处理器3处理3个，处理器4处理2个
result = multinomial(3,2,3,2)
print(f"任务分配方案数: {result}")
print("说明: 负载均衡中任务分配的组合数\n")
#任务分配方案数: 25200.0
#说明: 负载均衡中任务分配的组合数
3. 编码理论
编码中特定模式的排列数
result = multinomial(3,2,2,1,2)
print(f"编码模式数: {result}")
#编码模式数: 75600.0
#说明: 错误检测编码中的组合模式数
示例7. 经济学应用
1. 资源分配
将100单位资源分配到4个项目：项目1得30，项目2得25，项目3得25，项目4得20
result = multinomial(30,25,25,20)
print(f"资源分配方案数: {result}")
#资源分配方案数: 6.01073533124950E+56
#说明: 经济学中资源分配的组合方法数
2. 市场份额分析
有4个公司在市场中的订单分布
result = multinomial(5,3,4,3)
print(f"订单分布模式数: {result}")
#订单分布模式数: 12612600.0
#说明: 市场竞争中订单分布的组合分析
示例8. 生物学应用
1. 生态学 - 物种分布
在5个区域中观察到的物种个体数分布
result = multinomial(3,2,4,1,5)
print(f"物种分布模式数: {result}")
#物种分布模式数: 37837800.0
#说明: 生态学中种群分布的组合数
2. 基因组学
DNA序列中碱基的特定分布模式
result = multinomial(2,3,2,3)
print(f"碱基分布模式数: {result}")
#碱基分布模式数: 25200.0
#说明: 基因组序列分析的组合系数
示例9. 符号分析应用
1. 一般多项式系数
result1 = multinomial(n,k,m)
print(f"符号系数 (n,k,m): {result1}\n")
#符号系数 (n,k,m): factorial(k + m + n)
2. 参数化分布
result2 = multinomial(a,b,c)
print(f"参数化系数 (a,b,c): {result2}\n")
#参数化系数 (a,b,c): factorial(a + b + c)
示例10. 高等数学应用
1. 泰勒展开
多元函数的泰勒展开系数
result = multinomial(2,1,1)
print(f"三变量泰勒展开系数: {result}")
#三变量泰勒展开系数: 12.0
#说明: 多元函数泰勒展开中的组合系数
2. 微分方程
偏微分方程解的系数
result = multinomial(3,2,1)
print(f"偏微分方程解系数: {result}")
#偏微分方程解系数: 60.0
#说明: 某些偏微分方程级数解中的系数
示例11. 实际综合问题
1. 选举结果分析
有4个候选人在10个选区的得票分布
result = multinomial(3,2,3,2)
print(f"选举结果分布模式数: {result}")
#选举结果分布模式数: 25200.0
#说明: 政治学中选举结果的分析
2. 生产计划
工厂生产4种产品，每种产品的产量分配
result = multinomial(5,3,4,3)
print(f"生产计划方案数: {result}")
#生产计划方案数: 12612600.0
#说明: 运营管理中的生产调度
3. 交通流量
有4条道路的车辆分布
result = multinomial(4,3,5,3)
print(f"交通流分布模式数: {result}")
#交通流分布模式数: 12612600.0
#说明: 交通工程中的流量分析
示例12. 边界情况测试
1. 包含零的分布
result2 = multinomial(3,0,2)
print(f"包含零的分布 (3,0,2): {result2}\n")
#包含零的分布 (3,0,2): 10.0
2. 大数计算
result3 = multinomial(10,10,10)
print(f"大数计算 (10,10,10): {result3}\n")
#大数计算 (10,10,10): 5550996791340.00
3. 分数参数
result4 = multinomial(1/2,1/3,1/6)
print(f"分数参数 (1/2,1/3,1/6): {result4}\n")
#分数参数 (1/2,1/3,1/6): 1.36206225271963
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
def multinomial_coefficient(input_str):
"""
计算多项式系数（对标 MathStudio 的 Multinomial 函数）
参数:
input_str (str): 输入表达式，格式为 (k1, k2, ..., km)，其中 sum(ki) = n
返回:
float/sympy.Expr: 多项式系数 n!/(k1!k2!...km!)，若含符号则保留表达式
示例:
>>> multinomial_coefficient("(2,1,1)") # 计算 4!/(2!1!1!) = 12
12.0
>>> multinomial_coefficient("(1,1)") # 计算 2!/(1!1!) = 2
2.0
>>> multinomial_coefficient("(a,b)") # 符号计算 (a+b)!/(a!b!)
factorial(a + b)/(factorial(a)*factorial(b))
"""
try:
expr = sp.sympify(input_str)
error = False
result = None
def factorial_expansion(n, total_sum):
# 计算 (n + total_sum)!
return sp.factorial(n + total_sum)
if isinstance(expr, tuple):
n, *args = expr
total_sum = sum(args)
# 将 args 转换为元组，用于组合成 all_factors
args_tuple = tuple(args)
# 将 n 和 args 组合成一个元组 all_factors
all_factors = (n,) + args_tuple
# 计算总数的阶乘 (n + total_sum)!
factorial_total = factorial_expansion(n, total_sum)
# 计算需要除以的乘积，跳过符号变量
divisor_product = sp.prod(sp.factorial(k) for k in all_factors if not isinstance(k, sp.Symbol))
# 计算多项式系数 (n + total_sum)! / (除以的乘积)
coefficient = factorial_total / divisor_product
if n.is_symbol:
result = coefficient
else:
result = coefficient.evalf()
else:
error = True
return result if not error else f"输入格式错误: 应为类似 (2,1,1) 的元组形式"
except Exception as e:
return f"错误: {str(e)}"
print(multinomial_coefficient("x,1/3,3"))
# factorial(x + 10/3)/(6*factorial(1/3))
多元正态累积分布函数
p=mvncdf（X）返回具有零均值和单位协方差矩阵的多元正态分布的累积分布函数（cdf），在X的每一行进行评估。有关更多信息，请参阅多元正态分配。
p=mvncdf（X，mu，Sigma）返回具有均值mu和协方差Sigma的多元正态分布的cdf，在X的每一行进行评估。
当您只想指定Sigma时，为mu指定[]以使用其默认值零。
X —— 评价点,数值矩阵
mu —— 多元正态分布的均值向量, 零向量（默认）|数值向量|数值标量
Sigma —— 多元正态分布的协方差矩阵, 恒等矩阵（默认）|对称正定矩阵|对角条目的数值向量
示例1. 计算投资组合在给定阈值下的损失概率
假设两个资产：股票和债券，相关系数为0.3
阈值：股票下跌10%，债券下跌5%
result = mvncdf([-0.1, -0.05], [0.02, 0.01], [[0.04, 0.006], [0.006, 0.01]])
print(f"投资组合损失概率: {result}")
#投资组合损失概率: 0.110632
示例2. 计算产品多个尺寸同时满足规格的概率
假设长度和直径相关
规格：长度 10±0.2mm，直径 5±0.1mm
result = mvncdf([10.2, 5.1], [10, 5], [[0.01, 0.002], [0.002, 0.0025]])
print(f"产品合格率: {1 - result:.4f}")
#产品合格率: 0.0426
示例3. 计算系统中多个相关组件同时失效的概率
组件A寿命≤800小时，组件B寿命≤1000小时
result = mvncdf([800, 1000], [1000, 1200], [[10000, 3000], [3000, 25000]])
print(f"系统组件同时失效概率: {result:.6f}")
#系统组件同时失效概率: 0.004638
示例4. 额外的边界情况测试
一维标准正态:
print(mvncdf([0], [0], [[1]]))
#0.5
二维独立标准正态:
print(mvncdf([0, 0], [0, 0], [[1, 0], [0, 1]]))
#0.25
多个点同时计算:
print(mvncdf([[1, 2, 3], [1, 2, 3]], [0, 0], [[1, 0], [0, 1]]))
#[0.707861, 0.955017, 0.997302]
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.stats import norm, mvn
def multivariate_normal_cdf(input_str):
"""
计算多元正态累积分布函数（对标 MATLAB mvncdf）
参数:
input_str (str): 输入表达式，格式为 ([数据矩阵], [均值向量], [协方差矩阵])
返回:
float/list: CDF 值或错误信息
示例:
>>> multivariate_normal_cdf("([0,0])") # 2D标准正态
0.25
>>> multivariate_normal_cdf("([1.96], [0], [[1]])") # 1D 95%置信区间
0.975
"""
try:
expr = sp.sympify(input_str)
args = expr if isinstance(expr, tuple) else (expr,)
n_args = len(args)
if n_args < 1 or n_args > 3:
return "参数错误：需1-3个参数 (X, mu, Sigma)"
# 解析数据矩阵 X（转换为 NumPy 并检查维度）
X_sym = sp.Matrix(args[0]) if isinstance(args[0], list) else None
if X_sym is None:
return "X 应为二维矩阵（例如 [[x1,x2,...]]）"
X = np.array(X_sym.tolist()).astype(float) # 关键转换
if X.ndim != 2:
return "X 应为二维矩阵"
d, n_samples = X.shape[0], X.shape[1]
# 解析均值 mu（转换为 NumPy）
mu = np.zeros(d)
if n_args >= 2:
mu_sym = sp.Matrix(args[1]) if isinstance(args[1], list) else None
if mu_sym is None:
return f"mu 应为 ({d}, 1) 向量"
mu = np.array(mu_sym.tolist()).astype(float) # 关键转换
if mu.shape != (d, 1):
return f"mu 应为 ({d}, 1) 向量"
mu = mu.flatten()
# 解析协方差 Sigma（转换为 NumPy）
Sigma = np.eye(d)
if n_args >= 3:
Sigma_sym = sp.Matrix(args[2]) if isinstance(args[2], list) else None
if Sigma_sym is None:
return f"Sigma 应为 ({d}, {d}) 矩阵"
Sigma = np.array(Sigma_sym.tolist()).astype(float) # 关键转换
if Sigma.shape != (d, d):
return f"Sigma 应为 ({d}, {d}) 矩阵"
if not np.allclose(Sigma, Sigma.T):
return "Sigma 必须对称"
if np.any(np.linalg.eigvalsh(Sigma) <= 1e-8): # 考虑浮点误差
return "Sigma 必须正定"
# 计算 CDF
results = []
for i in range(n_samples):
x = X[:, i]
if d == 1:
cdf = norm.cdf(x[0], loc=mu[0], scale=np.sqrt(Sigma[0, 0]))
results.append(round(cdf, 6))
elif d == 2:
p, _ = mvn.mvnun([-np.inf, -np.inf], x, mu, Sigma)
results.append(round(p, 6))
else:
return f"暂不支持 {d} 维计算"
return results if n_samples > 1 else results[0]
except Exception as e:
return f"错误: {str(e)}"
# 测试用例
print(multivariate_normal_cdf("([0,0])"))
# 0.25
print(multivariate_normal_cdf("([1.96], [0], [[1]])"))
# 0.975002
print(multivariate_normal_cdf("([0,0])"))
# 0.25
print(multivariate_normal_cdf("([[-2,-2,-2,-2],[-1,-1,-1,-1]])"))
# [0.003609, 0.003609, 0.003609, 0.003609]
多元正态概率密度函数
y = mvnpdf(X) 返回一个 n×1 向量 y，其中包含具有零均值和单位协方差矩阵的 d 维多元正态分布的概率密度函数 (pdf) 值，在 n×d 矩阵 X 的每行处进行计算。有关详细信息，请参阅多元正态分布。
y = mvnpdf(X,mu) 返回 X 中的点处的 pdf 值，其中 mu 确定每个关联的多元正态分布的均值。
y = mvnpdf(X,mu,Sigma) 返回 X 中的点处的 pdf 值，其中 Sigma 确定每个关联的多元正态分布的协方差。
当您只想指定 Sigma 时，为 mu 指定 [] 以使用其默认值零。
X — 计算点,数值向量 | 数值矩阵
mu — 多元正态分布的均值,零向量 (默认) | 数值向量 | 数值矩阵
Sigma — 多元正态分布的协方差, 单位矩阵 (默认) | 对称正定矩阵 | 数值数组
示例1. 金融风险管理 - 资产收益联合分布
计算两种资产在特定收益水平下的联合概率密度
股票收益3%，债券收益2%的概率密度
result = mvnpdf([0.03, 0.02], [0.08, 0.04], [[0.04, 0.01], [0.01, 0.01]])
print(f"股票3%收益 & 债券2%收益的概率密度: {result}")
#股票3%收益 & 债券2%收益的概率密度: 8.87276945930973
示例2. 质量控制 - 产品尺寸分布
计算产品多个关键尺寸在规格中心点的概率密度
产品在标称尺寸处的概率密度
result = mvnpdf([10.0, 5.0], [10.0, 5.0], [[0.0025, 0.0005], [0.0005, 0.0009]])
print(f"产品在标称尺寸的概率密度: {result}")
#产品在标称尺寸的概率密度: 112.539539519638
示例3. 医学诊断 - 生物标志物分布
计算健康人群生物标志物的典型值概率密度
健康人群典型生物标志物值的概率密度
result = mvnpdf([120, 80, 90], [120, 80, 90], [[25, 8, 5], [8, 16, 4], [5, 4, 9]])
print(f"健康人群典型生物标志物概率密度: {result}")
#健康人群典型生物标志物概率密度: 0.00125884321651651
示例4. 气象学 - 多地区温度分布
计算多个城市在典型温度下的联合概率密度
北京15°C，上海20°C，广州25°C的概率密度
result = mvnpdf([15, 20, 25], [15, 20, 25], [[9, 4, 2], [4, 16, 6], [2, 6, 25]])
print(f"典型温度分布概率密度: {result}")
#典型温度分布概率密度: 0.00117742430464722
示例5. 教育评估 - 学生成绩分布
计算学生在多个科目上取得典型成绩的概率密度
学生数学75分，物理70分，化学72分的概率密度
result = mvnpdf([75, 70, 72], [75, 70, 72], [[64, 25, 20], [25, 49, 18], [20, 18, 36]])
print(f"学生典型成绩概率密度: {result}")
#学生典型成绩概率密度: 0.000243379629189959
示例6. 市场营销 - 消费者行为分布
计算典型消费者的购买行为概率密度
消费者在三个产品类别上的典型支出
result = mvnpdf(
[400, 250, 150], [400, 250, 150], [[10000, 2000, 1500], [2000, 3600, 1200], [1500, 1200, 2500]])
print(f"典型消费者支出概率密度: {result}")
#典型消费者支出概率密度: 2.49619240311505E-7
示例7. 工程可靠性 - 组件寿命分布
计算系统组件在典型寿命下的概率密度
两个相关组件的典型寿命
result = mvnpdf([1000, 1200], [1000, 1200], [[10000, 3000], [3000, 25000]])
print(f"组件典型寿命概率密度: {result}")
#组件典型寿命概率密度: 0.0000102520711217092
示例8. 生态学 - 物种分布
计算物种在特定环境条件下的分布概率密度
物种在特定温度、湿度、海拔下的概率密度
result = mvnpdf([20, 60, 500], [20, 60, 500], [[4, 1, -0.5], [1, 25, 2], [-0.5, 2, 10000]])
print(f"物种适宜环境概率密度: {result}")
#物种适宜环境概率密度: 0.0000638142867108645
示例9. 边界情况和特殊测试
一维标准正态在0点的密度:
print(mvnpdf([0], [0], [[1]]))
#0.398942280401433
二维独立标准正态在(0,0)点的密度:
print(mvnpdf([0, 0], [0, 0], [[1, 0], [0, 1]]))
#0.159154943091895
示例10. 异常检测应用
正常操作参数
normal_point = ([100, 50, 25], [100, 50, 25], [[4, 1, 0.5], [1, 9, 1], [0.5, 1, 1]])
normal_pdf = mvnpdf(normal_point)
异常操作参数
anomaly_point = ([110, 60, 40], [100, 50, 25], [[4, 1, 0.5], [1, 9, 1], [0.5, 1, 1]])
anomaly_pdf = mvnpdf(anomaly_point)
print(f"正常操作概率密度: {normal_pdf}")
print(f"异常操作概率密度: {anomaly_pdf}")
print(f"异常检测比率: {anomaly_pdf / normal_pdf if isinstance(normal_pdf, (int, float)) and normal_pdf != 0 else 'N/A'}")
#正常操作概率密度: 0.0116409041263550
#异常操作概率密度: 1.16782640781281E-52
#异常检测比率: N/A
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
from sympy.stats import MultivariateNormal, density
def multivariate_normal_pdf(input_str):
"""
计算多元正态分布的概率密度函数，支持多个数据点。
参数:
input_str (str): 输入的字符串表达式，格式为元组，包含以下元素：
- X: 数据点矩阵（每列为一个数据点）
- mu (可选): 均值向量（默认零向量）
- Sigma (可选): 协方差矩阵（默认单位矩阵）
返回:
List[sp.Expr] | sp.Expr | str: 概率密度值列表，若输入错误则返回错误信息。
"""
try:
expr = sp.sympify(input_str)
if isinstance(expr, tuple):
# 解析参数
n_args = len(expr)
X = sp.Matrix(expr[0]) if isinstance(expr[0], list) else None
else:
n_args = 1
X = sp.Matrix(expr) if isinstance(expr, list) else None
if n_args < 1 or n_args > 3:
return f"输入错误：需要1-3个参数，当前输入{n_args}个参数"
# 解析X, mu, Sigma
if X is None:
return "第一个参数必须为矩阵"
n_dim, n_samples = X.rows, X.cols
# 处理均值向量
mu = sp.Matrix(expr[1]) if n_args >= 2 else None
Mu = sp.Matrix.zeros(n_dim, 1) if mu is None else mu
if Mu.rows != n_dim or Mu.cols != 1:
return f"均值向量维度错误，应为({n_dim}, 1)，当前为({Mu.rows}, {Mu.cols})"
# 处理协方差矩阵
sigma = sp.Matrix(expr[2]) if n_args >= 3 else None
Sigma = sp.eye(n_dim) if sigma is None else sigma
if Sigma.shape != (n_dim, n_dim):
return f"协方差矩阵维度错误，应为({n_dim}, {n_dim})，当前为{Sigma.shape}"
if not Sigma.is_symmetric():
return "协方差矩阵必须对称"
# 创建分布
M = MultivariateNormal('M', Mu, Sigma)
results = []
# 计算每个样本的PDF
for i in range(n_samples):
x = X[:, i]
try:
pdf = density(M)(x)
results.append(pdf.evalf() if pdf.has(sp.Float) else pdf)
except Exception as e:
results.append(f"Error at sample {i + 1}: {str(e)}")
return results if n_samples > 1 else results[0]
except Exception as e:
return f"错误：{str(e)}"
# 示例测试
if __name__ == "__main__":
# 示例1：默认参数（单个样本）
print("示例1 - 默认参数:")
print(multivariate_normal_pdf("([0, 0])"))
# 1/(2*pi)
# 示例2：一维多样本
print("\n示例2 - 一维多样本:")
print(multivariate_normal_pdf("([[1.0, 2.0]], [0.0], [[1.0]])"))
# [0.241970724519143, 0.0539909665131881]
多元正态随机数
R = mvnrnd(mu,Sigma,n) 返回一个矩阵 R，该矩阵由从同一多元正态分布（具有均值向量 mu 和协方差矩阵 Sigma）中选择的 n 个随机向量组成。有关详细信息，请参阅多元正态分布。
R = mvnrnd(mu,Sigma) 返回一个 m×d 矩阵 R，该矩阵由从 m 个单独的 d 维多元正态分布（其均值和协方差分别由 mu 和 Sigma 指定）中采样的随机向量组成。R 的每行均为单个多元正态随机向量。
mu — 多元正态分布的均值, 数值向量 | 数值矩阵
Sigma — 多元正态分布的协方差, 对称半正定矩阵 | 数值数组
n — 多元随机数的数量, 正整数标量
import pandas as pd
import matplotlib.pyplot as plt
示例1. 金融模拟 - 资产收益随机生成
参数：股票期望收益8%，债券期望收益4%，正相关
samples = mvnrnd([0.08, 0.04], [[0.04, 0.01], [0.01, 0.01]], 1000)
if not isinstance(samples, str): # 检查是否出错
df = pd.DataFrame(samples, columns=['股票收益', '债券收益'])
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.scatter(df['股票收益'], df['债券收益'], alpha=0.5)
plt.xlabel('股票收益')
plt.ylabel('债券收益')
plt.title('资产收益联合分布')
plt.grid(True, alpha=0.3)
plt.subplot(1, 2, 2)
plt.hist(df['股票收益'], bins=30, alpha=0.7, label='股票')
plt.hist(df['债券收益'], bins=30, alpha=0.7, label='债券')
plt.xlabel('收益')
plt.ylabel('频数')
plt.title('资产收益分布')
plt.legend()
plt.tight_layout()
plt.show()
print(f"股票收益均值: {np.mean(df['股票收益']):.4f}")
print(f"债券收益均值: {np.mean(df['债券收益']):.4f}")
print(f"收益相关系数: {np.corrcoef(df['股票收益'], df['债券收益'])[0, 1]:.4f}")
示例2. 质量控制 - 产品尺寸变异模拟
模拟生产过程中产品关键尺寸的随机变异
模拟长度和直径的制造变异（相关）
samples = mvnrnd([10.0, 5.0], [[0.0025, 0.0005], [0.0005, 0.0009]], 500)
if not isinstance(samples, str):
df = pd.DataFrame(samples, columns=['长度(mm)', '直径(mm)'])
plt.figure(figsize=(10, 6))
plt.scatter(df['长度(mm)'], df['直径(mm)'], alpha=0.6)
plt.axhline(5.0, color='red', linestyle='--', alpha=0.8, label='标称直径')
plt.axvline(10.0, color='blue', linestyle='--', alpha=0.8, label='标称长度')
plt.xlabel('长度 (mm)')
plt.ylabel('直径 (mm)')
plt.title('产品尺寸制造变异')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
# 计算合格率
length_ok = ((df['长度(mm)'] >= 9.9) & (df['长度(mm)'] <= 10.1)).mean()
diameter_ok = ((df['直径(mm)'] >= 4.9) & (df['直径(mm)'] <= 5.1)).mean()
both_ok = (length_ok * diameter_ok)
print(f"长度合格率: {length_ok:.2%}")
print(f"直径合格率: {diameter_ok:.2%}")
print(f"综合合格率: {both_ok:.2%}")
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.stats import multivariate_normal
def multivariate_normal_rnd(input_str):
"""
生成多元正态分布随机样本 (对标 MATLAB mvnrnd)
参数:
input_str (str): 输入表达式，格式为 (mu, Sigma, n) 或 (mu, Sigma)
返回:
list/np.ndarray: 随机样本列表或错误信息
示例:
>>> mvnrnd("([0], [[1]], 5)") # 1D，生成5个样本
>>> mvnrnd("([0,0], [[1,0],[0,1]])") # 2D，生成1个样本
"""
try:
expr = sp.sympify(input_str)
args = expr if isinstance(expr, tuple) else (expr,)
n_args = len(args)
if n_args not in [2, 3]:
return "参数错误：需2或3个参数 (mu, Sigma, n)"
# 解析均值 mu
mu_sym = sp.Matrix(args[0]) if isinstance(args[0], list) else None
if mu_sym is None or mu_sym.shape[1] != 1:
return "mu 应为列向量（例如 [0; 0]）"
mu = np.array(mu_sym.tolist()).astype(float).flatten()
d = mu.shape[0]
# 解析协方差 Sigma
Sigma_sym = sp.Matrix(args[1]) if isinstance(args[1], list) else None
if Sigma_sym is None or Sigma_sym.shape != (d, d):
return f"Sigma 应为 {d}x{d} 矩阵"
Sigma = np.array(Sigma_sym.tolist()).astype(float)
if not np.allclose(Sigma, Sigma.T):
return "Sigma 必须对称"
if np.any(np.linalg.eigvalsh(Sigma) <= 1e-8):
return "Sigma 必须正定"
# 解析样本数 n
n = 1
if n_args == 3:
try:
n = int(args[2])
if n <= 0:
return "样本数 n 必须为正整数"
except:
return "n 应为整数"
# 生成随机样本
if d == 1: # 1D 情况
samples = np.random.normal(mu[0], np.sqrt(Sigma[0, 0]), size=n)
elif d == 2: # 2D 情况
rv = multivariate_normal(mean=mu, cov=Sigma)
samples = rv.rvs(size=n)
else:
return "暂不支持超过2维"
return samples.tolist() if n > 1 else samples.tolist()[0]
except Exception as e:
return f"错误: {str(e)}"
# 1D 生成5个样本
print(multivariate_normal_rnd("([0], [[1]], 5)"))
# [1.39499095753052, 0.10507847625026556, -1.2253721122950096, -0.6349617714285217, 0.6486533888429011]
# 2D 生成1个样本
print(multivariate_normal_rnd("([0,0], [[1,0],[0,1]])"))
# -0.5636153325082512