在科学和工程领域，理解和预测动态系统的行为至关重要。然而，很多时候我们并不完全知道描述这些系统的精确数学模型。近年来，一种名为SINDy（Sparse Identification of Nonlinear Dynamics，非线性动力学的稀疏识别）的系统辨识方法受到了广泛关注，它能够从数据中自动发现动态系统的数学模型。 SINDy的基本原理 SINDy方法的核心思想是利用稀疏回归技术从测量数据中识别出动态系统的非线性微分方程。简单来说，它假设系统的动态可以由一个函数库中的少数几个函数的线性组合来描述。通过优化算法，SINDy能够选择出对系统动态影响最大的几个函数，并确定它们的系数，从而构建出一个简洁而精确的数学模型。 一个简单的例子 假设我们有一个简单的一维动态系统，其动态方程为： 为了模拟这个系统，我们生成一些模拟数据，并使用SINDy方法来恢复这个方程。 步骤1：生成模拟数据 首先，我们使用北太天元来生成模拟数据。假设初始条件为(x(0) = 1)，我们使用欧拉方法或更高级的数值方法来求解这个微分方程，并收集一段时间内的数据。
% 模拟数据生成 dt = 0.1; % 时间步长 t = 0:dt:2; % 时间向量 t = t'; x = zeros(size(t)); % 初始化状态向量 x(1) = 0.6; % 初始条件 for i = 1:length(t)-1 x(i+1) = x(i) + dt*(-2*x(i) + x(i).^3 ); % 欧拉方法更新状态 end dxdt = diff(x) ./diff(t); % 计算速度步骤2：构建函数库 接下来，我们构建一个函数库，其中包含可能描述系统动态的候选函数。在这个例子中，我们可以选择多项式函数作为候选函数，例如(x), (x^2), (x^3), ...等。
% 构建函数库 x = x(1:end-1); Theta = [x, x.^2, x.^3]; % 候选函数库步骤3：稀疏回归 最后，我们使用稀疏回归方法来确定哪些函数对系统动态有显著影响。在这个例子中，我们可以使用LASSO（Least Absolute Shrinkage and Selection Operator）回归，它是一种能够产生稀疏解的线性回归方法。 针对这个具体的例子，LASSO回归的优化问题可以表示为： 使用迭代软阈值算法(ISTA)求解LASSO问题，其数学原理可以详细阐述如下： LASSO问题的定义 LASSO（Least Absolute Shrinkage and Selection Operator）是一种回归分析方法，它通过构造一个惩罚函数来得到一个较为精炼的模型。LASSO的基本思想是在回归系数的绝对值之和小于一个常数的约束条件下，使残差平方和最小化，从而能够产生某些严格等于0的回归系数，得到可以解释的模型。其数学形式通常表示为最小化以下目标函数： 其中，A 是设计矩阵，b 是观测向量，\beta 是待求的回归系数向量，λ 是正则化参数，控制惩罚项的影响。 2. 迭代软阈值算法(ISTA)的原理 迭代软阈值算法(ISTA)是用于求解LASSO问题的一种有效方法。它通过迭代的方式逐步逼近最优解。算法的基本步骤如下： 初始化 * 设定初始估计值 ，通常可以设为0向量或随机向量。 * 选择合适的步长 t_k 和正则化参数 λ 。 迭代过程 对于每一次迭代 ，执行以下步骤： 1. 梯度下降：首先，对当前估计值 进行梯度下降，以减小残差平方和的部分。这可以通过计算 并按一定步长 更新 β 来实现。 2. 软阈值处理：然后，对梯度下降后的结果进行软阈值处理，以实现L1正则化的效果。软阈值函数的形式为： 其中，\(\text{sgn}(x)\) 是符号函数，返回 \(x\) 的符号。这个软阈值函数会将绝对值小于 λ 的系数缩减到0，而对大于λ 的系数进行一定程度的缩减。 3. 更新估计值：将软阈值处理后的结果作为下一次迭代的起始点，即 。 终止条件 * 当 小于某个预设的阈值时，或者达到最大迭代次数时，算法停止。 3. 算法收敛性 在一定的条件下，ISTA算法可以保证收敛到LASSO问题的最优解。收敛速度取决于步长 的选择和问题的具体性质。通常，为了加速收敛，可以采用一些优化策略，如Nesterov加速或采用更复杂的步长选择策略。这里简单实现一个线性搜索确定步长方法
function [alpha] = linesearch(Theta, y, beta, grad, lambda) % 线性搜索寻找最佳步长 % 初始化参数 alpha = 1; % 初始步长 c = 1e-4; % Armijo条件中的常数 rho = 0.5; % 步长缩减因子 max_ls_iter = 50; % 线性搜索的最大迭代次数 % 计算梯度下降的初始方向 d = -grad; for ls_iter = 1:max_ls_iter % 尝试新的步长 beta_new = beta + alpha * d; % 应用软阈值操作 beta_new = sign(beta_new) .* max(0, abs(beta_new) - lambda); % 计算目标函数的值 f_new = 0.5 * norm(Theta * beta_new - y)^2 + lambda * norm(beta_new, 1); f_old = 0.5 * norm(Theta * beta - y)^2 + lambda * norm(beta, 1); % 检查Armijo条件是否满足 if f_new <= f_old + c * alpha * grad' * d break; % 如果满足，则退出线性搜索 end % 如果不满足，则缩减步长并继续搜索 alpha = rho * alpha; end % 如果达到最大迭代次数仍未找到满足条件的步长，则可能需要调整参数或检查算法实现 if ls_iter == max_ls_iter warning('线性搜索达到最大迭代次数，未找到满足条件的步长。'); end end北太天元中暂时没有内置的LASSO函数，可以用m函数实现一个简单的版本。
function [beta] = simpleLassoSolver(Theta, y, lambda, maxIter, tol) % simpleLassoSolver: 使用迭代软阈值算法(ISTA)求解LASSO问题 % 输入: % Theta - 设计矩阵 (n x p) % y - 响应向量 (n x 1) % lambda - LASSO正则化参数 % maxIter - 最大迭代次数 % tol - 收敛容忍度 % 输出: % beta - 回归系数 (p x 1) % 初始化 [n, p] = size(Theta); beta = zeros(p, 1); % 初始化解 x_old = zeros(p, 1); % 上一步的迭代解 for iter = 1:maxIter % 梯度下降步骤 grad = Theta' * (Theta * beta - y) ; % 使用线性搜索确定最佳步长 alpha = linesearch(Theta, y, beta, grad, lambda); % 更新解 x_new = beta - alpha * grad ; % 应用软阈值操作 beta = sign(x_new) .* max(0, abs(x_new) - lambda); % 检查收敛性 fprintf(['iter = ', num2str(iter), ' beta = ', num2str(beta(1)) ', ' num2str(beta(2)), ', ', num2str(beta(3)), ' norm = ' num2str( norm( Theta*beta - y) ), newline]); if norm( beta -x_old ) < tol && norm ( Theta*beta - y) < tol break; end x_old = beta; end end
% 稀疏回归 % 注意：在实际应用中，你可能需要对数据进行标准化处理 % LASSO回归 lambda = 1e-9; maxIter = 100000;tol = 1e-9; [coef] = simpleLassoSolver(Theta, dxdt, lambda, maxIter, tol) norm( Theta*coef - dxdt)结果分析 通过检查回归系数coef，我们可以发现哪些函数对系统动态有显著影响。在这个例子中，我们希望恢复出原始的微分方程(\dot{x} = -2x + x^3)，因此我们应该期望在回归系数中看到-2和1对应于(x)和(x^3)的系数，而其他函数的系数应该接近于0。 结论 SINDy是一种强大的方法，能够从数据中自动发现动态系统的数学模型。通过结合稀疏回归技术和适当的函数库，SINDy可以准确地识别出描述系统动态的关键函数和它们的系数。这种方法在科学和工程领域具有广泛的应用前景，特别是在对复杂系统进行建模和分析时。北太天元代码汇总
% 模拟数据生成 dt = 0.1; % 时间步长 t = 0:dt:2; % 时间向量 t = t'; x = zeros(size(t)); % 初始化状态向量 x(1) = 0.6; % 初始条件 for i = 1:length(t)-1 x(i+1) = x(i) + dt*(-2*x(i) + x(i).^3 ); % 欧拉方法更新状态 end dxdt = diff(x) ./diff(t); % 计算速度% 构建函数库 x = x(1:end-1); Theta = [x, x.^2, x.^3]; % 候选函数库% 稀疏回归 % 注意：在实际应用中，你可能需要对数据进行标准化处理 % LASSO回归 lambda = 1e-9; maxIter = 100000;tol = 1e-9; [coef] = simpleLassoSolver(Theta, dxdt, lambda, maxIter, tol) norm( Theta*coef - dxdt)function [beta] = simpleLassoSolver(Theta, y, lambda, maxIter, tol) % simpleLassoSolver: 使用迭代软阈值算法(ISTA)求解LASSO问题 % 输入: % Theta - 设计矩阵 (n x p) % y - 响应向量 (n x 1) % lambda - LASSO正则化参数 % maxIter - 最大迭代次数 % tol - 收敛容忍度 % 输出: % beta - 回归系数 (p x 1) % 初始化 [n, p] = size(Theta); beta = zeros(p, 1); % 初始化解 x_old = zeros(p, 1); % 上一步的迭代解 for iter = 1:maxIter % 梯度下降步骤 grad = Theta' * (Theta * beta - y) ; % 使用线性搜索确定最佳步长 alpha = linesearch(Theta, y, beta, grad, lambda); % 更新解 x_new = beta - alpha * grad ; % 应用软阈值操作 beta = sign(x_new) .* max(0, abs(x_new) - lambda); % 检查收敛性 fprintf(['iter = ', num2str(iter), ' beta = ', num2str(beta(1)) ', ' num2str(beta(2)), ', ', num2str(beta(3)), ' norm = ' num2str( norm( Theta*beta - y) ), newline]); if norm( beta -x_old ) < tol && norm ( Theta*beta - y) < tol break; end x_old = beta; end endfunction [alpha] = linesearch(Theta, y, beta, grad, lambda) % 线性搜索寻找最佳步长 % 初始化参数 alpha = 1; % 初始步长 c = 1e-4; % Armijo条件中的常数 rho = 0.5; % 步长缩减因子 max_ls_iter = 50; % 线性搜索的最大迭代次数 % 计算梯度下降的初始方向 d = -grad; for ls_iter = 1:max_ls_iter % 尝试新的步长 beta_new = beta + alpha * d; % 应用软阈值操作 beta_new = sign(beta_new) .* max(0, abs(beta_new) - lambda); % 计算目标函数的值 f_new = 0.5 * norm(Theta * beta_new - y)^2 + lambda * norm(beta_new, 1); f_old = 0.5 * norm(Theta * beta - y)^2 + lambda * norm(beta, 1); % 检查Armijo条件是否满足 if f_new <= f_old + c * alpha * grad' * d break; % 如果满足，则退出线性搜索 end % 如果不满足，则缩减步长并继续搜索 alpha = rho * alpha; end % 如果达到最大迭代次数仍未找到满足条件的步长，则可能需要调整参数或检查算法实现 if ls_iter == max_ls_iter warning('线性搜索达到最大迭代次数，未找到满足条件的步长。'); end end
以下内容为卢朓老师在B站进行的分享，可以作为学习参考： DSGE，即动态随机一般均衡（Dynamic Stochastic General Equilibrium）模型，是宏观经济学中用于解释不同经济变量之间关系的一种数学模型。Dynare是一个为MATLAB设计的工具箱，专门用于动态随机一般均衡（DSGE）模型的估计、模拟和政策分析。这个工具箱极大地简化了DSGE模型的编程和计算过程，使得研究人员和政策制定者能够更方便地构建、分析和理解复杂的宏观经济模型。Dynare工具箱的主要功能包括： 1. 模型设定：Dynare允许用户以简洁的符号形式定义DSGE模型，无需编写大量的底层代码。这使得模型设定更为直观和高效。 2. 参数估计：Dynare提供了多种参数估计方法，如极大似然估计、贝叶斯估计等，用于估计模型参数并评估模型的拟合优度。 3. 模拟分析：用户可以利用Dynare进行脉冲响应分析、方差分解等模拟实验，以研究模型在不同政策冲击下的动态响应。 4. 政策分析：Dynare还提供了丰富的政策分析工具，如政策函数、福利分析等，帮助用户评估不同经济政策的效果和成本。 Dynare工具箱与MATLAB的紧密结合使得用户可以充分利用MATLAB强大的数值计算和图形显示功能，从而更加深入地理解和分析DSGE模型。此外，Dynare还提供了丰富的文档和教程，帮助用户快速上手并充分利用工具箱的各项功能。 在《数值方法：原理、算法及应用》的课上，学生说要为北太天元开发一个具有类似Dynare功能的工具箱。开发这样的工具箱需要深入理解DSGE模型的原理和方法，以及北太天元编程技术。然而，一旦成功开发，这样的工具箱将为研究人员和政策制定者提供一个新的、可能更加灵活和适应性的工具，用于宏观经济建模和分析。 为了配合学生工作， 我自己也先了解了一下DSGE模型， 整理了一个简介来抛砖引玉。 以下是我整理的对DSGE模型的简单介绍： 一、基本定义与特点DSGE模型被视为宏观经济学的动态优化模型，它描述了家庭、企业和政府在一个稳定状态下进行决策的相互作用。该模型通过使用微积分和统计学的工具来建立一个精细的数学模型，以更好地描述各经济主体的组织和决策行为。DSGE模型假设所有人都是理性的，并在不同的环境中作出最优决策。它具有预测能力，可以用来预测在不同的条件下经济变量的变化，例如预测某种政策措施会对通货膨胀率、失业率和GDP等经济变量产生的影响。 DSGE模型的数学表达相当复杂，因为它涉及多个经济主体、市场以及均衡条件。这里提供一个简化的概述，以帮助大家理解模型的基本结构和关键概念。 变量时间：(t)消费：(C_t)劳动供给：(L_t)产出：(Y_t)产品价格：(P_t)工资率：(W_t)投资：(I_t)政府支出：(G_t)税收：(T_t)债券发行：(B_t)利率：(r_t)全要素生产率：(A_t) 目标函数家庭部门：最大化终身效用 其中， 是家庭的效用函数， β 是折现因子。企业部门：最大化利润 其中， 均衡条件商品市场均衡：总供给等于总需求劳动力市场均衡：劳动供给等于劳动需求（通过工资调整实现） 金融市场均衡：资金供给等于资金需求（通过利率调整实现） 政府预算约束：政府支出等于税收和债券发行之和 求解目标 在DSGE模型中，我们的目标是找到一组变量值，这组变量值在满足所有均衡条件和预算约束（商品市场均衡、劳动力市场均衡、金融市场均衡、政府预算约束）的同时，使得所有经济主体的目标函数（家庭效用和企业利润）达到在给定约束下的最大可能值。这通常涉及使用高级数学和统计学方法（如动态一般均衡分析、最优化理论）来求解模型的均衡解。这些均衡解提供了对经济系统在不同条件下运行状态的理解，并可用于政策分析和经济预测。 二、模型构成 从模型的基本构成和关键方程入手描述： 1. 家庭部门 效用最大化问题：家庭部门的目标是最大化其终身效用。设 U 为家庭的效用函数，C_t为消费，L_t 为劳动供给，则家庭的优化问题可以表示为： 其中，β 是折现因子。 在DSGE模型中，C_t 代表时期 t 的消费，它是以货币为单位来衡量家庭在该时期所消费的商品和服务的总价值，而不是以时间来计算。这个变量反映了代表性家庭或整个经济体在该时期的消费水平，是模型中重要的决策变量之一，用于分析家庭部门的消费行为和其对经济的影响。 劳动供给 L_t 在经济学模型中，特别是在DSGE（动态随机一般均衡）模型中，代表在时期 t 劳动者愿意并且能够提供的劳动量。这个变量是衡量劳动力市场状况的关键指标，它反映了劳动力资源的配置情况。可以给出如下的定义: 劳动供给 L_t 是指在时间 t（通常是一个特定的时间段，如一个季度或一年）内，劳动者愿意提供的劳动时间或劳动量。在DSGE模型的上下文中，当提到劳动供给 L_t 时，通常这个变量代表的是模型中代表性家庭（representative household）或代表性劳动者（representative worker）提供的劳动量。在宏观经济模型中，为了简化分析，经常会引入“代表性”的概念，即用一个代表性的经济主体来代表整个经济中的相似主体。所以，这里的 L_t 并不是指所有家庭的所有就业人员的总劳动供给，而是指模型中设定的那个“代表性家庭”或“代表性劳动者”在时期 t 提供的劳动量。这个代表性家庭被假定为能够反映整个经济中家庭或劳动者的平均或典型行为。 劳动供给受到多种因素的影响，包括但不限于工资水平、工作条件、非劳动收入、家庭责任、教育水平以及社会和文化因素。例如，工资水平的提高可能会激励劳动者提供更多的劳动，而恶劣的工作条件或家庭责任可能减少劳动供给。 在DSGE模型中，劳动供给是一个重要的决策变量。家庭部门作为劳动力的提供者，在模型中会基于效用最大化的原则来决定提供多少劳动。企业的生产函数也依赖于劳动供给，因为它直接影响产出水平。 由于DSGE模型强调经济的动态性质，劳动供给 L_t 也会随时间变化。这种变化可能受到经济周期、政策调整、技术进步等多种因素的影响。 在劳动力市场均衡的条件下，劳动供给将等于劳动需求。如果劳动供给大于需求，可能会导致失业率上升；如果劳动供给小于需求，则可能会导致工资上涨和劳动力短缺。 家庭的消费、投资和储蓄决策受到预算的限制。设 W_t 为工资收入，R_t 为资产收益，则预算约束可以表示为： 其中，C_t 是家庭在时期 t 的消费支出，I_t 是家庭在时期 t 的投资支出，而 S_t 代表家庭在时期t 的储蓄。 这个方程的含义是：- C_t：家庭在时期 t 的消费支出，涵盖购买商品和服务等费用。- I_t：家庭在时期 t 的投资支出，可能涉及购买股票、债券、房地产等能够产生未来收益的资产。- S_t：家庭在时期 \(t\) 的储蓄，即家庭选择不消费也不投资，留待未来使用的资金。- W_t：家庭在时期 \(t\) 的工资收入，这是家庭成员通过劳动所获得的报酬。- R_t：家庭在时期 \(t\) 从其所持有的资产中获得的收益，例如股息、利息、租金等。 这个方程揭示了，在任意时期 t，家庭的消费、投资和储蓄的总和必须等于其总收入（包括工资收入和资产收益）。这是家庭经济决策中的一个基本原则，同时也是宏观经济模型中一个关键的约束条件。 这个扩展的预算约束方程是DSGE模型中家庭部门优化问题的重要组成部分。家庭需要在满足这一预算约束的前提下最大化其终身效用。通过求解这个优化问题，我们能够得出家庭在不同经济环境下的最优消费、投资和储蓄决策。 为了简化分析，有时我们会假设 ，即家庭在时期 t 不进行储蓄。在这种情况下，预算约束简化为原来的形式：。然而，在实际的经济模型中，考虑储蓄 \(S_t\) 对于全面理解家庭的经济行为是非常重要的。 2. 企业部门 利润最大化问题：企业部门的目标是最大化利润。设Y_t 为产出，P_t 为产品价格，则利润 为 . 企业利润最大化的优化问题可以表示为： 在这个企业利润最大化的优化问题中， 表示的是企业试图找到一组变量值，这组变量值能够使得目标函数（这里是总利润达到最大值)。变量通常包括企业在每个时期t 的生产决策，具体来说可能包括：产出水平 (Y_t)：企业可以选择在每个时期生产多少产品来最大化利润。产出水平直接影响销售收入 以及可能的成本（包括变动成本，它随产出增加而增加）。企业的产出通常取决于劳动和资本等生产要素。一个简单的生产函数形式可以是：, 其中，K_t 是资本存量。劳动投入量 L_t：企业可以决定雇佣多少劳动力来进行生产。劳动投入量会影响劳动力成本 和生产效率。W_t 是单位劳动成本（例如，每小时的工资率），而L_t代表在时期t 所雇佣的劳动量（可以是工作小时数或员工数量等）。其他可控成本：虽然在这个表达式中用“其他成本”概括了，但在更详细的模型中，企业可能还可以控制如原材料采购量、广告投入、研发投入等，这些都是可以优化的变量。价格 P_t：在某些市场结构中，企业可能有一定的定价权，因此价格也可能成为优化问题中的一个变量。然而，在完全竞争的市场中，企业是价格的接受者，不能控制市场价格。资本投入和其他生产要素：根据模型的复杂程度，企业可能还需要决定资本投入量、技术选择等。 在求解这个优化问题时，企业需要找到这些变量的最优组合，以使得无限期内的总利润最大化。这通常涉及到对市场需求、成本结构、竞争环境等的深入分析和预测。 在这个特定的表达式中， 主要是对产出 (Y_t) 和劳动投入 (L_t)（以及其他可能的可控变量）来取的，因为这些是企业可以直接控制的决策变量。价格 (P_t) 在某些情况下也可能是优化变量，但在完全竞争市场中则不是。 除了劳动力成本外，企业还有其他多种成本，如原材料采购、设备维护、租金、利息支付、市场营销费用等，这些在这里被统称为“其他成本”。 因此，企业的优化问题可以表述为：寻找一个产出路径 ，使得从当前时期到无限远的未来，企业能够获得的总利润 最大化。这个总利润是企业各期收入 与相应时期的总成本（包括劳动力成本 和其他所有成本）之差的累计和。 这个问题是一个动态优化问题，因为它涉及到跨时期的决策，并且需要考虑未来可能的变化和不确定性。在实际操作中，企业会基于市场预测、成本分析和竞争策略来制定生产计划，以实现这个优化目标。 3. 政府部门 政府预算约束：政府的支出（如基础设施建设、社会福利等）通常通过税收和发行债券来融资。政府预算约束可以表示为： 其中，G_t 是政府支出，T_t 是税收收入，B_t是债券发行。G_t：表示政府在时期t的总支出。这包括了基础设施建设、社会福利支出、公共服务支出等所有政府需要承担的费用。T_t：表示政府在时期t通过税收获得的收入。税收是政府主要的财政收入来源之一，用于支持政府的各项开支。B_t：表示政府在时期t通过发行债券筹集的资金。发行债券是政府弥补财政赤字、进行长期投资或应对紧急情况的一种常见融资方式。 这个等式表达了政府预算的基本约束：在任何给定的时期t，政府的总支出必须等于其通过税收和债券发行所获得的总收入。这是确保政府财政可持续性的基本原则。 如果 G_t 大于 T_t + B_t，则意味着政府面临财政赤字，需要通过其他方式（如削减支出、增加税收或进一步举债）来平衡预算。相反，如果 G_t 小于 T_t + B_t，则政府会有财政盈余，可以用于偿还债务、增加储备或进行其他投资。 政府预算约束是宏观经济政策制定中的重要考虑因素，它要求政府在制定支出计划时必须考虑其财政收入状况，以确保财政的稳定和可持续性。 4. 市场均衡条件 在经济学中，市场均衡指的是市场上供求双方达到平衡的状态，即市场上的总供给与总需求相等。以下是您提到的三种市场均衡的详细解释： 商品市场均衡 在商品市场中，均衡意味着商品和服务的总供给 Y_t 等于商品和服务的总需求。总需求通常由三部分组成：消费 C_t、投资 I_t 和政府支出 G_t。因此，商品市场的均衡条件可以表示为：Y_t = C_t + I_t + G_t 这个等式表明，在均衡状态下，经济中生产的所有商品和服务（总供给 Y_t）必须正好等于消费者、投资和政府愿意购买的数量（总需求 C_t + I_t + G_t）。 劳动力市场均衡 劳动力市场的均衡是指劳动供给与劳动需求之间的平衡。这通常是通过工资水平的调整来实现的。当劳动力需求大于供给时，工资水平可能会上升以吸引更多的劳动力进入市场；反之，当劳动力供给过剩时，工资可能会下降。在均衡状态下，愿意并能够工作的劳动力数量正好等于企业愿意雇佣的劳动力数量。 设L_d为劳动力需求，L_s为劳动力供给，W为工资水平。劳动力市场均衡可以表达为： 即，在某一工资水平W下，劳动力的需求量L_d等于劳动力的供给量L_s。 金融市场均衡 金融市场的均衡，特别是资产市场如债券市场的均衡，是确保资金有效分配的关键。在金融市场上，资金的供给者（储蓄者）与资金的需求者（借款者）通过交易各种金融工具（如股票、债券等）来达到均衡。当金融市场的供给和需求相等时，就达到了均衡状态。这意味着储蓄被有效地转化为投资，支持了实体经济的发展。金融市场的均衡对于维护经济的稳定和持续增长至关重要。 M_s为资金供给，M_d为资金需求，r为利率水平（或其他金融资产的价格）。金融市场均衡可以表达为： 即，在某一利率水平r下，资金的供给量M_s等于资金的需求量M_d。 总的来说，这三种市场均衡是宏观经济稳定运行的基础。当这些市场达到均衡时，资源能够得到合理配置，经济能够实现稳定增长。 5. 随机冲击和外部因素 DSGE模型还考虑了随机冲击，如全要素生产率冲击、货币政策冲击等。这些冲击可以通过在模型中加入随机项来体现，例如： 其中，A_t 是全要素生产率， 是其平均值， 是随机冲击项。 随机冲击的概念： - 在DSGE模型中，随机冲击是指那些不可预测的、能够对经济系统产生显著影响的事件或因素。这些冲击可能来源于多个方面，如全要素生产率冲击、货币政策冲击等。 全要素生产率冲击： - 全要素生产率（TFP）是衡量单位总投入的总产量的生产率指标，即总产量与全部要素投入量之比。在DSGE模型中，全要素生产率冲击是指由于技术进步、管理效率提升等因素导致的生产率变化。这种冲击可以通过在模型中加入随机项来体现，如公式。 货币政策冲击： - 货币政策冲击是指中央银行通过调整利率、货币供应量等手段对经济产生的影响。在DSGE模型中，这种冲击被视为一种重要的随机因素，因为它能够显著改变经济的整体状况。模型通过引入货币政策变量及其随机变动来反映这种冲击。 随机冲击在模型中的作用： - DSGE模型通过引入随机冲击项，能够更真实地模拟现实经济的运行状况。这些冲击不仅影响经济的短期波动，还可能对长期经济增长趋势产生影响。因此，在模型中考虑这些冲击对于准确预测和制定经济政策具有重要意义。 三、DSGE模型的演化与分类DSGE模型在过去40年中在宏观经济建模领域取得了重要进展。按照发展阶段和模型设定的不同，DSGE模型大致可分为真实经济周期（RBC）模型和新凯恩斯主义（NK）模型。NK模型与RBC模型的主要区别在于是否包含垄断竞争和名义刚性或黏性。 四、应用与地位DSGE模型目前在宏观经济学研究中占重要地位，甚至可以说是主导地位。它主要用于讨论经济增长、经济周期以及政策工具（财政和货币政策）的效果。该模型能够解释经济中的不确定性，并继承了数理经济学的一般均衡理论。 - DSGE模型已被广泛应用于宏观经济政策分析、经济预测以及经济周期研究等领域。通过引入随机冲击，模型能够更全面地反映现实经济的复杂性和不确定性，从而为政策制定者提供更准确的决策依据。同时，DSGE模型还有助于我们深入理解经济体系中各种因素之间的相互作用和影响机制。 综上所述，DSGE模型通过一系列数学方程和约束条件来描述经济中不同部门的行为和市场均衡状态。这些方程可以根据具体的研究问题和数据可用性进行适当的调整和扩展。
以下内容为卢朓老师在B站进行的分享，可以作为学习参考： GARCH族模型，即广义自回归条件异方差模型（Generalized Autoregressive Conditional Heteroskedasticity Model），是一类用于估计时间序列数据波动率的统计模型。该模型由Bollerslev在1986年提出，是对ARCH（自回归条件异方差）模型的一种重要扩展。GARCH模型在金融时间序列分析中具有广泛的应用价值，尤其是在金融市场波动性的建模和预测方面。 一、GARCH模型的基本概念 GARCH模型主要用于描述时间序列数据（如股票价格、汇率、利率等）的波动性特征。传统计量经济学假设时间序列变量的波动幅度（方差）是固定的，但这往往不符合实际情况。例如，股票收益的波动幅度通常会随时间变化，表现出聚集性特征，即大的波动后面往往跟着大的波动，小的波动后面往往跟着小的波动。GARCH模型通过引入条件异方差来描述这种波动性聚集现象，从而更准确地捕捉时间序列数据的波动性特征。 二、GARCH模型的结构 GARCH模型通常由两部分组成：均值方程和方差方程。 均值方程：通常是一个ARMA（自回归移动平均）模型，用于描述时间序列数据的线性关系。它表示时间序列数据在某一时刻的期望值，即数据的均值部分。 方差方程：是GARCH模型的核心，用于描述时间序列数据的波动性。方差方程是一个自回归移动平均模型，但作用于时间序列的方差上，而不是直接作用于时间序列数据本身。通过考虑过去的波动率和误差项，方差方程能够预测未来的波动率。 三、GARCH模型的数学表达式 GARCH模型的均值方程通常用于描述时间序列数据的条件均值，即数据在给定信息集下的期望值。然而，值得注意的是，GARCH模型的核心在于其方差方程，该方程用于刻画时间序列数据的条件异方差性，即波动性。尽管如此，在构建GARCH模型时，通常也会同时指定一个均值方程。 对于GARCH模型的均值方程，其形式可以相对简单，也可以相对复杂，具体取决于数据的特性和研究目的。一个常见的均值方程认为时间序列数据的均值是恒定的。这种方程形式适用于数据围绕某个固定水平波动的情况。 其中，y_t 是t时刻的观测值， \mu 是常数均值，\epsilon_t 是残差项。 在GARCH模型中，残差项 $\epsilon_t$ 通常被表示为条件方差 $\sigma_t$ 和一个独立同分布（iid）的随机变量 $z_t$ 的乘积，即 这里，$z_t$ 通常被假设为标准正态分布 $N(0,1)$ 的随机变量，意味着它有一个均值为0和方差为1的正态分布。 条件方差 $\sigma_t$ 是通过GARCH模型的方差方程来估计的，一般形式的GARCH(p,q)模型的方差方程可以表示为： 其中， 是t时刻的条件方差，ϵ_t 是t时刻的残差项，α_0 是常数项，α_i 和 β_j 是模型的参数。p和q分别表示方差方程中自回归项和移动平均项的阶数。 四、GARCH族模型的求解 GARCH族模型作为一类用于精确估计金融时间序列数据波动率的统计模型，其求解过程涉及多个复杂步骤。 首先，准备和预处理金融时间序列的历史数据是至关重要的。这一步骤包括数据的清洗、转换和差分等，以确保数据的准确性、一致性和适用性。通过预处理，我们可以消除数据中的异常值、填补缺失值，并将其转换为适合模型分析的形式。 接下来，根据数据的特性和研究需求，我们需要选择合适的GARCH模型，并明确设定均值方程和方差方程的形式。这一步骤要求我们具备扎实的统计理论知识，以便能够准确地描述时间序列数据的动态特征。 在模型设定之后，构建似然函数成为关键。似然函数是基于残差项的分布（如标准正态分布）来构建的，它反映了模型参数与观测数据之间的匹配程度。通过最大化似然函数，我们可以找到最优的参数估计值，从而使模型更好地拟合观测数据。 参数估计是求解GARCH族模型的核心步骤。我们需要运用极大似然估计法（MLE）等高级统计参数估计方法，并通过求解似然函数的导数等于零的方程组或使用迭代算法（如Newton-Raphson算法、BFGS算法等）来找到最优的参数估计值。这一步骤要求我们具备深厚的数学功底和编程能力，以便能够准确地估计出模型的参数。 最后，对估计得到的模型进行严格的检验是必不可少的。我们需要运用残差检验、模型验证等统计方法，评估模型的拟合效果和预测能力。如果模型通过检验，我们就可以使用它来进行时间序列数据的波动率预测和进一步的分析。 五、GARCH族模型的衍生与发展 随着研究的深入，GARCH模型得到了不断的扩展和完善，形成了GARCH族模型。这些衍生模型包括但不限于： EGARCH：指数GARCH模型，用于解决GARCH模型中对正负扰动的对称性问题。 GJR-GARCH：Glosten-Jagannathan-Runkle GARCH模型，同样用于捕捉正负扰动对波动率的不对称影响。 APARCH：平均绝对偏差GARCH模型，通过引入绝对值项来改进模型。 IGARCH：积分GARCH模型，用于处理条件方差无限持久的情况。 GARCH-M：均值GARCH模型，在均值方程中引入方差项，以反映风险与收益之间的关系。 六、GARCH模型的应用 GARCH模型在金融市场中具有广泛的应用，主要包括： 波动性预测：通过对历史数据的分析，GARCH模型可以预测未来时间序列数据的波动性，为投资者提供决策支持。 风险管理：金融机构可以利用GARCH模型进行风险定价和风险管理，提高经营效率。 投资组合优化：投资者可以根据GARCH模型的预测结果调整投资组合，以降低投资风险并提高收益。 综上所述，GARCH族模型是一类强大的时间序列分析工具，特别适用于金融市场的波动性建模和预测。通过不断的研究和发展，GARCH模型在金融领域的应用前景将更加广阔。 七、北太天元的代码示例
% 示例数据（通常这里应该是金融时间序列数据，如股票收益率）
N = 1000; %数据点的总数
data = randn(N,1); % 使用随机数作为示例
% 设定初始值
mu = mean(data); % 样本均值
sigma2 = zeros(N,1); %和data 一样为列向量
sigma2(1) = var(data); % 样本方差
residuals = data - mu; % 计算残差
% 迭代计算GARCH(1,1)模型的条件方差（这里我们先使用简单的样本方差作为初始值）
omega = 0.01;
alpha = 0.1;
beta = 0.85;
for t = 2:length(data)
sigma2(t) = omega + alpha * residuals(t-1)^2 + beta * sigma2(t-1);
end
% 使用极大似然估计法（MLE）来估计参数
params0 = [omega, alpha, beta]; % 初始参数
options = optimset('Display','iter','TolFun',1e-8);
params_est = fminsearch(@(params) -logLikelihoodGARCH(params, residuals, sigma2), params0, options);
% 输出估计的参数
omega_est = params_est(1);
alpha_est = params_est(2);
beta_est = params_est(3);
fprintf('Estimated omega: %f\n', omega_est);
fprintf('Estimated alpha: %f\n', alpha_est);
fprintf('Estimated beta: %f\n', beta_est);
% 辅助函数：计算GARCH(1,1)模型的对数似然函数
function ll = logLikelihoodGARCH(params, residuals, sigma2)
omega = params(1);
alpha = params(2);
beta = params(3);
% 这里我们不重新计算残差和条件方差，而是使用传入的参数
% 计算对数似然函数
T = length(residuals);
ll = -0.5 * sum(log(2*pi*sigma2(2:T)) + residuals(2:T).^2 ./ sigma2(2:T));
% 注意：我们从第二个观测值开始计算对数似然，因为第一个观测值的条件方差是已知的（通常是样本方差）
end 代码可以copy 到文心一言里，让它给你解读一下， 这里仅仅解释一下对数似然函数。对数似然函数是根据GARCH(1,1)模型的定义和假设来推导的。GARCH(1,1)模型是一个时间序列模型，用于描述金融资产的波动率。它假设时间序列的残差（去均值化后的观测值）服从一个条件正态分布，其条件方差由过去的残差和过去的条件方差共同决定。 具体来说，GARCH(1,1)模型可以表示为： 其中，$\sigma_t^2$ 是时刻t的条件方差，$\omega$、$\alpha$ 和 $\beta$ 是模型的参数，$r_{t-1}$ 是时刻t-1的残差。 现在，我们假设残差 $r_t$ 服从条件正态分布 $N(0, \sigma_t^2)$。这意味着，给定过去的信息，$r_t$ 的概率密度函数是： 对数似然函数是概率密度函数的自然对数，即： 由于我们处理的是一系列观测值 $r_2, \ldots, r_T$（通常第一个观测值用于初始化），对数似然函数是这些观测值的对数概率密度之和： 将概率密度函数代入上式，并取对数，我们得到： 这就是我们在代码中使用的对数似然函数的形式。注意，这里有一个负号，因为我们在最大化对数似然函数时实际上是在最小化负对数似然函数（这是一个常见的做法，因为优化算法通常设计为最小化目标函数）。因此，当我们使用`fminsearch`函数时，我们实际上是在最小化负对数似然函数，这等价于最大化对数似然函数。 七、北太天元V3.5的一个小bug的修复 北太天元报错: 错误(文件 C:\baltamatica\plugins\optimization\scripts\fminsearch.m, 行474, 列1): 与function配对的end一旦出现一个以上，就必须做到一一对应。现在是没有做到一一对应。 C:\baltamatica\plugins\optimization\scripts\fminsearch.m解析出错 修复方法： 在 C:\baltamatica\plugins\optimization\scripts\fminsearch.m的471行后面增加一行，内容是end 增加后， 如下图所示相关视频讲解：
以下内容为卢朓老师在B站进行的分享，可以作为学习参考： 在科学与工程的广阔领域里，优化问题无处不在，它们挑战着我们的智慧和计算能力。其中，旅行商问题（TSP，Traveling Salesman Problem）便是一个经典而富有挑战性的优化难题。在这个问题中，一位旅行商需要访问一系列城市，并找到一条从起点出发，经过所有城市后返回起点的最短路径。随着城市数量的增加，问题的复杂度呈爆炸式增长，寻找最优解变得异常困难。旅行商问题的挑战 旅行商问题的难度在于其计算量的巨大。当城市数量较少时，我们还可以尝试通过穷举所有可能的路径来找到最短路径。但是，随着城市数量的增加，可能的路径数量呈指数级增长，穷举法变得不切实际。例如，当有10个城市时，可能的路径数量就已经达到了3628800条！对于更多的城市，计算量更是惊人。模拟退火算法：大自然的启示 为了应对这一挑战，科学家们开发了各种启发式搜索算法，其中模拟退火算法便是一种非常有效的方法。这种算法借鉴了物理学中固体退火的原理。正如金属在加热后再缓慢冷却的过程中，原子会逐渐排列成有序的结构，从而达到最低的能量状态，模拟退火算法也通过模拟这一过程来寻找问题的最优解。 在算法开始时，我们设定一个较高的“温度”，并随机生成一个可能的解作为当前解。然后，算法会在这个解的“邻域”内随机探索新的解。如果新解比当前解更优（即路径更短），我们就接受这个新解。如果新解不如当前解，我们也有一定的概率接受它，这个概率随着温度的降低而逐渐减小。这个过程就像金属在退火过程中，即使有时原子会排列成能量较高的状态，但随着温度的降低，它们最终还是会趋于最低能量状态。通过不断地探索和优化，模拟退火算法能够逐渐找到越来越接近最优解的解决方案。代码实现：模拟退火求解旅行商问题 接下来，我们将通过一段北太天元代码来展示如何使用模拟退火算法求解旅行商问题。这段代码首先生成了一系列随机分布的城市（或点），并计算了每对城市之间的距离。然后，它使用模拟退火算法来找到一条近似最短路径。
%北太天元 用模拟退火法 求解 旅行商 问题% 生成num_points个随机点的经纬度 num_points = 10; longitudes = rand(num_points, 1) * 360 - 180; % 随机经度在-180到180之间 latitudes = rand(num_points, 1) * 180 - 90; % 随机纬度在-90到90之间 % 将经纬度保存到sj.txt文件中 fileID = fopen('sj.txt', 'w'); for i = 1:num_points fprintf(fileID, '%f\t%f\n', longitudes(i), latitudes(i)); end fclose(fileID); disp('sj.txt文件已生成，包含10个随机点的经纬度。'); % 读取数据 sj0 = readmatrix('sj.txt'); x = sj0(:,[1:2:end]); y = sj0(:,[2:2:end]); sj = [x, y]; d1 = [70, 40]; sj = [d1; sj; d1]; sj = deg2rad(sj); % 初始化距离矩阵 d = zeros(num_points+2); % 计算距离矩阵 for i = 1:num_points+1 for j = i+1:num_points+2 d(i,j) = 6370 * acos(cos(sj(i,1)-sj(j,1)) .* cos(sj(i,2)) .* cos(sj(j,2)) + sin(sj(i,2)) .* sin(sj(j,2))); end end d = d + d'; % 初始化路径和路径长度 path = 1:num_points+2; long = sum(d(sub2ind(size(d),path(1:end-1),path(2:end)))); % 模拟退火算法参数 e = 1e-10; L = 20000; at = 0.8; T = 1; % 模拟退火过程 for k = 1:L for i = 1:num_points temp = path; m = randperm(num_points,2); temp( [m(1)+1,m(2)+1] ) = path([m(2)+1,m(1)+1]); % 交换路径中的两个点 temp_long = sum(d(sub2ind(size(d),temp(1:end-1),temp(2:end)))); if temp_long < long || (long - temp_long) / T > log(rand) path = temp; long = temp_long; end end T = T * at; if T < e break; end end % 输出结果 xx = sj(path,1); yy = sj(path,2); plot(rad2deg(xx), rad2deg(yy), '-o'); xlabel('经度'); ylabel('纬度'); hold on;scatter(rad2deg(xx(1)),rad2deg(yy(1)),'r')for i = 1:length(path) text(rad2deg(xx(i)), rad2deg(yy(i)), sprintf('%d', i), 'Color', 'red'); end hold off;disp(['总路径长度为：', num2str(long), ' km']); disp(['所需时间为：', num2str(long / 1000), ' 小时']); % 假设飞机速度为1000km/h这段代码通过模拟退火算法的不断迭代和优化，最终输出了一条近似最短路径以及总路径长度。尽管由于问题的复杂性和计算量的限制，我们无法保证找到的解是最优的，但模拟退火算法通常能够给出一个非常接近最优解的解决方案。结语 模拟退火算法作为一种启发式搜索算法，在解决旅行商问题等复杂优化问题方面展现出了强大的能力。通过模拟物理退火过程，它能够在庞大的解空间中高效地探索并找到近似最优解。尽管随着城市数量的增加，旅行商问题的计算量仍然是一个巨大的挑战，但模拟退火算法为我们提供了一种有效的应对策略。在未来的研究和应用中，我们可以进一步优化和改进这一算法，以应对更加复杂和庞大的优化问题。 p.s. 这里还假设飞机和目标点都是地球表面, 而且飞机总是贴着地球表面飞行,因还涉及到地球表面上的两个位置点之间的距离(不是三维空间的直线距离,而是贴着球面两点连线之间的最短距离, 也就是测地线的长度), 因此中学生朋友还请参考 北太天元科普: 测地线 - 哔哩哔哩 (bilibili.com) p.s. 上面的代码用文心一言进行解读 这段代码是一个使用模拟退火算法（Simulated Annealing）来解决旅行商问题（Traveling Salesman Problem, TSP）的MATLAB(或者北太天元)脚本。旅行商问题是一个经典的组合优化问题，目标是找到一条经过所有城市并返回起点的最短路径。下面是对代码的详细解读：生成随机点：num_points = 10; 定义了点的数量为10。longitudes 和 latitudes 分别生成了10个随机经度和纬度，经度范围在-180到180之间，纬度范围在-90到90之间。保存点到文件：将生成的经纬度保存到sj.txt文件中，每行一个点，经度在前，纬度在后。读取数据并初始化：读取sj.txt文件中的数据，并将其转换为弧度制（deg2rad），同时将出发点和终点设置为(70, 40)。计算距离矩阵：初始化一个num_points+2维的零矩阵d（包括出发点和终点）。使用球面距离公式计算每对点之间的距离，并填充到距离矩阵d中。绘图显示：使用scatter函数在地图上显示所有的点，并使用text函数在每个点旁边显示其编号（除了出发点和终点）。初始化路径和路径长度：初始化路径为按顺序访问每个点，并计算这个路径的总长度。模拟退火算法：设置模拟退火算法的参数，包括温度T、冷却率at、结束温度e和迭代次数L。在每次迭代中，随机选择两个点交换它们在路径中的位置，计算新路径的长度。如果新路径比当前路径短，或者满足一定的概率条件，则接受新路径。更新温度，如果温度低于结束温度，则结束迭代。输出结果：绘制最终路径的地图，并在每个点旁边显示其编号。显示总路径长度和假设飞机速度为1000km/h时的所需时间。 总的来说，这段代码展示了如何使用模拟退火算法来解决旅行商问题，包括生成随机点、计算距离矩阵、使用模拟退火算法寻找最短路径，并绘制结果地图和显示路径长度和所需时间。
数据age,sex,cp,trestbps,chol,fbs,restecg,thalach,exang,oldpeak,slope,ca,thal,target42,1,1,120,295,0,1,162,0,0,2,0,2,148,0,2,130,275,0,1,139,0,0.2,2,0,2,154,1,2,120,258,0,0,147,0,0.4,1,0,3,154,0,2,108,267,0,0,167,0,0,2,0,2,151,0,2,130,256,0,0,149,0,0.5,2,0,2,157,0,0,128,303,0,0,159,0,0,2,1,2,161,1,0,138,166,0,0,125,1,3.6,1,1,2,045,1,0,115,260,0,0,185,0,0,2,0,2,160,0,3,150,240,0,1,171,0,0.9,2,0,2,149,0,0,130,269,0,1,163,0,0,2,0,2,153,1,2,130,246,1,0,173,0,0,2,3,2,154,1,0,140,239,0,1,160,0,1.2,2,0,2,151,1,0,140,298,0,1,122,1,4.2,1,3,3,049,1,2,118,149,0,0,126,0,0.8,2,3,2,058,1,2,105,240,0,0,154,1,0.6,1,0,3,148,1,1,130,245,0,0,180,0,0.2,1,0,2,148,1,0,122,222,0,0,186,0,0,2,0,2,152,1,2,138,223,0,1,169,0,0,2,4,2,154,1,1,192,283,0,0,195,0,0,2,1,3,063,0,0,150,407,0,0,154,0,4,1,3,3,070,1,1,156,245,0,0,143,0,0,2,0,2,1这是一个关于心脏病筛查的数据集片段，每一行代表一个病人的信息。以下是各列字段的解释：age：年龄 - 病人的年龄（以岁为单位）。sex：性别 - 通常为0（女性）或1（男性）。cp：胸痛类型 - 这是一个分类变量，表示胸痛的类型。不同的数字可能代表不同的胸痛类型（例如，0可能表示无胸痛，1、2、3、4可能代表不同类型的胸痛）。trestbps：静息血压 - 病人在静息状态下的血压（单位可能是mmHg）。chol：血清胆固醇 - 病人的血清胆固醇水平（单位可能是mg/dL）。fbs：空腹血糖 - 病人的空腹血糖水平（0可能表示正常，1可能表示高于正常）。restecg：静息心电图结果 - 这可能是一个分类变量，表示静息状态下的心电图结果。thalach：最大心率 - 病人能达到的最大心率（单位可能是次/分钟）。exang：运动诱发的心绞痛 - 0可能表示没有，1可能表示有。oldpeak：ST段下降 - 这是心电图的一个参数，可能表示心肌缺血的程度（单位可能是mV或其他相关单位）。slope：ST段峰值斜率 - 这也是心电图的一个参数，可能与心脏的病理状况有关（例如，1、2、3可能表示不同的斜率类型）。ca：冠状动脉造影的主要血管数量 - 这可能是一个分类变量，表示在冠状动脉造影中检测到的主要血管数量（0可能表示无血管受损，其他数字可能表示受损血管的数量）。thal：缺陷类型 - 这可能是一个分类变量，表示某种心脏缺陷或病理状况的类型。target：目标变量 - 这是我们要预测或分类的变量。通常，1可能表示病人有心脏病，0可能表示没有。下面是北太天元软件的代码, 关于心脏病数据集的多类支持向量机（SVM）分类器的完整实现，包括数据预处理、主成分分析（PCA）降维、交叉验证划分、模型训练和测试、以及结果可视化。下面是对代码的详细解读：导入插件和数据北太天元代码
load_plugin("optimization");
data = readmatrix('heart disease.csv');load_plugin("optimization"): 加载优化工具箱，用于求解二次规划问题。readmatrix('heart disease.csv'): 读取心脏病数据集。数据预处理北太天元代码if any(any(ismissing(data))) data = rmmissing(data); end检查数据中是否有缺失值，并使用rmmissing函数删除包含缺失值的行。分离特征和标签北太天元代码
X = data(:,1:end-1); y = data(:,end);将数据集中的最后一列作为标签，其余列作为特征。主成分分析降维北太天元代码
[coeff,score,latent] = pca(X); X_pca = score(:,1:2);使用自定义的pca函数对特征进行降维，保留前两个主成分。划分训练集和测试集北太天元代码
cv = my_cvpartition(size(X_pca,1),'HoldOut',0.3); % ...（省略部分代码）使用自定义的my_cvpartition函数进行HoldOut交叉验证划分，将数据集分为70%的训练集和30%的测试集。训练SVM模型北太天元代码
SVMModel = my_fitcecoc(X_train,Y_train,C);使用自定义的my_fitcecoc函数训练多类SVM模型。该函数内部采用一对一（One-vs-One）策略，为每个类别对训练一个SVM分类器。预测测试集并计算准确率北太天元代码
[label,score] = predict(SVMModel,X_test); accuracy = sum(label-1 == Y_test) / length(Y_test);使用predict函数对测试集进行预测，并计算准确率。注意这里假设标签是从1开始的，所以计算准确率时需要将预测的标签减1。可视化matlab复制代码
gscatter(X_pca(:,1),X_pca(:,2),y); % ...（省略绘制决策边界的代码）使用gscatter函数绘制原始数据点的散点图，并使用轮廓线绘制决策边界。自定义函数解读pca: 实现主成分分析算法，对数据进行降维。my_cvpartition: 自定义的交叉验证划分函数，这里只实现了HoldOut方法。my_fitcecoc: 实现多类SVM分类器的训练，采用一对一策略。my_fitcsvm和my_fitcsvm_soft: 实现硬间隔和软间隔SVM分类器的训练，使用二次规划求解。predict: 对测试集进行预测，并返回预测的标签和得分。predictOneVsAll: 使用单个SVM模型进行预测，并返回预测结果和决策函数的值。注意事项代码中存在一些假设和简化，例如标签从1开始、直接计算决策函数的值作为得分等，这些在实际应用中可能需要根据具体情况进行调整。在使用二次规划求解SVM问题时，需要注意目标函数和约束条件的设置，以确保求解的正确性。可视化部分仅绘制了前两个主成分的散点图和决策边界，对于更高维度的数据可能需要进行其他形式的可视化。下面是全部代码
load_plugin("optimization");
% 导入数据
data = readmatrix('heart disease.csv');
% 检查并处理缺失值
if any(any(ismissing(data)))
data = rmmissing(data);
end
% 分离特征和标签
X = data(:,1:end-1);
y = data(:,end);
% 主成分分析降维
[coeff,score,latent] = pca(X);
X_pca = score(:,1:2);
% 划分训练集和测试集
cv = my_cvpartition(size(X_pca,1),'HoldOut',0.3);
idx = cv.test;
train_idx = setdiff(1: size(X_pca,1), idx);
X_train = X_pca(train_idx,:);
Y_train = y(train_idx,:);
X_test = X_pca(idx,:);
Y_test = y(idx,:);
% 训练SVM模型
C = 1e2; % 设置较大的C值以确保硬间隔分类（对于线性可分数据）
SVMModel = my_fitcecoc(X_train,Y_train,C);
% 预测测试集
[label,score] = predict(SVMModel,X_test);
% 计算准确率
accuracy = sum(label-1 == Y_test) / length(Y_test);
% 可视化
gscatter(X_pca(:,1),X_pca(:,2),y);
hold on;
% 绘制决策边界c
d = 0.15;
[x1Grid,x2Grid] = meshgrid(min(X_pca(:,1)):d:max(X_pca(:,1)),...
min(X_pca(:,2)):d:max(X_pca(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];
[~,scores] = predict(SVMModel,xGrid);
contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k');
% 输出准确率
fprintf('The accuracy of the SVM model is %.2f%%\n', accuracy * 100);
hold off;
function [coeff, score, latent] = pca(X)
% X: 数据矩阵，每一列是一个特征，每一行是一个样本
% coeff: 主成分系数
% score: 表示主成分得分
% latent: 主成分对应的特征值
% 标准化数据（均值为0，方差为1）
X = (X - mean(X)) ./ std(X);
% X = (X - mean(X));
% 计算协方差矩阵
CovMat = cov(X);
% 对协方差矩阵进行特征值分解
[V, D] = eig(CovMat);
% 将特征值按降序排序，并获取对应的特征向量
[latent, order] = sort(diag(D), 'descend');
coeff = V(:, order);
% 计算主成分得分
score = X * coeff;
end
function cv = my_cvpartition(n, method, param)
% n: 总样本数
% method: 分区方法，目前仅支持 'HoldOut'
% param: 分区方法的参数，对于 'HoldOut'，该参数为留出的比例
% 初始化输出结构体
cv = struct('train', [], 'test', []);
if strcmp(method, 'HoldOut')
% 生成一个从 1 到 n 的整数数组
idx = 1:n;
% 随机打乱 idx 数组
idx = idx(randperm(n));
% 计算留出样本的数量
numHoldOut = floor(n * param);
% 分配训练和测试索引
cv.test = idx(1:numHoldOut);
cv.train = idx(numHoldOut+1:end);
else
error('Unsupported partition method.');
end
end
function SVMModel = my_fitcecoc(X_train, Y_train, C)
% X_train: 特征矩阵 (n x d)，其中 n 是样本数，d 是特征维度
% Y_train: 标签向量 (n x 1)
% C: 正则化参数，控制对错分样本的惩罚程度
% 获取类别数量
uniqueClasses = unique(Y_train);
numClasses = length(uniqueClasses);
% 初始化 SVM 模型结构体数组
SVMModel = struct('Classifiers', cell(numClasses*(numClasses-1)/2, 1), ...
'ClassPairs', [], 'ClassPairsIndex', [], 'UniqueClasses', uniqueClasses);
% 训练一对一 SVM 分类器
classifierIndex = 1;
for i = 1:numClasses
for j = (i+1):numClasses
% 提取当前类别对的训练数据
classI = Y_train == uniqueClasses(i);
classJ = Y_train == uniqueClasses(j);
X_train_ij = [X_train(classI, :); X_train(classJ, :)];
Y_train_ij = [ones(sum(classI), 1); -1*ones(sum(classJ), 1)];
% 训练 SVM 分类器
SVMModel.Classifiers{classifierIndex} = my_fitcsvm_soft(X_train_ij, Y_train_ij, C);
% 记录当前分类器对应的类别对
SVMModel.ClassPairs(classifierIndex, :) = [uniqueClasses(i), uniqueClasses(j)];
SVMModel.ClassPairsIndex(classifierIndex, :) = [i,j];
classifierIndex = classifierIndex + 1;
end
end
end
function [wb] = my_fitcsvm(X, Y)
% X: 特征矩阵 (n x d)，其中 n 是样本数，d 是特征维度
% Y: 标签向量 (n x 1)，取值为 +1 或 -1
load_plugin("optimization");
% 假设X和Y已经定义，如之前的示例
N = size(X, 1); % 数据点的数量
D = size(X,2); % 数据的维度
% 将w和b组合成一个向量，以便使用quadprog
% 注意：这里我们将w放在前面，b放在最后
p = rand(D + 1, 1); % 初始猜测解（通常设为0）
Aeq = []; % 没有等式约束
beq = [];
% 构造二次规划的目标函数和线性不等式约束
H = eye(D+1); % Hessian矩阵（目标函数的二次项系数）
f = zeros(D+1,1); % 目标函数的一次项系数（对于SVM原始问题，通常为负）
% 线性不等式约束 Ax <= b
% 对于每个数据点 (x_i, y_i)，我们有 y_i * (w' * x_i + b) >= 1
A = [-Y.*X, -Y]; % 不等式约束的系数矩阵
b = -ones(N, 1); % 不等式约束的右侧向量
% 使用quadprog求解二次规划问题
options = optimoptions('quadprog','Algorithm','interior-point');
[w_b, fval, exitflag, output] = quadprog(H, f, A, b, Aeq, beq, [], [], p, options);
% 分离出w和b
w = w_b(1:D);
b = w_b(D+1);
wb = struct('w',w,'b',b);
end
function [wb] = my_fitcsvm_soft(X, Y, C)
% X: 特征矩阵 (n x d)，其中 n 是样本数，d 是特征维度
% Y: 标签向量 (n x 1)，取值为 +1 或 -1
% C: 正则化参数，控制对错分样本的惩罚程度
% 假设X和Y已经定义，如之前的示例
N = size(X, 1); % 数据点的数量
D = size(X, 2); % 数据的维度
% 初始化松弛变量（slack variables）
xi = zeros(N, 1);
% 将w和b以及松弛变量组合成一个向量，以便使用quadprog
p = zeros(D + 1 + N, 1); % 初始猜测解
% 构造二次规划的目标函数和线性不等式约束
H = [diag( [ones(1,D),0] ), zeros(D+1, N); zeros(N, D+1), zeros(N,N)]; % Hessian矩阵
f = [zeros(D+1, 1); C*ones(N, 1)]; % 目标函数的一次项系数
% 线性不等式约束 Ax <= b
% 对于每个数据点 (x_i, y_i)，我们有 y_i * (w' * x_i + b) >= 1 - xi_i
A = [-Y.*X, -Y, diag( -ones(1,N) )]; % 不等式约束的系数矩阵
b = -ones(N, 1); % 不等式约束的右侧向量
Aeq = []; % 没有等式约束
beq = [];
lb = [ -inf*ones(D+1,1); zeros(N,1)];
ub = inf*ones(D+1+N,1);
% 使用quadprog求解二次规划问题
options = optimoptions('quadprog','Algorithm','interior-point');
[w_b_xi, fval, exitflag, output] = quadprog(H, f, A, b, Aeq, beq, lb, ub, p, options)
% 分离出w, b和xi
w = w_b_xi(1:D);
b = w_b_xi(D+1);
xi = w_b_xi(D+2:end);
% 构造并返回结构体wb，包含w和b
wb = struct('w', w, 'b', b);
end
function [label, score] = predict(SVMModel, X_test)
% SVMModel: 训练好的多类SVM模型，由my_fitcecoc函数返回
% X_test: 测试集特征数据
% label: 预测的类别标签
% score: 预测的得分（可选，这里简化为投票得分）
% 初始化预测标签和得分
[numTest, ~] = size(X_test);
numClasses = length(SVMModel.Classifiers);
label = zeros(numTest, 1);
score = zeros(numTest, numClasses);
% 对每个测试样本进行预测
for i = 1:numTest
sample = X_test(i, :);
votes = zeros(1, numClasses+1);
% 使用每个SVM分类器进行预测，并进行投票
for j = 1:numClasses
svmModel = SVMModel.Classifiers{j};
[~, ~, ~, output] = predictOneVsAll(svmModel, sample);
if output == 1
classIdx1 = SVMModel.ClassPairsIndex(j, 1);
votes(classIdx1) = votes(classIdx1) + 1;
else
classIdx2 = SVMModel.ClassPairsIndex(j, 2);
votes(classIdx2) = votes(classIdx2) + 1;
end
end
% 找到得票最多的类别作为预测标签
[~, maxVoteIdx] = max(votes);
label(i) = SVMModel.UniqueClasses(maxVoteIdx);
score(i, maxVoteIdx) = votes(maxVoteIdx); % 将最高票数作为得分
end
end
function [predictedLabel, predictedScore, decisionValues, output] = predictOneVsAll(svmModel, sample)
% 使用单个SVM模型进行预测
% 这里简化为直接计算决策函数的值，实际应用中可能需要更复杂的处理
w = svmModel.w;
b = svmModel.b;
decisionValues = dot(w, sample) + b;
if decisionValues >= 0
predictedLabel = 1; % 属于正类
else
predictedLabel = -1; % 属于负类
end
predictedScore = abs(decisionValues); % 决策函数的绝对值作为得分
output = predictedLabel; % 用于投票的输出
end处理非连续或非标准化标签处理非连续或非标准化标签（例如，不是从1开始的整数或不是连续整数）是机器学习中的一个常见问题。在我的SVM分类器实现中，为了处理这种情况，我引入了SVMModel.UniqueClasses和SVMModel.ClassPairsIndex两个字段。下面我将详细解释这两个字段的作用以及如何在预测过程中使用它们。SVMModel.UniqueClassesSVMModel.UniqueClasses存储了数据集中所有唯一的类别标签。当训练数据中的标签不是从1开始的连续整数时，这个字段特别有用。例如，如果标签是[-1, 0, 2]，那么UniqueClasses就会是[-1, 0, 2]。在训练过程中，这个字段用于记录所有不同的标签，以便在预测时能够正确地将预测的类别索引映射回原始标签。SVMModel.ClassPairsIndexSVMModel.ClassPairsIndex存储了训练一对一SVM分类器时所使用的类别对的索引。由于我们可能不能直接使用原始标签作为索引（尤其是当它们不是从1开始的连续整数时），我们需要一个从1开始的连续整数索引来标识不同的类别。这个索引在训练SVM分类器和投票过程中都是必要的。例如，如果UniqueClasses是[-1, 0, 2]，我们可以为这些类别分配一个从1开始的索引：-1对应索引1，0对应索引2，2对应索引3。然后，ClassPairsIndex就会存储这些索引对应的类别对，例如[1, 2]、[1, 3]和[2, 3]。在预测过程中使用这些字段在预测过程中，SVMModel.UniqueClasses和SVMModel.ClassPairsIndex一起使用来确保正确的类别映射。以下是一个简化的流程：对于测试集中的每个样本，使用所有的一对一SVM分类器进行预测。根据每个分类器的预测结果（正类或负类）进行投票。但是，由于我们使用了从1开始的连续整数索引来标识类别，因此投票结果也是基于这些索引的。找到得票最多的索引（例如，索引2）。使用SVMModel.ClassPairsIndex（如果必要的话）和SVMModel.UniqueClasses将得票最多的索引转换回原始标签。在这个例子中，索引2对应UniqueClasses中的第二个元素，即标签0。返回预测的原始标签作为预测结果。通过这种方式，即使训练数据中的标签不是从1开始的连续整数，我的SVM分类器实现也能够正确地预测类别并返回原始标签。软间隔的SVM（支持向量机）算法软间隔的SVM（支持向量机）算法是为了解决线性不可分或存在噪声数据的问题而引入的。在硬间隔SVM中，所有样本都必须被正确地分类，这在实际应用中往往很难满足。软间隔SVM则允许部分样本被错误分类，通过在目标函数中引入一个惩罚项（或称为损失函数）来控制这种错误分类的程度。软间隔SVM的算法可以通过二次规划问题来描述。二次规划问题是一个优化问题，它涉及最小化或最大化一个二次函数，同时满足一组线性不等式约束。在软间隔SVM中，这个二次函数通常与超平面的参数（如权重向量w和偏置项b）有关，而约束条件则与样本点到超平面的距离和分类正确性有关。下面是软间隔SVM二次规划问题的标准形式：灵感来自卢朓老师《数值方法：原理、算法及应用》课堂 学生大作业~学生作业的原代码是用MATLAB写的，以上为卢朓老师修改之后可以在北太天元上运行的代码分享原视频见下方：
卢朓老师《数值方法：原理、算法及应用》课堂 学生大作业~2018年全国大学生数学建模竞赛A题：特种防护服设计供大家参考学习~
公众号