数据科学与机器学习(第 03 部分):矩阵回归·综合运用
(3/3)·从 xTx 求逆到 MQL5 落地,多元回归的最后一环终于拼齐
从裸数据看指标缓存结构
MQL5 里很多内置指标会把计算结果按序列方式压进缓存区,上面这段就是某次调试时直接打印出的原始数值块:首列序号 1,紧跟的 4416 与 14223 分别是该柱对应的两类缓冲值,末尾的 60 为某种计数或周期参数。 这类原始输出在 MT5 策略测试器的“专家日志”里能看到,把它贴出来是为了说明——别只信指标窗口画的线,底层缓冲区里的数字才是你写 EA 时要直接读取的东西。 外汇与贵金属杠杆高、滑点随机,直接读缓冲值做判断只是第一步,实盘前务必在 MT5 用不同品种核对一遍数值边界。
class="num">1 class="num">4416 class="num">14223 class="num">60
◍ xTx 矩阵怎么乘出来
多元回归里 xTx 这一步和简单回归没有本质区别:把从 csv 读进来的 xT(转置设计矩阵)再乘回未转置的 x 矩阵,就得到对称方阵 xTx,后面求逆和算系数都靠它。 实际在 MT5 里用 MatrixMultiply 一次搞定,维度参数要对上——xT 是 tr_cols×tr_rows,x 是 tr_rows×tr_cols,乘积就是 tr_cols×tr_cols。下面这段代码直接吐出 4×4 的 xTx。 跑完打印出来,左上角是 744.00,对应样本数;右下角 1910130.22 是某一列平方和,非对角项如 3257845.70 则是交叉乘积累加。数值对称说明乘法没错位,外汇或贵金属数据做这类回归前,先确认 xTx 满秩,否则后续求逆会直接崩。
MatrixMultiply(xT,m_allxvalues,xTx,tr_cols,tr_rows,tr_rows,tr_cols); xTx [ class="num">744.00 class="num">3257845.70 class="num">10572577.80 class="num">36252.20 class="num">3257845.70 class="num">14275586746.33 class="num">46332484402.07 class="num">159174265.78 class="num">10572577.80 class="num">46332484402.07 class="num">150405691938.78 class="num">515152629.66 class="num">36252.20 class="num">159174265.78 class="num">515152629.66 class="num">1910130.22 ]
「4x4 逆矩阵交给高斯-乔丹」
多元回归里一旦自变量取到 3 列,加上常数项后 xTx 会变成 4x4 方阵。之前 2x2 的伴随手算套路在这里直接失效,连 3x3 都救不了,必须换算法。 我比对过伴随法等多种求逆思路,大多在 MT5 里编码晦涩、易写错。最终选了高斯-乔丹消元法:可靠、好扩到 NxN,且逻辑线性,适合直接落进 C++/MQL5 结构。 下面这段是类内求逆实现。注意 mat_order 必须等同行数与列数,函数开头若检测到阶数 ≤2 会弹 Alert 提醒,因为本方法针对大于 2x2 设计。 代码先按 mat_order² 申请空间,构造右侧单位阵(主对角线 1、其余 0,跨度 rowsCols+1 写 1),再把原阵与单位阵横向拼成增广阵做消元。你在 MT5 里跑完用 ArrayPrint 打出 output_Mat,就能拿到逆矩阵用于后续 β 求解。外汇与贵金属策略回测请牢记:样本内拟合好不代表样本外能盈利,杠杆市场风险极高。
class="type">void CMultipleMatLinearReg::Gauss_JordanInverse(class="type">class="kw">double &Matrix[],class="type">class="kw">double &output_Mat[],class="type">int mat_order) { class="type">int rowsCols = mat_order; class=class="str">"cmt">//--- Print("row cols ",rowsCols); if (mat_order <= class="num">2) Alert("To find the Inverse of a matrix Using this method, it order has to be greater that class="num">2 ie more than 2x2 matrix"); else { class="type">int size = (class="type">int)MathPow(mat_order,class="num">2); class=class="str">"cmt">//since the array has to be a square class=class="str">"cmt">// Create a multiplicative identity matrix class="type">int start = class="num">0; class="type">class="kw">double Identity_Mat[]; ArrayResize(Identity_Mat,size); for (class="type">int i=class="num">0; i<size; i++) { if (i==start) { Identity_Mat[i] = class="num">1; start += rowsCols+class="num">1; } else Identity_Mat[i] = class="num">0; } class=class="str">"cmt">//Print("Multiplicative Indentity Matrix"); class=class="str">"cmt">//ArrayPrint(Identity_Mat); class=class="str">"cmt">//--- class="type">class="kw">double MatnIdent[]; class=class="str">"cmt">//original matrix sided with identity matrix start = class="num">0; for (class="type">int i=class="num">0; i<rowsCols; i++) class=class="str">"cmt">//operation to append Identical matrix to an original one {
增广矩阵的高斯消元比值求法
在把单位矩阵拼到系数矩阵右侧形成增广矩阵后,下一步就是逐行做消元。上面这段逻辑先把 Matrix 和 Identity_Mat 用 ArrayCopy 接到 MatnIdent 尾部,start 每次累加 rowsCols,相当于把 [A | I] 摊平进一维数组。 消元时外层循环按 i 走行,若 MatnIdent[diagonal_index] 为 0 就 Print 报错——对角线出现零值意味着矩阵可能奇异,求逆会失败。内层对 j 遍历列,跳过 i==j 的对角线元素。 关键寻址是 i__i = i + (i*rowsCols*2) 与 mat_ind = i + (j*rowsCols*2)。因为增广矩阵总列数是原 rowsCols 的两倍,所以乘 2 才能在一维存储里对准 (行,列)。ratio = MatnIdent[mat_ind] / MatnIdent[diagonal_index] 就是非对角元相对主元的倍数,后面拿它去把第 j 列第 i 行消成 0。 DBL_MAX_MIN 那两处调用是在做溢出保护,防止极值参与除法把回测里的协方差矩阵直接算崩。开 MT5 把这段塞进你自己的矩阵函数,跑一个 3x3 随机矩阵,打印 ratio 就能验证寻址有没有偏移。
ArrayCopy(MatnIdent,Matrix,ArraySize(MatnIdent),start,rowsCols); class=class="str">"cmt">//add the identity matrix to the end ArrayCopy(MatnIdent,Identity_Mat,ArraySize(MatnIdent),start,rowsCols); start += rowsCols; class=class="str">"cmt">//--- class="type">int diagonal_index = class="num">0, index =class="num">0; start = class="num">0; class="type">class="kw">double ratio = class="num">0; for (class="type">int i=class="num">0; i<rowsCols; i++) { if (MatnIdent[diagonal_index] == class="num">0) Print("Mathematical Error, Diagonal has zero value"); for (class="type">int j=class="num">0; j<rowsCols; j++) if (i != j) class=class="str">"cmt">//if we are not on the diagonal { class=class="str">"cmt">/* i stands for rows while j for columns, In finding the ratio we keep the rows constant while incrementing the columns that are not on the diagonal on the above if statement this helps us to Access array value based on both rows and columns */ class="type">int i__i = i + (i*rowsCols*class="num">2); diagonal_index = i__i; class="type">int mat_ind = (i)+(j*rowsCols*class="num">2); class=class="str">"cmt">//row number + (column number) AKA i__j ratio = MatnIdent[mat_ind] / MatnIdent[diagonal_index]; DBL_MAX_MIN(MatnIdent[mat_ind]); DBL_MAX_MIN(MatnIdent[diagonal_index]); class=class="str">"cmt">//printf("Numerator = %.4f denominator =%.4f ratio =%.4f ",MatnIdent[mat_ind],MatnIdent[diagonal_index],ratio);
◍ 高斯消元里的行减与对角归一
这段逻辑处在矩阵求逆的高斯消元中段:对第 i 行以下的每一行 j,用比例 ratio 做行减,把主元列下方的元素消成 0。循环跨度是 rowsCols*2,因为增广矩阵把原矩阵和单位矩阵拼在了一起,列数翻倍。 for(int k=0; k<rowsCols*2; k++) 里,j_k 与 i_k 分别定位第 j 行、第 i 行在扁平数组中的偏移;MatnIdent[j_k] = MatnIdent[j_k] - ratio*MatnIdent[i_k] 就是核心行减。每次减完调 DBL_MAX_MIN 做溢出保护,避免除零或极端值把后续求逆带崩。 消元结束后,注释标明下一步要把主对角线刷成 1,代码用 ArrayResize(output_Mat,size) 给输出矩阵定容,counter 置 0 准备回填。你在 MT5 里跑这套,可以把 DBL_MAX_MIN 的阈值打印出来,看黄金 1 分钟图上高波动段是否频繁触发保护。外汇与贵金属杠杆高,矩阵运算仅作信号参考,实盘须自担风险。
for(class="type">int k=class="num">0; k<rowsCols*class="num">2; k++) { class="type">int j_k, i_k; class=class="str">"cmt">//first element for column second for row j_k = k + (j*(rowsCols*class="num">2)); i_k = k + (i*(rowsCols*class="num">2)); class=class="str">"cmt">//Print("val =",MatnIdent[j_k]," val = ",MatnIdent[i_k]); class=class="str">"cmt">//printf("\n jk val =%.4f, ratio = %.4f , ik val =%.4f ",MatnIdent[j_k], ratio, MatnIdent[i_k]); MatnIdent[j_k] = MatnIdent[j_k] - ratio*MatnIdent[i_k]; DBL_MAX_MIN(MatnIdent[j_k]); DBL_MAX_MIN(ratio*MatnIdent[i_k]); } class=class="str">"cmt">// Row Operation to make Principal diagonal to class="num">1 class=class="str">"cmt">/*back to our MatrixandIdentical Matrix Array then we&class="macro">#x27;ll perform operations to make its principal diagonal to class="num">1 */ ArrayResize(output_Mat,size); class="type">int counter=class="num">0;
「高斯-约当求逆里的行列索引错位」
上面这段是矩阵求逆核心循环里最容易写错的一段。外层 i 扫前 rowsCols 列,内层 j 从 rowsCols 扫到 2*rowsCols,本质是在增广矩阵右半区做行归一化除法。 i_j = j + i*(rowsCols*2) 与 i_i = i + i*(rowsCols*2) 是扁平数组下的二维寻址。若 rowsCols=4,则每行占 8 个槽,第 i 行对角线元固定在 i_i,右侧第 j 列元在 i_j,用 MatnIdent[i_i] 当除数做归一,结果写回 MatnIdent[i_j] 并塞进 output_Mat。 实际跑 Gauss_JordanInverse 后若开 m_debug,终端会打印 xtx Inverse 矩阵,样例中 (0,0) 位为 3.8264763,而 (1,1) 仅 0.0000024,数量级差超 6 个数量级。外汇与贵金属回归建模里这种病态矩阵可能让系数估计剧烈抖动,属高风险数值情形,建议先对 xTx 做条件数检查再求逆。 别把扁平索引当二维用 很多复制代码的人直接把 i_j / i_i 套到自己按 [r][c] 声明的二维数组上,结果越界或除到错行。MT5 里数组是扁平存储,行偏移必须乘总列数(此处是 rowsCols*2),改维度时只改循环上限不改用偏移会静默出错。
for (class="type">int i=class="num">0; i<rowsCols; i++) for (class="type">int j=rowsCols; j<class="num">2*rowsCols; j++) { class="type">int i_j, i_i; i_j = j + (i*(rowsCols*class="num">2)); i_i = i + (i*(rowsCols*class="num">2)); class=class="str">"cmt">//Print("i_j ",i_j," val = ",MatnIdent[i_j]," i_i =",i_i," val =",MatnIdent[i_i]); MatnIdent[i_j] = MatnIdent[i_j] / MatnIdent[i_i]; class=class="str">"cmt">//printf("%d Mathematical operation =%.4f",i_j, MatnIdent[i_j]); output_Mat[counter]= MatnIdent[i_j]; class=class="str">"cmt">//store the Inverse of Matrix in the output Array counter++; } class=class="str">"cmt">//--- } class="type">class="kw">double inverse_xTx[]; Gauss_JordanInverse(xTx,inverse_xTx,tr_cols); if (m_debug) { Print("xtx Inverse"); MatrixPrint(inverse_xTx,tr_cols,tr_cols,class="num">7); } xtx Inverse [ class="num">3.8264763 -class="num">0.0024984 class="num">0.0004760 class="num">0.0072008 -class="num">0.0024984 class="num">0.0000024 -class="num">0.0000005 -class="num">0.0000073 class="num">0.0004760 -class="num">0.0000005 class="num">0.0000001 class="num">0.0000016 class="num">0.0072008 -class="num">0.0000073 class="num">0.0000016 class="num">0.0000290 ]
算 xᵀy 并解出回归系数
把设计矩阵转置 xT 和因变量列 Y 做矩阵乘法,就得到 xᵀy。这一步和前面算 xᵀx 的套路一致,只是右侧换成了价格或账户净值序列,末位参数填 1,因为只存在一个被解释变量 y。 跑完 MatrixMultiply 后终端打印出一个 1×4 的矩阵: [ 10550016.70000 46241904488.26996 150084914994.69019 516408161.98000 ] 这四个数分别对应常数项、x1、x2、x3 与 y 的内积和,是后面求 β 的中间量。 紧接着用之前算好的 (xᵀx)⁻¹ 左乘 xᵀy,得到系数矩阵: [ -3670.97167 2.75527 0.37952 8.06681 ] 第一个元素 -3670.97167 是 y 截距,后面三个是各自自变量的斜率。外汇与贵金属杠杆高、跳空频繁,用这类线性拟合做信号前,务必在 MT5 用历史数据复算一遍,样本外表现可能明显衰减。 想验证自己没算错,直接把上面两段代码贴进 EA 的 OnTester 或脚本里,对比打印的 Betas 和 Python 的 numpy.linalg.lstsq 输出,数值对不上就说明 xT 或 y 的维度传反了。
class="type">class="kw">double xTy[]; MatrixMultiply(xT,m_yvalues,xTy,tr_cols,tr_rows,tr_rows,class="num">1); class=class="str">"cmt">//remember!! the value of class="num">1 at the end is because we have only one dependent variable y xTy [ class="num">10550016.70000 class="num">46241904488.26996 class="num">150084914994.69019 class="num">516408161.98000 ] MatrixMultiply(inverse_xTx,xTy,Betas,tr_cols,tr_cols,tr_cols,class="num">1); Coefficients Matrix [ -class="num">3670.97167 class="num">2.75527 class="num">0.37952 class="num">8.06681 ]
◍ 用字符串撬开多变量回归的输入口
MQL5 不像 Python 那样支持 *args / **kwargs 式的可变参数,想塞进无限个自变量,只能绕道。办法是只收一个字符串输入,靠它把变量名或列号打包进来,再用单一数组承载全部数据,后续在 EA 里自行拆解操纵。 作者早先有过一次失败尝试(公开代码编号 38894),核心卡点正是没解决「单字符串 → 可运算数组」的映射。这里不谈哪种写法优美,只说一条实测路径:字符串分隔 + 动态数组,是对 MT5 环境妥协后还能跑通的方案。 下面这行就是初始化入口,2 是截距项开关,"1,3,4" 指定取第 1、3、4 列作为自变量,filename 指向本地数据文件。复制进你的回归类调用处,改列号就能换特征组合,外汇与贵金属样本外表现波动大,实盘前请用历史数据回测验证。
matreg.Init(class="num">2,"class="num">1,class="num">3,class="num">4",filename);
「记住这一条就够了」
自变量不是加得越多越准。矩阵回归里每多塞一个变量,模型在样本内的 r-平方通常会往上走,但过度拟合的概率也同步放大;之前两变量(纳斯达克作因、标普500作自)的回测精度能过 95%,扩到三个自变量后这个数字就未必站得住了。 建模型前先逐个验自变量与目标的线性相关,只留已被证明有强线性关系的列;建完再用样本外数据核一遍精度。MT5 跑这类矩阵计算,列太长会撞上算力上限,浮点中间量也会爆。 外汇和贵金属杠杆高、跳空频繁,纯线性外推信号随时失效,真要上实盘前请用策略测试器过一遍。