统计估计·进阶篇
（2/3）· 从异常值筛除到正态概率图，把教科书里的统计假设变成盘前可执行的诊断步骤
接上篇铺垫的随机序列概念，这一篇直接钻进估计环节。很多交易者把指标当黑箱，却从没确认过自己的价格样本到底是不是正态、有没有被极端跳空带偏，参数一错，后面模型全歪。
「四分位均值在MT5里的归一化写法」
MQL5里算完原始IQM后，常需要按样本规模做归一，否则不同品种或不同窗口长度之间数值不可比。上面这行把IQM除以(n-2*m)，本质是剔除上下各m个极端值后的有效样本量归一，让跨周期对比有意义。 在MT5里把这段直接贴进自定义指标，改n和m就能看不同分位区间的均值漂移；外汇与贵金属波动跳变多，极端值剔除比例调大一点，曲线会更抗噪，但滞后也会上来，属于高风险环境下的参数权衡。
IQM=IQM/(n-class="num">2*m); class=class="str">"cmt">// Interquartile mean(IQM)
用 Gnuplot 把 MT5 脚本结果画出来
想在 MT5 里把脚本算出的分布直接可视化，最省事的办法是借免费的 Gnuplot 画图。它能在独立窗口实时出图，也能按格式写文件，跨平台且带源码和编译版。本文示例只在 Windows XP SP3 下用 4.2.2 版跑过，安装包 gp442win32.zip 解压后会得到 \gnuplot 文件夹，必须整个丢进 MetaTrader 5 客户端根目录才能和脚本互动。 验证步骤很直接：进 \gnuplot\binary 跑 wgnuplot.exe，在 gnupplot> 提示符敲 plot sin(x)，弹窗画出正弦曲线就说明环境通了。之后启动 erremove.mq5，就会在单独窗口绘出分布图。 程序侧的交互极简：先写 gplot.txt 放 Gnuplot 命令和数据，再用 shell32.dll 导出的 ShellExecuteW() 调 wgnuplot.exe 并把这个文件当参数传进去，所以客户端要开「允许导入外部 dll」。 终端选型有个坑。wxt 终端抗锯齿、画质好，但关图窗口时 wgnuplot.exe 进程不会自动退；反复调用会堆出一堆僵尸进程。本文示例刻意用 windows 终端，命令行写 wgnuplot.exe -p MQL5\Files\gplot.txt，关窗口即杀进程，不留垃圾。图窗支持鼠标和键盘，活跃时按 h 能在文本窗列出默认热键。外汇与贵金属行情高波动，这类可视化仅用于辅助概率判断，不构成方向承诺。
◍ 用 h 估计量算清分布的四项特征
拿到一段行情采样，总体参数未知时只能靠无偏估计。采样平均值（数学期望估计）是分布中心，其他离差、标准方差、偏度、峰度都绕着它转。 峰度这块容易踩坑：正态分布的理论峰度是 3，但很多工具减掉 3 后报“超出峰度”。Excel 用 k 估计算超出量，而交易者统计书里常用 h 估计做无偏峰度，两者第四动差算法不同，表达式不一样，别混用。 下面 dStat() 直接按 h 估计落地。数组少于 4 个元素会返回 -1，因为中心动差到三阶以上样本量太小估计失效。中位数顺手排个序取中值，偶数长度取中间两数平均。 方差用 sum2/(n-1) 而非 n，偏度 skew=n*sum3/(n-2)/sum2/stdev，峰度那行乘了一串 (n-1)/(n-2)/(n-3) 修正项——这就是 h 估计的无偏处理。外汇和贵金属波动肥尾常见，算出来峰度显著大于 3 时，正态假设可能不成立，杠杆品种高风险。
class="kw">struct statParam { class="type">class="kw">double mean; class="type">class="kw">double median; class="type">class="kw">double var; class="type">class="kw">double stdev; class="type">class="kw">double skew; class="type">class="kw">double kurt; }; class=class="str">"cmt">//---------------------------------------------------------------------------- class="type">int dStat(const class="type">class="kw">double &x[],statParam &sP) { class="type">int i,m,n; class="type">class="kw">double a,b,sum2,sum3,sum4,y[]; ZeroMemory(sP); class=class="str">"cmt">// Reset sP n=ArraySize(x); if(n<class="num">4) class=class="str">"cmt">// Error { Print("Function dStat() error!"); class="kw">return(-class="num">1); } sP.kurt=class="num">1.0; ArrayResize(y,n); ArrayCopy(y,x); ArraySort(y); m=(n-class="num">1)/class="num">2; sP.median=y[m]; class=class="str">"cmt">// Median if((n&0x01)==class="num">0)sP.median=(sP.median+y[m+class="num">1])/class="num">2.0; sP.mean=class="num">0; for(i=class="num">0;i<n;i++)sP.mean+=x[i]; sP.mean/=n; class=class="str">"cmt">// Mean sum2=class="num">0;sum3=class="num">0;sum4=class="num">0; for(i=class="num">0;i<n;i++) { a=x[i]-sP.mean; b=a*a;sum2+=b; b=b*a;sum3+=b; b=b*a;sum4+=b; } if(sum2<class="num">1.e-150)class="kw">return(class="num">1); sP.var=sum2/(n-class="num">1); class=class="str">"cmt">// Variance sP.stdev=MathSqrt(sP.var); class=class="str">"cmt">// Standard deviation sP.skew=n*sum3/(n-class="num">2)/sum2/sP.stdev; class=class="str">"cmt">// Skewness sP.kurt=((n*n-class="num">2*n+class="num">3)*sum4/sum2/sum2-(class="num">6.0*n-class="num">9.0)/n)* (n-class="num">1.0)/(n-class="num">2.0)/(n-class="num">3.0); class=class="str">"cmt">// Kurtosis class="kw">return(class="num">1);
「用峰度算直方图分段数」
看采样分布不能只盯均值和标准差，有限样本的长相往往藏在直方图里。把数值范围切成若干段、数每段落进去的个数，再按段宽标准化，就得到经验分布密度的近似——这一步能直观判断采样是否贴近正态。 分段数不是拍脑袋定的。按交易者统计文献里的经验式，段数 L = (峰度+1.5) × N^0.4 / 6，N 是样本量、峰度来自已算好的 statParam。峰度越高尾巴越肥，分段就会自动变多，避免把厚尾误画成单峰。 下面的 dHist() 就是照这个公式落地的：传进来的 x[] 只读不改，histo[] 是动态数组接收结果，元素数 = 段数+2，头尾各补一个零值元素当边界。若 histo[] 不是动态数组、或 x[] 少于 4 个元素，函数直接返回 -1 并打印错误，不污染原序列。 别把正态当圣经 直方图画出来若明显比正态钟形胖尾或偏一边，外汇和贵金属这种高波动品种里往往意味着极端行情概率被低估，实盘前最好再叠一张正态概率图交叉看。
class="kw">struct statParam { class="type">class="kw">double mean; class="type">class="kw">double median; class="type">class="kw">double var; class="type">class="kw">double stdev; class="type">class="kw">double skew; class="type">class="kw">double kurt; }; class=class="str">"cmt">//---------------------------------------------------------------------------- class="type">int dHist(const class="type">class="kw">double &x[],class="type">class="kw">double &histo[],const statParam &sp) { class="type">int i,k,n,nbar; class="type">class="kw">double a[],max,s,xmin; if(!ArrayIsDynamic(histo)) class=class="str">"cmt">// Error { Print("Function dHist() error!"); class="kw">return(-class="num">1); } n=ArraySize(x); if(n<class="num">4) class=class="str">"cmt">// Error { Print("Function dHist() error!"); class="kw">return(-class="num">1); } nbar=(sp.kurt+class="num">1.5)*MathPow(n,class="num">0.4)/class="num">6.0; if((nbar&0x01)==class="num">0)nbar--; if(nbar<class="num">5)nbar=class="num">5; class=class="str">"cmt">// Number of bars ArrayResize(a,n); ArrayCopy(a,x); max=class="num">0.0; for(i=class="num">0;i<n;i++) { a[i]=(a[i]-sp.mean)/sp.stdev; class=class="str">"cmt">// Normalization if(MathAbs(a[i])>max)max=MathAbs(a[i]); } xmin=-max; s=class="num">2.0*max*n/nbar; ArrayResize(histo,nbar+class="num">2); ArrayInitialize(histo,class="num">0.0); histo[class="num">0]=class="num">0.0;histo[nbar+class="num">1]=class="num">0.0; for(i=class="num">0;i<n;i++) { k=(a[i]-xmin)/max/class="num">2.0*nbar; if(k>(nbar-class="num">1))k=nbar-class="num">1; histo[k+class="num">1]++; } for(i=class="num">0;i<nbar;i++)histo[i+class="num">1]/=s; class="kw">return(class="num">1); }
用 Rankit 法把收益序列拉成一条直线看正态
正态概率图的核心目的，是把样本序列做拉伸变换后画到坐标里——若价格收益率真服从正态，点会近似落在一条直线上。肉眼扫偏离程度，比跑一个 p 值更适合盯盘时快速判断分布假设是否成立。 计算这类图的坐标，靠的是 dRankit() 函数。它吃进原始数组 x[]，把 Y 轴用的标准化值和 X 轴用的理论分位值分别写进 resp[] 与 xscale[]，同时要求传入已算好的 statParam 结构（均值、标准差等）。样本数少于 4 会直接返回 -1，实盘里小样本画这图没意义。 X 轴分位值由 ltqnorm() 算逆正态累积分布得到，算法源自一篇计算逆 CDF 的论文。下面这段代码里，xscale 首尾用 0.5^(1/n) 对称处理，中间走 (i+1-0.3175)/(n+0.365) 的 Rankit 近似——换不同系数（比如 Blom 用 3/8）会让尾部点轻微偏移，MT5 里改那两个常数就能比对。 外汇与贵金属杠杆高，分布厚尾是常态，正态假设被拒是概率事件而非异常，别拿直线近似当风控依据。
class="kw">struct statParam { class="type">class="kw">double mean; class="type">class="kw">double median; class="type">class="kw">double var; class="type">class="kw">double stdev; class="type">class="kw">double skew; class="type">class="kw">double kurt; }; class=class="str">"cmt">//---------------------------------------------------------------------------- class="type">int dRankit(const class="type">class="kw">double &x[],class="type">class="kw">double &resp[],class="type">class="kw">double &xscale[],const statParam &sp) { class="type">int i,n; class="type">class="kw">double np; if(!ArrayIsDynamic(resp)||!ArrayIsDynamic(xscale)) class=class="str">"cmt">// Error { Print("Function dHist() error!"); class="kw">return(-class="num">1); } n=ArraySize(x); if(n<class="num">4) class=class="str">"cmt">// Error { Print("Function dHist() error!"); class="kw">return(-class="num">1); } ArrayResize(resp,n); ArrayCopy(resp,x); ArraySort(resp); for(i=class="num">0;i<n;i++)resp[i]=(resp[i]-sp.mean)/sp.stdev; ArrayResize(xscale,n); xscale[n-class="num">1]=MathPow(class="num">0.5,class="num">1.0/n); xscale[class="num">0]=class="num">1-xscale[n-class="num">1]; np=n+class="num">0.365; for(i=class="num">1;i<(n-class="num">1);i++)xscale[i]=(i+class="num">1-class="num">0.3175)/np; for(i=class="num">0;i<n;i++)xscale[i]=ltqnorm(xscale[i]); class="kw">return(class="num">1); } class=class="str">"cmt">//---------------------------------------------------------------------------- class="type">class="kw">double A1 = -class="num">3.969683028665376e+01, A2 = class="num">2.209460984245205e+02, A3 = -class="num">2.759285104469687e+02, A4 = class="num">1.383577518672690e+02, A5 = -class="num">3.066479806614716e+01, A6 = class="num">2.506628277459239e+00; class="type">class="kw">double B1 = -class="num">5.447609879822406e+01, B2 = class="num">1.615858368580409e+02, B3 = -class="num">1.556989798598866e+02, B4 = class="num">6.680131188771972e+01, B5 = -class="num">1.328068155288572e+01; class="type">class="kw">double C1 = -class="num">7.784894002430293e-03, C2 = -class="num">3.223964580411365e-01, C3 = -class="num">2.400758277161838e+00, C4 = -class="num">2.549732539343734e+00, C5 = class="num">4.374664141464968e+00, C6 = class="num">2.938163982698783e+00; class="type">class="kw">double D1 = class="num">7.784695709041462e-03, D2 = class="num">3.224671290700398e-01,
◍ 分尾区用对数变换兜底
中心区之外（p<0.02425 或 p>0.97575）直接套有理式会飘，代码改用 sqrt(-2*log(p)) 做自变量变换再喂进另一组系数。 上下尾的区别只在一个符号：下尾 s=1、上尾 s=-1，其余共用 C1~C6 与 D1~D4 的分子分母结构。D3 写死为 2.445134137142996、D4 为 3.754408661907416，是这组近似的固定常数。 把 ltqnorm 贴进 MT5 脚本，传入 0.001 和 0.999 各跑一次，返回约 -3.09 与 +3.09；外汇与贵金属波动带常落在这类极端分位，用前先认清杠杆下的高风险。
class="type">class="kw">double ltqnorm(class="type">class="kw">double p) { class="type">int s=class="num">1; class="type">class="kw">double r,x,q=class="num">0; if(p<=class="num">0||p>=class="num">1){Print("Function ltqnorm() error!");class="kw">return(class="num">0);} if((p>=class="num">0.02425)&&(p<=class="num">0.97575)) class=class="str">"cmt">// Rational approximation for central region { q=p-class="num">0.5; r=q*q; x=(((((A1*r+A2)*r+A3)*r+A4)*r+A5)*r+A6)*q/(((((B1*r+B2)*r+B3)*r+B4)*r+B5)*r+class="num">1); class="kw">return(x); } if(p<class="num">0.02425) class=class="str">"cmt">// Rational approximation for lower region { q=sqrt(-class="num">2*log(p)); s=class="num">1; } else class=class="str">"cmt">//if(p>class="num">0.97575) // Rational approximation for upper region { q = sqrt(-class="num">2*log(class="num">1-p)); s=-class="num">1; } x=s*(((((C1*q+C2)*q+C3)*q+C4)*q+C5)*q+C6)/((((D1*q+D2)*q+D3)*q+D4)*q+class="num">1); class="kw">return(x); }