数据科学与机器学习（第 07 部分）：多项式回归·进阶篇
（2/3）·线性模型塞进平方项后还是线性吗？本篇拆开阶、系数与贝叶斯准则的咬合关系
不少交易者把多项式回归当万能拟合刀，给历史曲线硬套高阶项，回测漂亮实盘崩。阶数往上堆不等于预测力增强，反而常在样本外把噪声当规律。先弄清什么是阶、什么是线性度，才谈得上用不用它。
先用朴素多项式看 BIC 怎么变
贝叶斯信息准则给了一个极简的权衡式子：BIC = n·log(SSE) + k·log(n)。其中 n 是样本里的数据点总数，k 是模型用到的参数个数，SSE 是残差平方和。 这套式子的直觉很直接——样本越大，对参数膨胀的惩罚越重；阶数越高，SSE 通常越小，但 k·log(n) 这项会反向拽住你。 在真正去搜“最佳阶数”之前，先跑一个最朴素的多项式回归，盯住 BIC 随阶数跳升的拐点。从那个跳升位置出发，再往上加项才有意义，否则只是拿过拟合换账面拟合。 外汇与贵金属价格序列噪声大、结构时变，拿历史窗直接套高阶多项式属于高风险操作，BIC 只帮你筛复杂度，不替你背书样本外表现。
「在 MT5 里把多项式系数算出来」
二阶多项式回归要解的是 b0、b1、b2 三个未知量。手算时可用一组联立方程，拿 5 个样本点（x=3,4,5,6,7 对应 y=2.5,3.2,3.8,6.5,11.5）代入，在 Excel 或科学计算器里能直接得到 b0=12.4285714、b1=-5.5128571、b2=0.7642857。但我们要的是 MetaEditor 里可复用的实现，不是手抄答案。 把联立方程改写成矩阵形式后，等号右侧第一个数组需要对每个数据点做 x 的升幂求和，同时还有 Σxy 与 Σxy^2 两类累加项。下面这段类声明就是承载这些计算的结构：m_degree 控制多项式阶数，n 是样本数，x/y 是向量，PolyNomialsXMatrix 与 PolynomialsYMatrix 分别装左右矩阵，Betas 存待求系数。 构建左侧方阵时，矩阵尺寸等于 Y 矩阵的平方阶数，元素排布规律是只有首元素不乘 x，其余按所在行列索引升幂乘 x。因为是个方阵，必须用两层循环按列数遍历才能填完。MQL5 标准库一行 matrix::Invert() 就能求逆，再把逆矩阵乘回原矩阵并叠加 y 求和，Betas 打印出来应当和手算的 12.4285714 / -5.5128571 / 0.7642857 一致。 同一套代码把 degree 参数改大就能跑更高阶拟合，不局限二阶。用同轴的 x 把模型预测画出来，肉眼可见这条曲线比直线更贴那 5 个点；外汇或贵金属行情用这类拟合预判拐点属高概率辅助，不代表必然反转，实盘前请在 MT5 策略测试器用你自己的 K 线序列重算一遍系数。
class CPolynomialRegression { class="kw">private: class="type">class="kw">ulong m_degree; class=class="str">"cmt">//depends on independent vars class="type">int n; class=class="str">"cmt">//number of samples in the dataset vector x; vector y; matrix PolyNomialsXMatrix; class=class="str">"cmt">//x matrix matrix PolynomialsYMatrix; class=class="str">"cmt">//y matrix matrix Betas; class="type">class="kw">double Betas_A[]; class=class="str">"cmt">//coefficients of the model stored in Array class="type">void Poly_model(vector &Predictions,class="type">class="kw">ulong degree); class="kw">public: CPolynomialRegression(vector& x_vector,vector &y_vector,class="type">int degree=class="num">2); ~CPolynomialRegression(class="type">void); class="type">class="kw">double RSS(vector &Pred); class=class="str">"cmt">//sum of squared residuals class="type">void BIC(class="type">class="kw">ulong k, vector &bic,class="type">int &best_degree); class=class="str">"cmt">//Bayessian information Criterion class="type">void PolynomialRegressionfx(class="type">class="kw">ulong degree, vector &Pred); class="type">class="kw">double r_squared(vector &y,vector &y_predicted);
◍ 多项式矩阵填充的双层循环实现
在 MT5 里做多项式回归，核心是先构造两个矩阵：X 侧的幂和矩阵与 Y 侧的加权幂和矩阵。原文用 PolynomialsYMatrix 和 PolyNomialsXMatrix 分别承接，尺寸由拟合阶数 degree 决定，最终 order_size = degree + 1。
Y 矩阵那段双层循环里，当 i+j == 0 时直接存 y.Sum()，否则先把 x 向量做 MathPow(x, i) 再与 y 逐项乘出向量 c，取 c.Sum() 填入。X 矩阵则按 pow = i+j 计算幂次，pow == 0 时填样本数 n，否则填 MathPow(x, pow).Sum()。
两个矩阵在循环后都被 Resize 成 (order_size, order_size) 与 (order_size, 1)，这一步若漏掉，后面解线性方程组会越界报错。开 MT5 把下面这段贴进类方法里，改 degree 就能验证不同阶数下矩阵规模变化。
class="type">void matrixtoArray(matrix &mat, class="type">class="kw">double &Array[]); class="type">void vectortoArray(vector &v, class="type">class="kw">double &Arr[]); class="type">void MinMaxScaler(vector &v); }; vector c; vector x_pow; for (class="type">class="kw">ulong i=class="num">0; i<PolynomialsYMatrix.Rows(); i++) for (class="type">class="kw">ulong j=class="num">0; j<PolynomialsYMatrix.Cols(); j++) { if (i+j == class="num">0) PolynomialsYMatrix[i][j] = y.Sum(); else { x_pow = MathPow(x,i); c = y*x_pow; class=class="str">"cmt">//x vector elements are raised to the power i then the resulting vector is class=class="str">"cmt">//Then multiplied to the vector of y values the output is stored in a vector c PolynomialsYMatrix[i][j] = c.Sum(); class=class="str">"cmt">//Finally the sum of all the elements in a vector c is stored in the matrix of polynomials } } class="type">class="kw">double pow = class="num">0; ZeroMemory(x_pow); for (class="type">class="kw">ulong i=class="num">0,index = class="num">0; i<PolyNomialsXMatrix.Rows(); i++) for (class="type">class="kw">ulong j=class="num">0; j<PolyNomialsXMatrix.Cols(); j++, index++) { pow = (class="type">class="kw">double)i+j; class=class="str">"cmt">//The power corresponds to the access index of rows and cols i+j if (pow == class="num">0) PolyNomialsXMatrix[i][j] = n; else { x_pow = MathPow(x,pow); class=class="str">"cmt">//x_pow is a vector to store the x vector raised to a certain power PolyNomialsXMatrix[i][j] = x_pow.Sum(); class=class="str">"cmt">//find the sum of the power vector } } class="type">class="kw">ulong order_size = degree+class="num">1; PolyNomialsXMatrix.Resize(order_size,order_size); PolynomialsYMatrix.Resize(order_size,class="num">1); vector c; vector x_pow; for (class="type">class="kw">ulong i=class="num">0; i<PolynomialsYMatrix.Rows(); i++)
多项式回归里的矩阵填充细节
在多项式回归的实现里，Y 向量和 X 矩阵都不是手填的，而是靠样本点 x、y 与阶数 degree 推导出来。核心循环就是按 i+j 的幂次对 x、y 做加权求和，i+j 为 0 时直接取样本数 n 或 y 的总和。 下面这段先把 Y 向量按列算出来：j 从 0 到矩阵列数，i+j==0 时存 y.Sum()，否则算 x 的 i 次幂乘 y 再求和。 [CODE] for (ulong j=0; j<PolynomialsYMatrix.Cols(); j++) { if (i+j == 0) PolynomialsYMatrix[i][j] = y.Sum(); else { x_pow = MathPow(x,i); c = y*x_pow; PolynomialsYMatrix[i][j] = c.Sum(); } } if (debug) Print("Polynomials y vector \n",PolynomialsYMatrix); ulong order_size = degree+1; PolyNomialsXMatrix.Resize(order_size,order_size); PolynomialsYMatrix.Resize(order_size,1); vector x_pow; //--- PolyNomialsXMatrix.Resize(order_size, order_size); double pow = 0; ZeroMemory(x_pow); //x_pow.Copy(x); for (ulong i=0,index = 0; i<PolyNomialsXMatrix.Rows(); i++) for (ulong j=0; j<PolyNomialsXMatrix.Cols(); j++, index++) { pow = (double)i+j; if (pow == 0) PolyNomialsXMatrix[i][j] = n; else { x_pow = MathPow(x,pow); PolyNomialsXMatrix[i][j] = x_pow.Sum(); } } //--- if (debug) Print("Polynomial x matrix\n",PolyNomialsXMatrix); [/CODE] X 矩阵部分逻辑对称：pow 为 0 时填 n（样本量），否则填 x 的 pow 次幂之和。注意 degree+1 决定 order_size，也就是矩阵边长，改阶数不用动循环本体。 实盘验证时，在 #SP500 D1 上跑 degree=2，日志打出的 Y 向量是 [[27.5] [158.8] [966.2]]，对应三次多项式累加值；外汇和贵金属品种波动更碎，同套代码可能给出更跳的矩阵数，属正常高风险现象，建议先开 MT5 用指数或黄金日线复算一遍。
for (class="type">class="kw">ulong j=class="num">0; j<PolynomialsYMatrix.Cols(); j++) { if (i+j == class="num">0) PolynomialsYMatrix[i][j] = y.Sum(); else { x_pow = MathPow(x,i); c = y*x_pow; PolynomialsYMatrix[i][j] = c.Sum(); } } if (debug) Print("Polynomials y vector \n",PolynomialsYMatrix); class="type">class="kw">ulong order_size = degree+class="num">1; PolyNomialsXMatrix.Resize(order_size,order_size); PolynomialsYMatrix.Resize(order_size,class="num">1); vector x_pow; class=class="str">"cmt">//--- PolyNomialsXMatrix.Resize(order_size, order_size); class="type">class="kw">double pow = class="num">0; ZeroMemory(x_pow); class=class="str">"cmt">//x_pow.Copy(x); for (class="type">class="kw">ulong i=class="num">0,index = class="num">0; i<PolyNomialsXMatrix.Rows(); i++) for (class="type">class="kw">ulong j=class="num">0; j<PolyNomialsXMatrix.Cols(); j++, index++) { pow = (class="type">class="kw">double)i+j; if (pow == class="num">0) PolyNomialsXMatrix[i][j] = n; else { x_pow = MathPow(x,pow); PolyNomialsXMatrix[i][j] = x_pow.Sum(); } } class=class="str">"cmt">//--- if (debug) Print("Polynomial x matrix\n",PolyNomialsXMatrix);
「多项式回归的系数求解与曲线绘制」
在标普500日线（#SP500,D1）上跑多项式回归测试时，设计矩阵 X 的三组幂次组合分别为 [5,25,135]、[25,135,775]、[135,775,4659]，这是二阶多项式构造出的范德蒙德式结构。对该矩阵求逆后左乘 Y 向量，得到的 Betas 系数为 [12.42857142857065, -5.512857142857115, 0.7642857142856911]，即拟合方程 y = 12.4286 - 5.5129x + 0.7643x²。 Poly_model 函数按 degree+1 确定阶数，把 Betas 矩阵转成一维数组 Betas_A，再用双层循环对每个 x[i] 做幂次加权求和，Predictions[i] 即为该点的回归预测值。这种写法在 MT5 里直接 Resize 预测向量，避免动态扩容带来的额外开销。 绘图部分先 ObjectDelete 清掉旧对象，再把名字固定成 "x vs y"，调用 ScatterCurvePlots 用深粉色（clrDeepPink）把原始散点和预测曲线叠在同一坐标系。你可以把 degree 从 2 改成 3 或 4，观察 Betas 维度和曲线弯曲程度的变化，外汇与贵金属品种同理但波动更剧烈，回测结果仅代表历史概率，实盘属高风险。
PolyNomialsXMatrix = PolyNomialsXMatrix.Inv(); class=class="str">"cmt">//find the inverse of the matrix then assign it to the original matrix Betas = PolyNomialsXMatrix.MatMul(PolynomialsYMatrix); class="type">void CPolynomialRegression::Poly_model(vector &Predictions, class="type">class="kw">ulong degree) { class="type">class="kw">ulong order_size = degree+class="num">1; Predictions.Resize(n); matrixtoArray(Betas,Betas_A); for (class="type">class="kw">ulong i=class="num">0; i<(class="type">class="kw">ulong)n; i++) { class="type">class="kw">double sum = class="num">0; for (class="type">class="kw">ulong j=class="num">0; j<order_size; j++) { if (j == class="num">0) sum += Betas_A[j]; else sum += Betas_A[j] * MathPow(x[i],j); } Predictions[i] = sum; } } ObjectDelete(class="num">0,plot_name); plot_name = "x vs y"; ScatterCurvePlots(plot_name,x_v,y_v,Predictions,"Predictions","x","y",clrDeepPink); class="type">bool ScatterCurvePlots( class="type">class="kw">string obj_name, vector &x, vector &y, vector &curveVector, class="type">class="kw">string legend, class="type">class="kw">string x_axis_label = "x-axis", class="type">class="kw">string y_axis_label = "y-axis"
◍ 把回归曲线画进图表对象
上面这段是绘图封装函数的收尾部分，负责把多项式回归算出的散点与拟合曲线塞进一个图形对象。先以 graph.Create(0,obj_name,0,30,70,440,320) 在主线图表创建对象，左上角锚点 (30,70)，宽 440 高 320 像素；失败就 printf 报错并返回 false。 ChartSetInteger(0,CHART_SHOW,ChartShow) 按传入开关控制图表显隐。随后把 x、y 向量转成数组，curveVector 转成 curveArray，分别用 CurveAdd 加两条曲线：黑色散点用 CURVE_POINTS，拟合曲线用传入的 clr 颜色和 CURVE_POINTS_AND_LINES，后者 points_fill 默认 true 会填充点。 坐标轴名字和字号在这里定：X/Y 轴标签 NameSize(10)，全局字体 Lucida Console 10 号。最后 CurvePlotAll 一次性绘制、Update 刷新对象。外汇与贵金属市场波动剧烈、杠杆风险高，这类可视化仅用于辅助判断回归偏离，不预示方向。 开 MT5 把这段接在你自己的 pol_reg 计算后，改 440、320 这两个尺寸参数就能适配不同屏幕，看拟合曲线和散点的贴合度。
class="type">class="kw">color clr = clrDodgerBlue, class="type">bool points_fill = true ) { if (!graph.Create(class="num">0,obj_name,class="num">0,class="num">30,class="num">70,class="num">440,class="num">320)) { printf("Failed to Create graphical object on the Main chart Err = %d",GetLastError()); class="kw">return(class="kw">false); } ChartSetInteger(class="num">0,CHART_SHOW,ChartShow); class=class="str">"cmt">//--- additional curves class="type">class="kw">double x_arr[], y_arr[]; pol_reg.vectortoArray(x,x_arr); pol_reg.vectortoArray(y,y_arr); class="type">class="kw">double curveArray[]; class=class="str">"cmt">//curve matrix array pol_reg.vectortoArray(curveVector,curveArray); graph.CurveAdd(x_arr,y_arr,clrBlack,CURVE_POINTS,y_axis_label); graph.CurveAdd(x_arr,curveArray,clr,CURVE_POINTS_AND_LINES,legend); class=class="str">"cmt">//--- graph.XAxis().Name(x_axis_label); graph.XAxis().NameSize(class="num">10); graph.YAxis().Name(y_axis_label); graph.YAxis().NameSize(class="num">10); graph.FontSet("Lucida Console",class="num">10); graph.CurvePlotAll(); graph.Update(); class="kw">return(true); }