概率论与数理统计示例(第一部分):基础与初级理论·进阶篇
统计与概率到底差在哪
概率论和数理统计常被混为一谈,但边界其实很清楚:概率论默认模型完全已知,往后是推结论;数理统计则是模型只认得一半,靠实验数据去补剩下的认知。前面章节里聊的那些问题,严格说都还停在概率论层面。 还有个更现代的说法,把数理统计塞进决策论里——核心变成构造决策规则,目标是让误差的平均代价最小。机器学习和它在这里有交集也有分叉:统计里模型类型基本框死了(比如未知参数给个精度范围),机器学习连模型长什么样都经常不确定。 传统口径下,数理统计真正要啃的硬骨头是样本问题。样本不是随便抓几个数,而是用有限观测去反推那个不完整的模型,这一步才是后面所有回测和分布检验的地基。
「用极大似然给离散行情估参数」
把资产价格按固定百分比步长离散成 1 和 -1 的阶梯序列,就能直接套伯努利方案做数理统计。原文取 EURUSD 从 2020 年 5 月 20 日到 6 月 20 日、增量 0.5% 的分钟开盘价,得到 27 步足迹:1, -1, 1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, 1, 1, 1, -1, -1, 1, 1, 1, -1, -1, -1, 1, -1, -1。这一步先在 MT5 里跑通离散脚本,后续所有估算才有样本。 参数点估值最实用的是最大似然法(MLE):找使似然函数最大的模型参数。对简单伯努利序列,上移概率 p 的估值就是频率 pe = n_up / n_total。上面那段 EURUSD 离散序列里上移 16 次、下移 11 次,算得 pe = 0.59(四舍五入),意味着该月样本里阶梯向上的概率略高于随机游走的 0.5。 把序列切两段分别估 p1、p2,可得 n1e=21、p1e=0.71、p2e=0.17,模型在尾部“捕捉”到方向切换;再换最简两状态马尔可夫链,用前一步状态决定当前上移概率,估得 p1e=0.56、p2e=0.60,两值接近且都 >0.5,说明相邻走势几乎无相关、整体偏上行。复杂模型没带来额外信息,这种品种周期里用简单伯努利就够了。 检验假设时,空猜想 p0=0.5(公平硬币)在显著水平下未被拒绝:k_up=16 落在 kl=9 与 kr=18 之间。但把序列拆两截做费舍尔精确检验,第二段上移次数 k2=2 ≤ kl=2,空猜想 p1=p2 被拒绝,印证前半段与后半段确实可能来自不同偏置的子过程。外汇与贵金属离散化统计仅描述历史样本特征,实盘应用属高风险,参数估值不等于未来概率。
from scipy.stats class="kw">import hypergeom n = class="num">1000 k = class="num">20 x = class="num">1 lhx = class="num">0.0 be = class="num">0 for b in range(x, n - k + x): lh = hypergeom.pmf(x, n, b, k) if lh > lhx: be = b lhx = lh print("be =",be)
◍ 把连续价格压成离散涨跌步
做价格行为统计前,先把 MT5 的连续报价压成只含 +1 / -1 的离散序列,否则后面跑超几何分布检验会没法对齐样本空间。上面 Python 片段用 scipy 的 hypergeom.pmf 做反向试探:给定 b=50、k=55、命中 x=3,从 ne=b+k-x-1 起循环累加总体数,直到概率密度比上一步小才停,打印出的 ne 就是该参数下的最小总体边界,实测落在 107 附近。 MQL5 这边用 setmv() 干实事:以 dpr 百分比步长(默认 0.5%)把指定时间段 M1 开盘价取对数后分层,每跨一层就往 mv[] 里补相应数量的 1 或 -1。代码里 ArrayResize 预分配 1000 冗余位,CopyOpen 拉不到 2 根以上 K 线会直接 Print 报错退出,避免空数组炸脚本。 主脚本 OnStart 把 2020.5.20–6.20 的 XAUUSD 类标的离散化后,先 Print 出整条 1/-1 串供肉眼看噪点,再用 CGraphic 关掉原图把累加路径画成 750×350 的窗。外汇和贵金属杠杆高、滑点跳层频繁,0.5% 步长对黄金可能偏粗,建议先拿历史数据跑一遍看 mv 总数是否过千再调 dprcnt。
class=class="str">"cmt">// 构造走势的离散 mv[] 数组 class=class="str">"cmt">// 以指定的百分比步长,在指定的时间区间 class="type">void setmv(class="type">int& mv[], class="type">class="kw">datetime t1, class="type">class="kw">datetime t2, class="type">class="kw">double dpr) { class="type">int ND = class="num">1000; ArrayResize(mv, class="num">0, ND); class=class="str">"cmt">// 获取价格历史记录 class="type">class="kw">double price[]; class="type">int nprice = CopyOpen(Symbol(), PERIOD_M1, t1, t2, price); if(nprice < class="num">2) { Print("not enough price history"); class="kw">return; } class=class="str">"cmt">// 构造 mv[] class="type">int lvl = class="num">0, dlvl, nmv = class="num">0, dmv; class="type">class="kw">double lp0 = log(price[class="num">0]), lstep = log(class="num">1 + class="num">0.01 * dpr); for(class="type">int i = class="num">1; i < nprice; ++i) { dlvl = (class="type">int)((log(price[i]) - lp0) / lstep - lvl); if(dlvl == class="num">0) class="kw">continue; lvl += dlvl; dmv = class="num">1; if(dlvl < class="num">0) { dmv = -class="num">1; dlvl = -dlvl; } ArrayResize(mv, nmv + dlvl, ND); for(class="type">int j = class="num">0; j < dlvl; ++j) mv[nmv + j] = dmv; nmv += dlvl; } } class="macro">#include <Discr.mqh> class="macro">#include <Graphics\Graphic.mqh> class="macro">#class="kw">property script_show_inputs class=class="str">"cmt">//+------------------------------------------------------------------+ class="kw">input class="type">class="kw">datetime tstart = D&class="macro">#x27;class="num">2020.05.class="num">20 class="num">00:class="num">00&class="macro">#x27;; class="kw">input class="type">class="kw">datetime tstop = D&class="macro">#x27;class="num">2020.06.class="num">20 class="num">00:class="num">00&class="macro">#x27;; class="kw">input class="type">class="kw">double dprcnt = class="num">0.5; class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void OnStart() { class="type">int mv[], nmv; setmv(mv, tstart, tstop, dprcnt); nmv = ArraySize(mv); if(nmv < class="num">1) { Print("not enough moves"); class="kw">return; } class=class="str">"cmt">// 显示 mv[] 为 class="num">1 和 -class="num">1 的序列 class="type">class="kw">string res = (class="type">class="kw">string)mv[class="num">0]; for(class="type">int i = class="num">1; i < nmv; ++i) res += ", " + (class="type">class="kw">string)mv[i]; Print(res); class=class="str">"cmt">// 显示 mv[] 累加汇总图表 ChartSetInteger(class="num">0, CHART_SHOW, false); CGraphic graphic; graphic.Create(class="num">0, "G", class="num">0, class="num">0, class="num">0, class="num">750, class="num">350);
从离散步幅推算上行概率
把价格按固定百分比离散化后,每一段涨跌就变成整数数组 mv[],正数为向上步幅、负数为向下步幅。直接数出其中大于 0 的元素个数,再除以总步数,就能得到一个朴素的上行概率估值 pe,作为后续统计的基准。 下面这段函数只做一件事:遍历 mv 统计上行步数 nup,返回 nup/nmv。若传入数组长度不足 1,直接返回 0.0 避免除零。 double calcpe(int& mv[]) { int nmv = ArraySize(mv); if(nmv < 1) return 0.0; int nup = 0; for(int i = 0; i < nmv; ++i) if(mv[i] > 0) ++nup; return ((double)nup) / nmv; } 实际脚本里,输入参数把研究区间锁在 2020.05.20 到 2020.06.20、步长百分比 dprcnt=0.5,跑完 setmv 后若 ArraySize(mv)<1 会打印 not enough moves 并退出。外汇与贵金属波动剧烈,该 pe 仅反映历史片段中的方向偏好,样本外可能明显偏移,切勿当作方向确定性依据。 进一步可拆出 n1e(连续上行段数)与 p1e、p2e(不同条件概率),当 ArraySize(mv)<2 时脚本同样拒绝计算。这种分层估值比单 pe 更耐看,但本质仍是离散重构,参数 dprcnt 从 0.5 调到 0.3 会让步数翻倍、pe 序列更抖,建议开 MT5 自己改参验证。
class="type">class="kw">double calcpe(class="type">int& mv[]) { class="type">int nmv = ArraySize(mv); if(nmv < class="num">1) class="kw">return class="num">0.0; class="type">int nup = class="num">0; for(class="type">int i = class="num">0; i < nmv; ++i) if(mv[i] > class="num">0) ++nup; class="kw">return ((class="type">class="kw">double)nup) / nmv; }
「用对数似然扫描变点分割样本」
上面这段把整段行情拆成离散涨跌步序列 mv[] 后,核心任务就是找出一个分割点 n1,使前后两段各自上涨概率 p1、p2 的极大似然估计最合理。外层循环从 n1=2 扫到 nmv-1,每次调用 llhx_n1 算当前分割下的对数似然值,只保留更大的那一组。 llhx_n1 内部先按 n1 把 mv 切成两份,分别数出上涨次数 nu1、nu2,再用 nu/n 算出 p1、p2,拼出标准二项对数似然 l = Σ nu·log(p) + (n-nu)·log(1-p)。当某段无上涨时跳过对应项,避免 log(0) 崩脚本。 OnStart 里写死了 tstart=D'2020.05.20'、tstop=D'2020.06.20'、dprcnt=0.5,跑完直接 Print 出 p1e、p2e。你在 MT5 里把这段挂上 Discr.mqh,改一下那三个 input 就能看不同月份 EURUSD 的涨跌结构是否出现概率断层——外汇和贵金属杠杆高,离散化结论只反映历史样本,实盘切换周期可能完全失效。 calcpes 另算了一种相邻步转移概率(nuu/nud/ndu/ndd),和外层变点扫描互补:前者找「哪天起规律变了」,后者看「涨后更易涨还是跌」。两者结合,才可能把一段行情从「随机游走」里抠出可交易的偏态。
class="type">int nmv = ArraySize(mv); if(nmv < class="num">2) class="kw">return; n1e = class="num">1; class="type">class="kw">double llhx = llhx_n1(mv, class="num">1, p1e, p2e), llh, p1, p2; for(class="type">int n1 = class="num">2; n1 < nmv; ++n1) { llh = llhx_n1(mv, n1, p1, p2); if(llh > llhx) { llhx = llh; n1e = n1; p1e = p1; p2e = p2; } } } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">// 对数似然函数最大值取决于 n1 class="type">class="kw">double llhx_n1(class="type">int& mv[], class="type">int n1, class="type">class="kw">double& p1, class="type">class="kw">double& p2) { p1 = p2 = class="num">0.0; class="type">int nmv = ArraySize(mv); if(nmv < class="num">2 || n1 < class="num">1 || n1 >= nmv) class="kw">return class="num">0.0; class="type">int nu1 = class="num">0, nu2 = class="num">0; for(class="type">int i = class="num">0; i < n1; ++i) if(mv[i] > class="num">0) ++nu1; for(class="type">int i = n1; i < nmv; ++i) if(mv[i] > class="num">0) ++nu2; class="type">class="kw">double l = class="num">0.0; if(nu1 > class="num">0) { p1 = ((class="type">class="kw">double)nu1) / n1; l += nu1 * log(p1); } if(nu1 < n1) l += (n1 - nu1) * log(((class="type">class="kw">double)(n1 - nu1)) / n1); if(nu2 > class="num">0) { p2 = ((class="type">class="kw">double)nu2) / (nmv - n1); l += nu2 * log(p2); } if(nu2 < nmv - n1) l += (nmv - n1 - nu2) * log(((class="type">class="kw">double)(nmv - n1 - nu2)) / (nmv - n1)); class="kw">return l; } class="macro">#include <Discr.mqh> class="macro">#class="kw">property script_show_inputs class=class="str">"cmt">//+------------------------------------------------------------------+ class="kw">input class="type">class="kw">datetime tstart = D&class="macro">#x27;class="num">2020.05.class="num">20 class="num">00:class="num">00&class="macro">#x27;; class="kw">input class="type">class="kw">datetime tstop = D&class="macro">#x27;class="num">2020.06.class="num">20 class="num">00:class="num">00&class="macro">#x27;; class="kw">input class="type">class="kw">double dprcnt = class="num">0.5; class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void OnStart() { class="type">int mv[]; setmv(mv, tstart, tstop, dprcnt); if(ArraySize(mv) < class="num">2) { Print("not enough moves"); class="kw">return; } class="type">class="kw">double p1e, p2e; calcpes(mv, p1e, p2e); Print("p1e=", p1e, " p2e=", p2e); } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">// 计算概率估值 class="type">void calcpes(class="type">int& mv[], class="type">class="kw">double& p1e, class="type">class="kw">double& p2e) { p1e = p2e = class="num">0; class="type">int nmv = ArraySize(mv); if(nmv < class="num">2) class="kw">return; class="type">int nuu = class="num">0, nud = class="num">0, ndu = class="num">0, ndd = class="num">0, nu, nd; for(class="type">int i = class="num">0; i < nmv - class="num">1; ++i) if(mv[i] > class="num">0) {