尧图网站建设 
简介本资源为哈尔滨工业大学2023年秋季《计算建模》课程配套实验包面向计算机、数学、自动化等相关专业本科生及建模初学者聚焦数学建模→算法设计→编程实现的全链路实践训练。压缩包共17个文件1.3MB含10个Python源码文件覆盖Viterbi解码、HMM建模、EM算法、FFT/DCT图像处理、中值滤波去噪等典型实验、2张PNG/TIFF格式示例图像、1份Excel实验数据、1个CSV结果文件、1份Markdown说明书README.md及1个TIFF原始图像类型分布体现“理论-代码-数据-结果-文档”闭环结构。已有237人学习下载资源开放可修改提供完整实验任务框架、可运行代码与清晰说明便于读者理解隐马尔可夫模型、随机模拟、参数拟合、图像频域处理等核心知识点并通过动手调试深化对算法原理与工程实现间关系的认知。1. 这不是一份“交完就扔”的课程实验包而是哈工大计算建模课里真正能跑通HMM参数学习与序列解码的实操入口如果你在搜索“哈工大2023秋计算建模实验”时点开这个压缩包却发现解压后只有几个.py文件和一份PDF说明书——别急着关掉。它不是模板作业而是一套完整可复现的隐马尔可夫模型HMM教学实现从观测序列生成、初始参数设定到用EM算法迭代估计转移/发射概率再用Viterbi算法回溯最优隐状态路径。整个流程不依赖任何黑盒库如hmmlearn所有矩阵运算、对数似然计算、前向-后向递推都手写Python实现。适合两类人一是刚学完《数值分析》想验证HMM中矩阵迭代收敛性的哈工大学生二是需要快速搭建可调试HMM基线、理解EM收敛行为或Viterbi剪枝逻辑的算法工程师。它不追求工程封装但每行代码都对应教材公式——比如EM_algorithm.py里gamma和xi的更新逻辑直接映射到《统计学习方法》第10章的推导步骤。2. 为什么用纯NumPy重写HMM核心算法从哈工大课程设计目标看选型逻辑哈工大计算建模课程强调“可解释性优先于封装性”这决定了本实验不采用scikit-learn或PyTorch等高层框架。当学生需要调试EM算法中E步的后向概率溢出、M步的归一化失效或Viterbi路径回溯时索引越界黑盒API只会返回ValueError而手写实现能让你在print()里看到每一帧的alpha值、每轮的log_likelihood变化。这种设计直指HMM教学中的三个关键断点数值稳定性陷阱原始概率连乘易下溢必须转为log-space运算边界条件混淆Viterbi初始化时delta[0]应取log(π_i) log(b_i(o_0))而非直接乘EM收敛判据误设用参数差值而非对数似然增量判断收敛会导致过早终止。本实验包通过HMM.py定义基础类、Viterbi.py专注解码、EM_algorithm.py分离E/M步形成清晰职责边界。所有函数签名强制接收log_spaceTrue参数默认启用对数运算——这是哈工大数值分析课反复强调的“避免浮点灾难”实践。2.1 HMM类的结构设计为什么forward_log比forward更适合作为教学基线HMM.py中forward_log函数是整个流程的数值锚点。它不返回原始α_t(i)而是返回log_alpha[t][i]其递推式为log_alpha[t][i] logsumexp( [log_alpha[t-1][j] np.log(A[j][i]) for j in range(N)] ) np.log(B[i][O[t]])提示logsumexp是关键——它先减去最大值再指数求和再加回避免exp(-1000)导致的0.0。哈工大数值分析课中“防止下溢的补偿技巧”在此直接落地。对比传统forward实现需处理0.0除零forward_log天然规避了三类错误初始时刻t0时log(0)报错因π_i或B_i(o_0)为0t0时sum(α_{t-1}·A)结果为0导致后续log(0)多次迭代后α值趋近机器精度下限~1e-308引发NaN传播。实验说明书中明确要求所有概率运算必须经过np.log和logsumexp封装。这不是代码洁癖而是哈工大保研面试中常被追问的“如何保证HMM训练数值鲁棒性”的标准答案。2.2 Viterbi算法的路径回溯陷阱为什么psi数组必须用整数索引而非状态名Viterbi.py中viterbi_decode函数的psi数组定义为psi[t][i] argmax_j (delta[t-1][j] log(A[j][i]))其数据类型为int而非str。这个细节常被忽略却直接影响哈工大实验验收——当隐状态集为[Sunny, Rainy]时若psi存字符串回溯时无法用psi[t][i]作为delta[t-1]的索引因delta是数值数组。正确做法是将状态映射为{0:Sunny, 1:Rainy}psi[t][i]存储整数j即上一时刻最优前驱状态编号回溯时用path[t] psi[t1][path[t1]]最后用映射表转换为可读名。# 正确回溯逻辑摘自Viterbi.py path [0] * T path[-1] np.argmax(delta[-1]) for t in range(T-2, -1, -1): path[t] psi[t1][path[t1]] # psi[t1]是t1时刻记录的t时刻最优前驱 return [state_names[i] for i in path]注意psi维度为(T, N)但psi[0]无意义首时刻无前驱实际只用psi[1:T]。实验说明书中要求打印psi中间值验证正是为排查此索引偏移错误。3. 用真实观测序列跑通EM-Viterbi全流程从数据生成到参数收敛验证本实验包自带generate_data.py脚本可按指定π、A、B生成带标签的观测序列。但真正体现哈工大计算建模深度的是参数学习闭环验证用生成数据训练模型再对比学习出的A/B与真值的Frobenius范数误差。以下是在本地复现的最小命令链3.1 生成带标签的训练数据含隐状态真值python generate_data.py \ --n_states 3 \ --n_obs 4 \ --seq_len 100 \ --n_samples 50 \ --output_dir ./data/该命令生成50条长度为100的观测序列数字0~3同时保存对应隐状态序列0~2到./data/true_states.npy。关键参数说明--n_states隐状态数影响HMM类初始化时的N--n_obs观测符号数决定发射矩阵B的列数--seq_len单条序列长度过短会导致EM收敛震荡哈工大实验要求≥50--n_samples训练样本数少于20时EM易陷入局部最优。生成的数据结构为numpy.ndarray形状(50, 100)直接喂给EM_algorithm.py的train函数。3.2 执行EM算法并监控收敛过程from EM_algorithm import EMTrainer from HMM import HMM # 加载数据 X np.load(./data/observed_sequences.npy) # shape: (50, 100) true_A np.load(./data/true_transition.npy) # 真值用于对比 # 初始化HMM随机π,A,B hmm HMM(n_states3, n_obs4, log_spaceTrue) trainer EMTrainer(hmmhmm, max_iter100, tol1e-4) # 训练并获取每轮log_likelihood log_likelihoods, A_history, B_history trainer.train(X) # 验证收敛检查最后10轮log_likelihood波动是否tol converged np.std(log_likelihoods[-10:]) 1e-5 print(fEM converged: {converged}, final LL: {log_likelihoods[-1]:.4f})提示tol1e-4是哈工大实验默认阈值但若log_likelihoods出现平台期后突然下降说明E步中xi计算有误常见于未用logsumexp处理分母。3.3 用Viterbi解码并评估状态预测准确率训练完成后对同一批数据做解码# 加载真值隐状态 true_states np.load(./data/true_states.npy) # shape: (50, 100) # 解码预测 pred_states [] for seq in X: pred hmm.viterbi_decode(seq) # 返回list[int] pred_states.append(pred) pred_states np.array(pred_states) # (50, 100) # 计算整体准确率哈工大验收硬指标 accuracy np.mean(pred_states true_states) print(fViterbi accuracy: {accuracy:.4f})若准确率低于0.7需检查Viterbi.py中delta[0]是否用了log(π_i) log(B[i][o_0])psi数组是否在t0时未赋值应跳过logsumexp是否对delta[t-1] log(A[:,i])整体操作而非逐元素。4. 哈工大保研面试高频考点EM算法中E步的xi矩阵为何必须用log-space重写在EM_algorithm.py的e_step函数中xi[t][i][j]表示在t时刻处于状态i、t1时刻转移到状态j的概率。其原始公式为xi[t][i][j] alpha[t][i] * A[i][j] * B[j][o_{t1}] * beta[t1][j] / P(O|λ)但直接实现会因alpha/beta下溢导致xi全零。哈工大计算建模课要求将其转为log-spacelog_xi ( log_alpha[t][i] np.log(A[i][j]) np.log(B[j][O[t1]]) log_beta[t1][j] - log_likelihood ) xi[t][i][j] np.exp(log_xi) # 最终才转回概率这里log_likelihood是forward_log返回的总对数似然。面试官常追问“如果log_xi小于-700np.exp(log_xi)会是0.0此时xi矩阵稀疏M步更新会失效——如何避免”答案是在M步中不直接用xi而用其log-sum形式更新A# M步更新转移矩阵A[i][j] log_numerator logsumexp([ log_xi[t][i][j] for t in range(T-1) ]) log_denominator logsumexp([ log_gamma[t][i] for t in range(T-1) ]) A[i][j] np.exp(log_numerator - log_denominator)其中log_gamma[t][i] logsumexp([log_xi[t][i][j] for j in range(N)])。这种写法将数值问题完全隔离在log-space内直到最后一步才指数化——这正是哈工大数值分析课强调的“延迟指数化”原则。4.1 三个必调参数max_iter、tol、log_space的协同影响参数典型值调整逻辑哈工大实验约束max_iter50~100过小导致未收敛过大增加耗时但不提升精度≥80说明书明确要求tol1e-4~1e-6过大易早停过小在噪声数据中引发震荡1e-4保研面试常考此值合理性log_spaceTrue关闭则必然失败测试集已预设下溢场景必须为True否则forward_log报错当tol1e-6且max_iter200时若log_likelihoods在第150轮后波动1e-7说明模型已稳定但若第180轮突降大概率是logsumexp实现有缺陷如未减去最大值。5. 修改源码的实操技巧如何安全添加高斯发射概率以适配连续观测原实验包仅支持离散观测B[i][k]为状态i发射符号k的概率。若需处理哈工大数值分析课中的温度序列连续值需将HMM.py中的发射矩阵B替换为高斯分布参数每个状态i对应均值mu[i]和方差sigma2[i]。修改要点如下5.1 在HMM类中扩展_gaussian_emission_log方法def _gaussian_emission_log(self, state_i, obs_val): 计算状态i发射连续观测obs_val的log概率 return ( -0.5 * np.log(2 * np.pi * self.sigma2[state_i]) - 0.5 * ((obs_val - self.mu[state_i]) ** 2) / self.sigma2[state_i] )注意self.mu和self.sigma2需在__init__中初始化为np.random.randn(n_states)和np.random.rand(n_states)。5.2 改写forward_log中的观测概率项原离散版log_alpha[t][i] np.log(self.B[i][O[t]])改为连续版log_alpha[t][i] self._gaussian_emission_log(i, O[t])提示O[t]此时为float而非int需确保输入数据为np.float64。哈工大实验说明书中要求“连续观测需先标准化”即O (O - O.mean()) / O.std()否则mu/sigma2更新会发散。5.3 EM算法中M步的高斯参数更新公式在EM_algorithm.py的m_step中新增# 更新高斯均值mu[i] numerator sum( gamma[t][i] * O[t] for t in range(T) ) denominator sum(gamma[t][i] for t in range(T)) self.hmm.mu[i] numerator / denominator # 更新方差sigma2[i] numerator sum( gamma[t][i] * (O[t] - self.hmm.mu[i]) ** 2 for t in range(T) ) self.hmm.sigma2[i] numerator / denominator此处gamma[t][i]仍来自log-space的log_gamma但sum()操作在概率空间进行需先np.exp(log_gamma[t][i])。哈工大保研面试曾以此题考察“如何在log-space中安全执行加权平均”。本文还有配套的精品资源点击获取 哈工大HMM教学实现：纯NumPy手写EM与Viterbi算法
简介本资源为哈尔滨工业大学2023年秋季《计算建模》课程配套实验包面向计算机、数学、自动化等相关专业本科生及建模初学者聚焦数学建模→算法设计→编程实现的全链路实践训练。压缩包共17个文件1.3MB含10个Python源码文件覆盖Viterbi解码、HMM建模、EM算法、FFT/DCT图像处理、中值滤波去噪等典型实验、2张PNG/TIFF格式示例图像、1份Excel实验数据、1个CSV结果文件、1份Markdown说明书README.md及1个TIFF原始图像类型分布体现“理论-代码-数据-结果-文档”闭环结构。已有237人学习下载资源开放可修改提供完整实验任务框架、可运行代码与清晰说明便于读者理解隐马尔可夫模型、随机模拟、参数拟合、图像频域处理等核心知识点并通过动手调试深化对算法原理与工程实现间关系的认知。1. 这不是一份“交完就扔”的课程实验包而是哈工大计算建模课里真正能跑通HMM参数学习与序列解码的实操入口如果你在搜索“哈工大2023秋计算建模实验”时点开这个压缩包却发现解压后只有几个.py文件和一份PDF说明书——别急着关掉。它不是模板作业而是一套完整可复现的隐马尔可夫模型HMM教学实现从观测序列生成、初始参数设定到用EM算法迭代估计转移/发射概率再用Viterbi算法回溯最优隐状态路径。整个流程不依赖任何黑盒库如hmmlearn所有矩阵运算、对数似然计算、前向-后向递推都手写Python实现。适合两类人一是刚学完《数值分析》想验证HMM中矩阵迭代收敛性的哈工大学生二是需要快速搭建可调试HMM基线、理解EM收敛行为或Viterbi剪枝逻辑的算法工程师。它不追求工程封装但每行代码都对应教材公式——比如EM_algorithm.py里gamma和xi的更新逻辑直接映射到《统计学习方法》第10章的推导步骤。2. 为什么用纯NumPy重写HMM核心算法从哈工大课程设计目标看选型逻辑哈工大计算建模课程强调“可解释性优先于封装性”这决定了本实验不采用scikit-learn或PyTorch等高层框架。当学生需要调试EM算法中E步的后向概率溢出、M步的归一化失效或Viterbi路径回溯时索引越界黑盒API只会返回ValueError而手写实现能让你在print()里看到每一帧的alpha值、每轮的log_likelihood变化。这种设计直指HMM教学中的三个关键断点数值稳定性陷阱原始概率连乘易下溢必须转为log-space运算边界条件混淆Viterbi初始化时delta[0]应取log(π_i) log(b_i(o_0))而非直接乘EM收敛判据误设用参数差值而非对数似然增量判断收敛会导致过早终止。本实验包通过HMM.py定义基础类、Viterbi.py专注解码、EM_algorithm.py分离E/M步形成清晰职责边界。所有函数签名强制接收log_spaceTrue参数默认启用对数运算——这是哈工大数值分析课反复强调的“避免浮点灾难”实践。2.1 HMM类的结构设计为什么forward_log比forward更适合作为教学基线HMM.py中forward_log函数是整个流程的数值锚点。它不返回原始α_t(i)而是返回log_alpha[t][i]其递推式为log_alpha[t][i] logsumexp( [log_alpha[t-1][j] np.log(A[j][i]) for j in range(N)] ) np.log(B[i][O[t]])提示logsumexp是关键——它先减去最大值再指数求和再加回避免exp(-1000)导致的0.0。哈工大数值分析课中“防止下溢的补偿技巧”在此直接落地。对比传统forward实现需处理0.0除零forward_log天然规避了三类错误初始时刻t0时log(0)报错因π_i或B_i(o_0)为0t0时sum(α_{t-1}·A)结果为0导致后续log(0)多次迭代后α值趋近机器精度下限~1e-308引发NaN传播。实验说明书中明确要求所有概率运算必须经过np.log和logsumexp封装。这不是代码洁癖而是哈工大保研面试中常被追问的“如何保证HMM训练数值鲁棒性”的标准答案。2.2 Viterbi算法的路径回溯陷阱为什么psi数组必须用整数索引而非状态名Viterbi.py中viterbi_decode函数的psi数组定义为psi[t][i] argmax_j (delta[t-1][j] log(A[j][i]))其数据类型为int而非str。这个细节常被忽略却直接影响哈工大实验验收——当隐状态集为[Sunny, Rainy]时若psi存字符串回溯时无法用psi[t][i]作为delta[t-1]的索引因delta是数值数组。正确做法是将状态映射为{0:Sunny, 1:Rainy}psi[t][i]存储整数j即上一时刻最优前驱状态编号回溯时用path[t] psi[t1][path[t1]]最后用映射表转换为可读名。# 正确回溯逻辑摘自Viterbi.py path [0] * T path[-1] np.argmax(delta[-1]) for t in range(T-2, -1, -1): path[t] psi[t1][path[t1]] # psi[t1]是t1时刻记录的t时刻最优前驱 return [state_names[i] for i in path]注意psi维度为(T, N)但psi[0]无意义首时刻无前驱实际只用psi[1:T]。实验说明书中要求打印psi中间值验证正是为排查此索引偏移错误。3. 用真实观测序列跑通EM-Viterbi全流程从数据生成到参数收敛验证本实验包自带generate_data.py脚本可按指定π、A、B生成带标签的观测序列。但真正体现哈工大计算建模深度的是参数学习闭环验证用生成数据训练模型再对比学习出的A/B与真值的Frobenius范数误差。以下是在本地复现的最小命令链3.1 生成带标签的训练数据含隐状态真值python generate_data.py \ --n_states 3 \ --n_obs 4 \ --seq_len 100 \ --n_samples 50 \ --output_dir ./data/该命令生成50条长度为100的观测序列数字0~3同时保存对应隐状态序列0~2到./data/true_states.npy。关键参数说明--n_states隐状态数影响HMM类初始化时的N--n_obs观测符号数决定发射矩阵B的列数--seq_len单条序列长度过短会导致EM收敛震荡哈工大实验要求≥50--n_samples训练样本数少于20时EM易陷入局部最优。生成的数据结构为numpy.ndarray形状(50, 100)直接喂给EM_algorithm.py的train函数。3.2 执行EM算法并监控收敛过程from EM_algorithm import EMTrainer from HMM import HMM # 加载数据 X np.load(./data/observed_sequences.npy) # shape: (50, 100) true_A np.load(./data/true_transition.npy) # 真值用于对比 # 初始化HMM随机π,A,B hmm HMM(n_states3, n_obs4, log_spaceTrue) trainer EMTrainer(hmmhmm, max_iter100, tol1e-4) # 训练并获取每轮log_likelihood log_likelihoods, A_history, B_history trainer.train(X) # 验证收敛检查最后10轮log_likelihood波动是否tol converged np.std(log_likelihoods[-10:]) 1e-5 print(fEM converged: {converged}, final LL: {log_likelihoods[-1]:.4f})提示tol1e-4是哈工大实验默认阈值但若log_likelihoods出现平台期后突然下降说明E步中xi计算有误常见于未用logsumexp处理分母。3.3 用Viterbi解码并评估状态预测准确率训练完成后对同一批数据做解码# 加载真值隐状态 true_states np.load(./data/true_states.npy) # shape: (50, 100) # 解码预测 pred_states [] for seq in X: pred hmm.viterbi_decode(seq) # 返回list[int] pred_states.append(pred) pred_states np.array(pred_states) # (50, 100) # 计算整体准确率哈工大验收硬指标 accuracy np.mean(pred_states true_states) print(fViterbi accuracy: {accuracy:.4f})若准确率低于0.7需检查Viterbi.py中delta[0]是否用了log(π_i) log(B[i][o_0])psi数组是否在t0时未赋值应跳过logsumexp是否对delta[t-1] log(A[:,i])整体操作而非逐元素。4. 哈工大保研面试高频考点EM算法中E步的xi矩阵为何必须用log-space重写在EM_algorithm.py的e_step函数中xi[t][i][j]表示在t时刻处于状态i、t1时刻转移到状态j的概率。其原始公式为xi[t][i][j] alpha[t][i] * A[i][j] * B[j][o_{t1}] * beta[t1][j] / P(O|λ)但直接实现会因alpha/beta下溢导致xi全零。哈工大计算建模课要求将其转为log-spacelog_xi ( log_alpha[t][i] np.log(A[i][j]) np.log(B[j][O[t1]]) log_beta[t1][j] - log_likelihood ) xi[t][i][j] np.exp(log_xi) # 最终才转回概率这里log_likelihood是forward_log返回的总对数似然。面试官常追问“如果log_xi小于-700np.exp(log_xi)会是0.0此时xi矩阵稀疏M步更新会失效——如何避免”答案是在M步中不直接用xi而用其log-sum形式更新A# M步更新转移矩阵A[i][j] log_numerator logsumexp([ log_xi[t][i][j] for t in range(T-1) ]) log_denominator logsumexp([ log_gamma[t][i] for t in range(T-1) ]) A[i][j] np.exp(log_numerator - log_denominator)其中log_gamma[t][i] logsumexp([log_xi[t][i][j] for j in range(N)])。这种写法将数值问题完全隔离在log-space内直到最后一步才指数化——这正是哈工大数值分析课强调的“延迟指数化”原则。4.1 三个必调参数max_iter、tol、log_space的协同影响参数典型值调整逻辑哈工大实验约束max_iter50~100过小导致未收敛过大增加耗时但不提升精度≥80说明书明确要求tol1e-4~1e-6过大易早停过小在噪声数据中引发震荡1e-4保研面试常考此值合理性log_spaceTrue关闭则必然失败测试集已预设下溢场景必须为True否则forward_log报错当tol1e-6且max_iter200时若log_likelihoods在第150轮后波动1e-7说明模型已稳定但若第180轮突降大概率是logsumexp实现有缺陷如未减去最大值。5. 修改源码的实操技巧如何安全添加高斯发射概率以适配连续观测原实验包仅支持离散观测B[i][k]为状态i发射符号k的概率。若需处理哈工大数值分析课中的温度序列连续值需将HMM.py中的发射矩阵B替换为高斯分布参数每个状态i对应均值mu[i]和方差sigma2[i]。修改要点如下5.1 在HMM类中扩展_gaussian_emission_log方法def _gaussian_emission_log(self, state_i, obs_val): 计算状态i发射连续观测obs_val的log概率 return ( -0.5 * np.log(2 * np.pi * self.sigma2[state_i]) - 0.5 * ((obs_val - self.mu[state_i]) ** 2) / self.sigma2[state_i] )注意self.mu和self.sigma2需在__init__中初始化为np.random.randn(n_states)和np.random.rand(n_states)。5.2 改写forward_log中的观测概率项原离散版log_alpha[t][i] np.log(self.B[i][O[t]])改为连续版log_alpha[t][i] self._gaussian_emission_log(i, O[t])提示O[t]此时为float而非int需确保输入数据为np.float64。哈工大实验说明书中要求“连续观测需先标准化”即O (O - O.mean()) / O.std()否则mu/sigma2更新会发散。5.3 EM算法中M步的高斯参数更新公式在EM_algorithm.py的m_step中新增# 更新高斯均值mu[i] numerator sum( gamma[t][i] * O[t] for t in range(T) ) denominator sum(gamma[t][i] for t in range(T)) self.hmm.mu[i] numerator / denominator # 更新方差sigma2[i] numerator sum( gamma[t][i] * (O[t] - self.hmm.mu[i]) ** 2 for t in range(T) ) self.hmm.sigma2[i] numerator / denominator此处gamma[t][i]仍来自log-space的log_gamma但sum()操作在概率空间进行需先np.exp(log_gamma[t][i])。哈工大保研面试曾以此题考察“如何在log-space中安全执行加权平均”。本文还有配套的精品资源点击获取