数据科学与机器学习(第 03 部分):矩阵回归(基础篇)
动态多元回归的输入扩展瓶颈
在 MT5 上做矩阵回归时,最现实的约束不是算法本身,而是模型能否动态吃进更多自变量。前两篇里作者卡了很久的点就在这里:写死的输入维度一多,策略就没法覆盖成百上千条数据。 外汇与贵金属行情里,可用的驱动因子往往几十到上百个(利差、波动率、跨品种价差等),若回归模型只能硬编码 5~10 个输入,实盘里基本没法用。高杠杆下这类多因子建模误差会被放大,属高风险尝试。 破局思路是让程序在运行时动态申请输入列,而不是编译期定死。这样同一套回归内核,既能跑 3 个因子,也能直接接 200 个,不需要改结构重编译。
◍ 先搞懂矩阵长什么样
给没碰过线性代数的人补一句:矩阵就是按行和列排起来的矩形数字(或对象)表,用来装一批同构数据及其关系。 看个最直白的例子——假设说「房间里有一头大象」这种离散状态,也能映射成 2 行 3 列的格子,写作 2x3 矩阵:行数在前、列数在后,读法是「行 x 列」。 现代计算机啃大数据基本都靠这种结构,因为矩阵底层是连续数组存储,CPU 和内存寻址能直接批处理,不用一条条解释。后面接机器学习时,这一层数组表示正是特征张量的起点。
「用矩阵写出简单线性回归」
线性回归的本质是找一条拟合数据的直线 y = β0 + β1x + ε,其中 ε 是误差项,β0 为 y 轴截距,β1 为斜率。把这个方程写成向量形式,就能借矩阵运算一口气算出所有系数,而不必手工逐点描线。 在单变量场景下,设计矩阵 x 的第一列全为 1(对应截距项),第二列是样本的自变量观测值。模型可表达为 y = x·β + ε,其中 β 是包含 β0 和 β1 的二维向量。 我们真正关心的是普通最小二乘估计下的 β 向量,其闭式解为 β = (xᵀx)⁻¹xᵀy。xᵀx 必然是一个 2×2 对称矩阵,因为 xᵀ 的列数等于 x 的行数,转置相乘后主对角线两侧元素相同。 在 MT5 里你可以直接构造这类矩阵验证:随便拉一段 XAUUSD 的收盘序列作为 y,把 [1, 序号] 排成 x,用矩阵库算一遍 (xᵀx)⁻¹xᵀy,就能拿到这段样本里肉眼难辨的斜率。外汇与贵金属杠杆高、跳空频繁,回归系数仅反映历史样本,下一段走势仍可能明显偏离。
给回归矩阵首列灌入常数 1
做矩阵回归的第一步,是把每一行最前面的那一列全部填成 1,一直铺到末行。这个动作在库里的 Init() 函数完成,目的是给后续的最小二乘计算留出截距项的位置,等真正算权重时你会看到它的用处。 从打印出来的设计矩阵能直观确认:前 714 行(索引 0 到 714)每行前 21 个值全是 1.0,直到索引 735 那一行,第 10 个位置起才出现真实 x 数据,例如 4173.8、4179.2、4182.7,而该行前 9 个仍是 1.0。也就是说填充值结束的地方,恰好是 x 值开始的地方。 x 转置(xT)本质上是把矩阵的行和列对调。由于我们在采集数据时本来就是按转置形态存的,所以乘矩阵时可以跳过对整块的转置,但要把 x 值本身取消转置,才能和已经转置好的 x 矩阵相乘。那个未转置的 nx2 形态实际上是 [1 x1 1 x2 … 1 xn] 这样的交错排列。 下面这段 Init() 就是上述逻辑的代码落地,注意 ArrayFill 那一行才是真正灌 1 的操作。
class="type">void CSimpleMatLinearRegression::Init(class="type">class="kw">double &x[],class="type">class="kw">double &y[], class="type">bool debugmode=true) { ArrayResize(Betas,class="num">2); class=class="str">"cmt">//since it is simple linear Regression we only have two variables x and y if (ArraySize(x) != ArraySize(y)) Alert("There is variance in the number of independent variables and dependent variables \n Calculations may fall class="type">short"); m_rowsize = ArraySize(x); ArrayResize(m_xvalues,m_rowsize+m_rowsize); class=class="str">"cmt">//add one row size space for the filled values ArrayFill(m_xvalues,class="num">0,m_rowsize,class="num">1); class=class="str">"cmt">//fill the first row with one(s) here is where the operation is performed ArrayCopy(m_xvalues,x,m_rowsize,class="num">0,WHOLE_ARRAY); class=class="str">"cmt">//add x values to the array starting where the filled values ended ArrayCopy(m_yvalues,y); m_debug=debugmode; }
◍ 把转置矩阵掰回原样再乘
从 csv 读进来的设计矩阵在打印时往往是一行全 1、一行观测值交替铺开,这种布局其实是转置态。原始片段里转置矩阵行 0 全是 1,行 720 末段出现 4297、4321、4402、4416 这类价格观测值,说明维度是观测数 × 2,而不是回归要的 2 × 观测数。 取消转置就是行列互换,逻辑和普通转置完全一致。把列当行、行当列之后,xT 变成 2 × n,x 变成 n × 2,两者相乘得到 2 × 2 的 xTx。矩阵乘法硬性前提是左矩阵列数等于右矩阵行数,否则 MT5 里跑矩阵运算会直接报错。 具体乘出来四项:行1×列1 因两边都是 1,结果就是观测数 n;行1×列2 与行2×列1 都等于 x 的总和;行2×列2 是 x 的平方和。所以 xTx 左上角那个数字,就是你能直接读取的观测样本量,不用再单独数 csv 行数。 下面这段 MQL5 把取消转置落了地。tr_cols 设成 1+1,是因为单一自变量外还要留一列全 1 的截距项空间;MatrixUnTranspose 跑完再 Print,就能看到 1 配 4248、4201、4352、4402…4416 的竖排结构。
class="type">int tr_rows = m_rowsize, tr_cols = class="num">1+class="num">1; class=class="str">"cmt">//since we have one independent variable we add one for the space created by those values of one MatrixUnTranspose(m_xvalues,tr_cols,tr_rows); Print("UnTransposed Matrix"); MatrixPrint(m_xvalues,tr_cols,tr_rows); UnTransposed Matrix [ class="num">1 class="num">4248 class="num">1 class="num">4201 class="num">1 class="num">4352 class="num">1 class="num">4402 ... ... class="num">1 class="num">4405 class="num">1 class="num">4416 ] class=class="str">"cmt">//inside MatrixRegTest.mq5 script class="macro">#include "MatrixRegression.mqh"; class="macro">#include "LinearRegressionLib.mqh"; CSimpleMatLinearRegression matlr; CSimpleLinearRegression lr; class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Script program start function | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void OnStart() { class=class="str">"cmt">//class="type">class="kw">double x[] = {class="num">651,class="num">762,class="num">856,class="num">1063,class="num">1190,class="num">1298,class="num">1421,class="num">1440,class="num">1518}; //stands for sales class=class="str">"cmt">//class="type">class="kw">double y[] = {class="num">23,class="num">26,class="num">30,class="num">34,class="num">43,class="num">48,class="num">52,class="num">57,class="num">58}; //money spent on ads class=class="str">"cmt">//--- class="type">class="kw">double x[], y[]; class="type">class="kw">string file_name = "NASDAQ_DATA.csv", delimiter = ","; lr.GetDataToArray(x,file_name,delimiter,class="num">1); lr.GetDataToArray(y,file_name,delimiter,class="num">2); } CSimpleLinearRegression lr;
「把一维数组硬算成 xTx 矩阵」
前面那段 xT[] 其实就是把 csv 读进来的 x 值原样复制了一份,纯粹是为了把转置这件事说清楚——真正参与运算的是库里的全局数组 m_xvalues[],它早就存好了 x 值。拿 xT 去乘 m_xvalues,这一步发生在 MatrixMultiply() 内部,逻辑上就是 x 的转置乘 x,得到 xTx。 之所以代码看起来绕,是因为作者刻意没用二维数组 Matrix[i][k],MQL5 多维数组在传参和 resize 上有不少坑,所以他用一维数组加索引偏移来模拟行列。比如 mat1_index = k + (i*row2) 就是把第 i 行第 k 列拍平成一维下标,懂了这个,后面逆矩阵才好接。 跑完乘法后用 MatrixPrint() 把 xTx 打出来,注意输出时传的行列数都是 tr_cols,精度参数给的 5 位。你会在终端看到 xTx 左上角第一个元素等于样本观测总数——这正是设计矩阵第一列全填 1 的回报,截距项就靠它估计。 下一步就是求 xTx 的逆,外汇和贵金属行情噪声大,这类线性回归矩阵若接近奇异,逆算出来可能失真,实盘前务必在 MT5 用历史数据验证条件数。
MatrixMultiply(xT,m_xvalues,xTx,tr_cols,tr_rows,tr_rows,tr_cols); Print("xTx"); MatrixPrint(xTx,tr_cols,tr_cols,class="num">5); class=class="str">"cmt">//remember?? the output of the matrix will be the row1 and col2 marked in red class="type">void CSimpleMatLinearRegression::MatrixMultiply(class="type">class="kw">double &A[],class="type">class="kw">double &B[],class="type">class="kw">double &output_arr[],class="type">int row1,class="type">int col1,class="type">int row2,class="type">int col2) { class=class="str">"cmt">//--- class="type">class="kw">double MultPl_Mat[]; class=class="str">"cmt">//where the multiplications will be stored if (col1 != row2) Alert("Matrix Multiplication Error, \n The number of columns in the first matrix is not equal to the number of rows in second matrix"); else { ArrayResize(MultPl_Mat,row1*col2); class="type">int mat1_index, mat2_index; if (col1==class="num">1) class=class="str">"cmt">//Multiplication for 1D Array { for (class="type">int i=class="num">0; i<row1; i++) for(class="type">int k=class="num">0; k<row1; k++) { class="type">int index = k + (i*row1); MultPl_Mat[index] = A[i] * B[k]; } class=class="str">"cmt">//Print("Matrix Multiplication output"); class=class="str">"cmt">//ArrayPrint(MultPl_Mat); } else { class=class="str">"cmt">//if the matrix has more than class="num">2 dimensionals for (class="type">int i=class="num">0; i<row1; i++) for (class="type">int j=class="num">0; j<col2; j++) { class="type">int index = j + (i*col2); MultPl_Mat[index] = class="num">0; for (class="type">int k=class="num">0; k<col1; k++) { mat1_index = k + (i*row2); class=class="str">"cmt">//k + (i*row2) mat2_index = j + (k*col2); class=class="str">"cmt">//j + (k*col2) } } } } }