阿里云服务器添加网站,域名推荐,做网站的哪个好,阿里云服务器价格表1. 灰色预测#xff1a;当数据少得可怜时#xff0c;我们还能做什么#xff1f; 你好#xff0c;我是老张#xff0c;一个在数据分析和算法领域摸爬滚打了十多年的工程师。今天想跟你聊聊一个特别“接地气”的预测方法——灰色预测。你可能听说过各种高大上的预测模型clc year [1995:1:2004]; % 年份列向量从1995到2004 x0 [174,179,183,189,207,234,220.5,256,270,285]; % 原始废水排放数据 n length(x0); % 获取数据长度 n10 % 1. 累加生成得到新序列 x1 x1 cumsum(x0); % cumsum是累加函数比用for循环更简洁高效 % 2. 计算光滑比 rho rho zeros(n,1); % 初始化光滑比数组 for i 2:n % 从第二期开始计算因为第一期没有“之前的总高度” rho(i) x0(i) / x1(i-1); end % 3. 可视化画出光滑比变化曲线 figure(1) plot(year(2:n), rho(2:n), o-, LineWidth, 1.5); % 画折线图 hold on; plot([year(2), year(n)], [0.5, 0.5], r--, LineWidth, 1.5); % 画0.5的临界线 grid on; % 显示网格 xlabel(年份); ylabel(光滑比 \rho(k)); title(原始数据光滑比检验); legend(光滑比, 临界线 (0.5)); set(gca, xtick, year(2:1:end)); % 设置x轴刻度为每一年 % 4. 计算两个关键指标 % 指标1全部光滑比除第一期小于0.5的占比 prop1 sum(rho(2:end) 0.5) / (n-1); % 指标2剔除前两期后光滑比小于0.5的占比 prop2 sum(rho(3:end) 0.5) / (n-2); fprintf(指标1全部光滑比0.5的比例: %.2f%%\n, prop1*100); fprintf(指标2后%d期光滑比0.5的比例: %.2f%%\n, n-2, prop2*100); if prop1 0.6 prop2 0.9 disp(恭喜数据通过级比检验适合使用GM(1,1)模型); else disp(警告数据未通过级比检验直接使用GM(1,1)模型风险较大。); disp(可能需要对原始数据进行平移、对数处理等预处理或考虑其他模型。); end运行这段代码你会看到图形和命令行输出。对于我们的示例数据输出结果应该是两个指标都远超标准线稳稳通过检验。这个步骤至关重要它决定了我们后续所有工作的基础是否牢固。我见过不少朋友跳过检验直接建模结果预测曲线离奇古怪问题往往就出在这第一步。3. GM(1,1)模型核心从累加到微分方程的求解之旅通过了数据体检我们就可以正式请出GM(1,1)模型了。别被“灰色模型”、“一阶微分方程”这些词吓到我会用“搭积木”的方式把整个建模过程拆解给你看。整个过程就像侦探破案先制造线索累加再找到线索背后的规律拟合微分方程最后还原真相预测。3.1 第一步累加生成——让数据“显形”原始数据x0 [174, 179, 183, ...]看起来波动不大但也没什么明显的函数规律。GM(1,1)的第一招叫一次累加生成记作1-AGO。操作非常简单新序列x1的第m项等于原始序列x0前m项的和。x1(1) x0(1) 174 x1(2) x0(1)x0(2) 174179 353 x1(3) x0(1)x0(2)x0(3) 174179183 536 ...你可以用Matlab的cumsum函数一键完成。累加之后x1序列会呈现出明显的单调递增趋势看起来非常接近一条指数曲线或者直线。这一步的目的就是弱化原始数据的随机性强化其内在的规律性。好比把一堆散乱的珠子用线串起来它的整体走向就清晰了。3.2 第二步构建灰微分方程——找到规律的“数学描述”现在我们面对的是看起来有指数增长趋势的x1序列。在数学上什么函数能描述这种趋势指数函数。而指数函数的导数与其自身成正比这自然引出了一阶常微分方程dx1/dt a * x1 u这里a被称为发展系数u被称为灰色作用量。a反映了x1的增长速度u可以理解为系统内的“驱动力”。我们的目标就是根据已知的x1数据把这个方程里的a和u这两个参数给求出来。但我们的数据是离散的没有连续的导数dx1/dt。怎么办用差分来近似微分。在离散时间点t导数近似等于相邻数据的差值dx1/dt ≈ x1(t) - x1(t-1) x0(t)。同时为了减少随机波动的影响我们用x1的相邻两项均值z1(t) 0.5*x1(t) 0.5*x1(t-1)来代替方程中的x1(t)。于是高大上的微分方程就被我们“离散化”成了一个非常亲切的一元线性方程x0(t) -a * z1(t) u看是不是变成了y kx b的形式这里y就是x0(t)x就是z1(t)k是-ab是u。3.3 第三步参数求解——最小二乘法的妙用现在问题简化成了我们有一堆(z1(t), x0(t))的数据点要找一条直线y kx b最好地拟合它们。这就是经典的线性回归问题而最小二乘法正是求解它的利器。我们可以直接套用最小二乘法的公式来求解k和b进而得到a-k, ubn length(x0); y x0(2:end); % 因变量从第二期开始 x z1; % 自变量z1序列 % 最小二乘法公式求解 k ((n-1)*sum(x.*y) - sum(x)*sum(y)) / ((n-1)*sum(x.*x) - sum(x)^2); b (sum(x.*x)*sum(y) - sum(x)*sum(x.*y)) / ((n-1)*sum(x.*x) - sum(x)^2); a -k; u b;当然在Matlab里你也可以直接用polyfit函数或者regress函数需要统计工具箱来求解但理解这个公式能让你更清楚背后的原理。我个人的习惯是在核心算法上自己实现一遍这样调试和排查问题时心里更有底。3.4 第四步预测公式推导——从方程解出未来求出a和u后代回最初的微分方程dx1/dt a*x1 u并求解这个微分方程这是一个标准的一阶线性常微分方程就能得到累加序列x1的时间响应函数也就是它的拟合与预测公式x1_hat(m1) [x0(1) - u/a] * exp(-a*m) u/a其中m 0, 1, 2, ...。当m取0到n-1时算出来的是对历史累加值x1的拟合值当m n时算出来的就是对未来累加值的预测值。但我们最终需要的是原始序列x0的预测值。这好办因为x0是x1的“差分”累加的逆运算。所以对上述公式的结果再做一次累减生成1-IAGOx0_hat(m1) x1_hat(m1) - x1_hat(m) (1 - exp(a)) * [x0(1) - u/a] * exp(-a*m)这个x0_hat就是我们千辛万苦最终得到的——原始数据的拟合值m对应历史期和预测值m对应未来期。整个推导过程从累加到微分方程再到最小二乘和最终解逻辑链条非常完整体现了灰色预测“用已知信息生成规律用规律推测未知”的核心思想。4. 代码封装与实战一个稳健的GM(1,1)预测函数理解了原理我们把上面的步骤打包成一个健壮的、可重复使用的Matlab函数。这个函数我称之为gm11它不仅能预测还会输出拟合效果的评价指标方便我们判断模型好坏。function [result, x0_hat, relative_residuals, eta] gm11(x0, predict_num) % GM(1,1)灰色预测模型 % 输入 % x0: 原始数据序列列向量 % predict_num: 需要向后预测的期数 % 输出 % result: 预测值列向量 % x0_hat: 对原始数据的拟合值列向量 % relative_residuals: 相对残差用于模型评价 % eta: 级比偏差用于模型评价 n length(x0); % 1. 累加生成(1-AGO) x1 cumsum(x0); % 2. 计算紧邻均值生成序列 z1 z1 0.5 * x1(2:end) 0.5 * x1(1:end-1); % 3. 构造数据矩阵利用最小二乘法求解参数 a, u Y x0(2:end); B [-z1, ones(size(z1))]; % 构造系数矩阵 B [-z1, 1] % 最小二乘法求解 (B*B) \ (B*Y) parameters (B * B) \ (B * Y); a parameters(1); % 发展系数 u parameters(2); % 灰色作用量 % 4. 计算拟合值 x0_hat x0_hat zeros(n, 1); x0_hat(1) x0(1); % 第一期的拟合值就是原始值 for m 1:n-1 x0_hat(m1) (1 - exp(a)) * (x0(1) - u/a) * exp(-a * m); end % 5. 计算预测值 result result zeros(predict_num, 1); for i 1:predict_num result(i) (1 - exp(a)) * (x0(1) - u/a) * exp(-a * (n i - 1)); end % 6. 模型评价计算相对残差 absolute_residuals x0(2:end) - x0_hat(2:end); relative_residuals abs(absolute_residuals) ./ x0(2:end); % 7. 模型评价计算级比偏差 % 首先计算原始序列的级比 class_ratio x0(2:end) ./ x0(1:end-1); % 然后计算级比偏差 eta abs(1 - (1 - 0.5*a) / (1 0.5*a) * (1 ./ class_ratio)); end这个函数是我经过多个项目打磨后的版本比网上很多简单的实现要更完整。它有几个亮点一是使用矩阵运算(B*B)\(B*Y)来求解最小二乘比直接套公式更简洁且数值稳定性更好二是同时输出了拟合值和两个重要的模型评价指标三是结构清晰每一步都有对应方便你调试和学习。如何使用它假设我们想用全部10年数据预测未来3年的废水排放量并查看拟合效果% 假设x0和year已经定义好 predict_years 3; % 预测未来3年 [future_values, fitted_values, rel_res, eta] gm11(x0, predict_years); % 打印预测结果 fprintf(未来%d年的预测值为\n, predict_years); for i 1:predict_years fprintf(第%d年: %.2f\n, i, future_values(i)); end % 计算平均评价指标 avg_rel_res mean(rel_res); avg_eta mean(eta); fprintf(\n模型评价\n); fprintf(平均相对残差: %.4f (%.2f%%)\n, avg_rel_res, avg_rel_res*100); fprintf(平均级比偏差: %.4f\n, avg_eta);运行后你不仅能得到预测值还能立刻看到模型对历史数据的拟合精度这对判断预测结果是否可信至关重要。5. 模型好不好让残差与级比偏差告诉你模型建好了预测值也出来了但我们不能“王婆卖瓜自卖自夸”。必须有一套客观的评价标准来检验这个GM(1,1)模型到底靠不靠谱。这里主要看两个指标残差检验和级比偏差检验。它们就像模型的“体检报告”一份给拟合效果一份给模型结构。5.1 残差检验拟合值与真实值差多少残差说白了就是“预测值”和“真实值”的差距。我们主要看相对残差因为它消除了数据本身量纲的影响。计算公式是相对残差 |真实值 - 拟合值| / 真实值我们对所有历史期通常从第二期开始算计算这个值然后求平均得到平均相对残差。怎么判断这是经验准则如果平均相对残差 20%认为模型拟合达到一般要求可以接受。如果平均相对残差 10%认为模型拟合效果非常好预测结果可信度高。如果超过20%甚至更大那就要小心了模型可能并不适合你的数据或者数据需要预处理。在我们的废水排放案例中计算出的平均相对残差约为2.6%远低于10%这说明模型对历史数据的还原能力非常强。5.2 级比偏差检验模型结构与数据本身匹配吗残差检验看的是“点对点”的误差而级比偏差检验则更深入地考察模型揭示的规律是否与数据内在的规律一致。它比较的是由模型参数a推导出的理论级比和原始数据实际计算出的级比两者之间的偏差。计算步骤稍多先计算原始数据x0的级比σ(k) x0(k) / x0(k-1)。利用求出的发展系数a计算理论级比(1 - 0.5a) / (1 0.5a)。计算级比偏差η(k) |1 - (理论级比) / σ(k)|。最后计算平均级比偏差η_avg。判断标准与残差检验类似η_avg 0.2模型可用。η_avg 0.1模型精度很高。我们的案例中平均级比偏差约为0.047同样远小于0.1。两个检验同时通过且指标优秀这给了我们使用该模型进行预测的充分信心。在实际项目中我强烈建议你把这两个检验的结果可视化出来像下面这样figure(Position, [100, 100, 800, 600]) % 子图1绘制相对残差 subplot(2,1,1) plot(year(2:end), rel_res*100, s-, LineWidth, 1.5, MarkerSize, 8) % 转换为百分比 grid on ylabel(相对残差 (%)) title(模型残差检验) % 在图上添加平均线 hold on plot([year(2), year(end)], [avg_rel_res*100, avg_rel_res*100], r--, LineWidth, 1.2) legend(各期相对残差, sprintf(平均线 (%.2f%%), avg_rel_res*100), Location, best) % 子图2绘制级比偏差 subplot(2,1,2) plot(year(2:end), eta, o-, LineWidth, 1.5, MarkerSize, 8) grid on xlabel(年份) ylabel(级比偏差 \eta) title(模型级比偏差检验) hold on plot([year(2), year(end)], [avg_eta, avg_eta], r--, LineWidth, 1.2) legend(各级比偏差, sprintf(平均线 (%.4f), avg_eta), Location, best)这样一张图就能清晰、直观地展示模型在所有历史点上的表现有没有异常值一目了然。如果发现某一期的残差或级比偏差突然变得很大就需要回头检查那期的原始数据是否有问题或者思考是否发生了影响系统的特殊事件。6. 完整预测流程与结果可视化从数据到洞察现在让我们把前面所有步骤串起来形成一个从数据检验、模型构建、预测到结果展示的完整工作流。这个过程也是我日常分析的标准流程。首先我们进行级比检验确认数据可用。接着为了更严谨我习惯将最后2-3期数据留作“试验组”不参与模型训练用来事后检验模型的预测能力。用前7-8期数据训练模型去预测最后2-3期然后与真实值对比这能很好地模拟“预测未来”的场景。% 假设 year 和 x0 已定义 % 1. 划分训练集和试验集 test_num 3; % 取最后3年作为试验集 train_x0 x0(1:end-test_num); test_x0 x0(end-test_num1:end); train_year year(1:end-test_num); test_year year(end-test_num1:end); % 2. 使用训练集训练模型并预测试验集年份 [pred_test, fit_train, ~, ~] gm11(train_x0, test_num); % 3. 可视化对比试验集预测值与真实值 figure(Position, [100, 100, 900, 400]) subplot(1,2,1) plot(test_year, test_x0, bd-, LineWidth, 2, MarkerSize, 10, MarkerFaceColor, b) hold on plot(test_year, pred_test, rs--, LineWidth, 2, MarkerSize, 10, MarkerFaceColor, r) grid on xlabel(年份) ylabel(废水排放量) title(试验集预测效果对比) legend(真实值, GM(1,1)预测值, Location, best) set(gca, xtick, test_year) % 计算试验集上的预测误差 test_error abs(test_x0 - pred_test) ./ test_x0; fprintf(试验集预测相对误差\n); for i 1:test_num fprintf( %d年: 真实值%.1f, 预测值%.1f, 误差%.2f%%\n, ... test_year(i), test_x0(i), pred_test(i), test_error(i)*100); end % 4. 使用全部数据重新训练模型进行真正未来预测 future_num 5; % 预测未来5年 [future_pred, all_fit, rel_res_all, eta_all] gm11(x0, future_num); all_year_pred [year; (year(end)1:year(end)future_num)]; all_data_pred [all_fit; future_pred]; % 5. 绘制完整的拟合与预测图 subplot(1,2,2) plot(year, x0, ko-, LineWidth, 2, MarkerSize, 8, MarkerFaceColor, k) hold on plot(year, all_fit, b^--, LineWidth, 1.5, MarkerSize, 8) plot(year(end)1:year(end)future_num, future_pred, r*--, LineWidth, 1.5, MarkerSize, 12) grid on xlabel(年份) ylabel(废水排放量) title(GM(1,1)模型完整拟合与预测) legend(原始数据, 拟合数据, 预测数据, Location, best) % 添加预测区间标识 plot([year(end), year(end)1], [x0(end), future_pred(1)], r--) set(gca, xtick, [year(1):2:year(end)future_num])运行这段代码你会得到两张非常有说服力的图。左边是模型在“未知”试验集上的表现直观看到预测值与真实值的接近程度。右边则是完整的趋势图包含了历史拟合和未来预测。从我们的案例图可以看到预测曲线平滑地延续了历史增长趋势。最后别忘了输出关键结论fprintf(\n GM(1,1)模型预测报告 \n); fprintf(发展系数 a %.6f\n, a_from_full_model); % 需要从完整模型获取a fprintf(灰色作用量 u %.6f\n, u_from_full_model); fprintf(平均相对残差 %.4f%% 10%%模型拟合优度优秀。\n, mean(rel_res_all)*100); fprintf(平均级比偏差 %.4f 0.1模型结构合理性优秀。\n, mean(eta_all)); fprintf(未来%d年预测结果\n, future_num); for i 1:future_num fprintf( %d年: %.2f\n, year(end)i, future_pred(i)); end fprintf(\n);这样一份从数据预处理、模型检验、参数求解、效果评估到最终预测的完整分析报告就生成了。整个过程在Matlab中实现了自动化你只需要准备好符合格式的原始数据就能快速得到预测结果和可信度评估。在实际应用中比如预测产品销量、能源消耗、用户增长等场景这套流程非常高效实用。