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 — 解,向量,满矩阵,稀疏矩阵
# 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 — 多项式系数,矩阵
# 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 —— 样本数据和审查信息,向量,二维矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
from scipy.stats import binom, norm
from functools import partial
from scipy.optimize import minimize
def maximum_likelihood_estimate(input_str):
"""
执行最大似然估计，针对二项分布数据估计成功概率p。
支持两种输入格式：
1. 向量输入：每个元素表示成功次数k，默认试验次数n=1（伯努利试验）
2. 两列矩阵：第一列k为成功次数，第二列n为试验次数
参数:
input_str: 字符串形式的输入，可解析为向量或两列矩阵
返回:
元组 (phat, pci) 或错误信息:
- phat: 估计的成功概率
- pci: 95%置信区间
示例:
>>> maximum_likelihood_estimate("[1,0,1,1]")
(0.75, [0.5283..., 0.9716...])
"""
try:
# SymPy 表达式解析
expr = sp.sympify(input_str)
if isinstance(expr, list):
# 转换为numpy数组并预处理
data = np.array(expr, dtype=float)
else:
return f"输入解析错误: 无法将 {input_str} 转换为矩阵"
'''
elif (matrix := is_func(expr)) is not None:
data = np.array(matrix, dtype=float)
'''
# 维度处理：确保二维结构
if data.ndim == 1:
data = data.reshape(-1, 1)
# 列数校验与处理
if data.shape[1] == 1: # 向量输入
# 添加试验次数列（默认n=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"
# 定义负对数似然函数
def negative_log_likelihood(params, k, n):
p = params[0]
if not 0 < p < 1: # 边界检查
return np.inf
return -np.sum(binom.logpmf(k, n, p))
# 优化求解
result = minimize(
partial(negative_log_likelihood, k=k, n=n_values),
x0=[0.5],
bounds=[(1e-6, 1 - 1e-6)], # 数值稳定性处理
method='L-BFGS-B'
)
if not result.success:
return f"优化失败: {result.message}"
p_hat = result.x[0]
# 计算标准误差和置信区间
total_trials = np.sum(n_values)
se = np.sqrt(p_hat * (1 - p_hat) / total_trials)
z = norm.ppf(0.975) # 95% CI
ci = (p_hat - z * se, p_hat + z * se)
return (round(p_hat, 4), [round(ci[0], 4), round(ci[1], 4)])
except Exception as e:
return f"处理过程中发生错误: {str(e)}"
# 测试用例
if __name__ == "__main__":
# 伯努利试验（向量输入）
print(maximum_likelihood_estimate("[1,0,1,1,1]"))
# (0.8, [0.4494, 1.1506])
# 常规两列矩阵
print(maximum_likelihood_estimate("[[5,10],[6,10],[7,10]]"))
# (0.6, [0.4247, 0.7753])
# 混合试验次数
print(maximum_likelihood_estimate("[[3,5],[5,10],[2,8]]"))
# (0.4348, [0.2322, 0.6374])
多项式概率密度函数
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的每一行,或者在需要时复制它们.
# 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个多项式区间中每个区间的计数.
# 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 — 除数,标量,向量,矩阵
# 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"
# 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"
# 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是方阵,计算使用特征值分解.
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
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 — 解,向量,满矩阵,稀疏矩阵
# 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--产品,标量,向量,矩阵
# Copyright 2025 小塔立软件有限公司及其旗下网站:www.qikjik.com
# Licensed under the MIT License.
import sympy as sp
import numpy as np
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 -- 标量, 符号变量, 复数.
# 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 —— 多元正态分布的协方差矩阵, 恒等矩阵（默认）|对称正定矩阵|对角条目的数值向量
# 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 — 多元正态分布的协方差, 单位矩阵 (默认) | 对称正定矩阵 | 数值数组
# 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 — 多元随机数的数量, 正整数标量
# 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