未知概率密度函数的核密度估计·进阶篇
（2/3）· 面对 10 到 10000 个值的非平稳行情片段，正态假设为何会误导你的判断
◍ 用核密度把随机序列估出概率分布
在 MT5 里做非参数密度估计，最省事的办法是借 CDens 类直接算。下面这段把输入先填成均匀随机量，再交给核密度过程跑出平滑的 PDF，外汇与贵金属价格序列虽非均匀，但同样可用此框架替换 X[] 来源做分布体检，注意杠杆品种回测失败率可能偏高。 代码先以 MathRand() 造 ndata 个原始样本，NTpoints 设测试点数量，Density 以带宽 0.22 做估计；随后把测试坐标 T[] 与密度值 Y[] 拷出，类实例即释放。
for(i=class="num">0;i<ndata;i++)X[i]=MathRand();class=class="str">"cmt">// Create class="kw">input data CDens *kd=new CDens; kd.NTpoints(npoint); class=class="str">"cmt">// Setting number of test points kd.Density(X,class="num">0.22); class=class="str">"cmt">// Density estimation with h=class="num">0.22 ArrayCopy(T,kd.T); class=class="str">"cmt">// Copy test points ArrayCopy(Y,kd.Y); class=class="str">"cmt">// Copy result(pdf) class="kw">delete(kd); class=class="str">"cmt">// Result: T[]-test points, Y[]-density estimation } class=class="str">"cmt">//--------------------------------------------------------------------
for(i=class="num">0;i<ndata;i++)X[i]=MathRand();class=class="str">"cmt">// Create class="kw">input data CDens *kd=new CDens; kd.NTpoints(npoint); class=class="str">"cmt">// Setting number of test points kd.Density(X,class="num">0.22); class=class="str">"cmt">// Density estimation with h=class="num">0.22 ArrayCopy(T,kd.T); class=class="str">"cmt">// Copy test points ArrayCopy(Y,kd.Y); class=class="str">"cmt">// Copy result(pdf) class="kw">delete(kd); class=class="str">"cmt">// Result: T[]-test points, Y[]-density estimation } class=class="str">"cmt">//--------------------------------------------------------------------
带宽 h 选错，密度估计就废了
图 1 用 CDens 类跑了不同 h 的密度估计，叠了真实高斯曲线做对照。h=0.22 时估计最贴真值；h 再大就过度平滑，再小则平滑不足，曲线明显发散或移位。 Silverman 拇指法则够轻量：取序列标准差与四分位差/1.34 的较小值作 A，再乘 0.9/N^0.2。原文序列标准差为 1，所以 A 直接取 min(IQR/1.34, 1)，代码里一句话就算完。 但正态假设一破，这法子就偏。好在它看四分位差能把 h 往小调，当个初始值够用，不烧 CPU。 真要准，Sheather-Jones 插件法（SJPI）更顶。它用高斯核（无限支集），直接算要 O(N*N) 次操作，SJPI 还得先估两轮密度导数、再用 Newton-Raphson 逼零，耗时约是普通估计的十倍。 加速靠三招：自变量超 4 的高斯核值直接舍（图上看不出差别）；估计点 M 从 N 降到 200（N=100000 时操作数从 O(N*N) 掉到 O(M*N)）；再不济上快速高斯变换或分箱。外汇和贵金属价格序列高波动、非正态倾向强，核密度带宽选错会误导形态判断，实盘前务必在 MT5 里跑一遍验证。
ArraySort(X); i=(class="type">int)((N-class="num">1.0)/class="num">4.0+class="num">0.5); a=(X[N-class="num">1-i]-X[i])/class="num">1.34; class=class="str">"cmt">// IQR/class="num">1.34 a=MathMin(a,class="num">1.0); h=class="num">0.9*a/MathPow(N,class="num">0.2); class=class="str">"cmt">// Silverman&class="macro">#x27;s rule of thumb
「SJPI 插件怎么算出更细的带宽」
CSJPlugin 类只暴露一个公共方法：double SelectH(double &px[],double h,double stdev=1)。调用时传原始序列引用、初始 h 和标准化标准差，成功返回最优 h，失败则退回传入的初始 h。 初始 h 越接近真实解越好，通常用 Silverman 拇指法则先算一个垫底值，能压住 Newton-Raphson 跑偏的概率。序列在类里先被压到 [0,1]，h 随尺度正态化，算完再反正态化回去。 构造函数里写死了 eps=1e-4（密度与导数精度）、P_UL=500（最大迭代）；SelectH 内 IMAX=20 限制 Newton-Raphson 步数，PREC=1e-4 控方程求解精度。导数不走解析，而是自变量加个小增量做估计，省掉每次迭代的求导开销。 拿 10000 个元素的随机序列实测：Silverman 给 h=0.14，SJPI 给 h=0.07。图 2 里红色是真密度，h=0.07 时尖顶塑形更利，左倾凹陷处离散没明显放大，比 0.14 更接近真相。 SJPI 不是万能。长序列即使走了快速算法仍可能磨时间；短序列（10–30 点）它还容易把 h 估大，超过 Silverman。落地规则：默认用 Silverman，SJPI 需显式开启，且最终 h 取两种方法的最小值来兜住高估风险。外汇与贵金属波动剧烈，带宽选错会放大噪声，回测前务必在 MT5 用真实 tick 复算一遍。
◍ 核密度估计里的边界裁剪与反射校正
用 Epanechnikov 这类有限支集核时，核中心一旦靠近数据序列边界，能覆盖的样本最多掉到 50%，估计值会明显移位、离差放大。即便换高斯核（无限支集），自变量大了核值近乎 0，边界效应表现和有限支集核一样。 拿平坦序列 X=1,2,3,…,n 试一下最直观：理论密度应是一条水平直线，但原始 kdens() 画出来在边界处塌陷出明显凹陷（见图 4）。四种常规解法——数据反射、数据变换、伪数据、边界核——里挑了反射法，原因是它不产出负密度、数学轻量，只是总运算量因序列虚增而上升。 实现上没去物理扩数组，而是在循环里用镜像坐标 c、d 把边界外的虚拟点算进来。下面对带反射的 kdens() 逐行拆： hh=h/MathSqrt(0.5) —— 带宽按反射规则缩放； s=sqrt(M_PI+M_PI)*N*h —— 归一化分母，含反射后的等效样本量； c=(X[0]+X[0])/hh 与 d=(X[N-1]+X[N-1])/hh —— 左右边界点的镜像中心，用于后续 g=c-e-b、g=d-e-b 的虚拟样本； 循环里对每个网格点 e=T[i]/hh，先算与首、尾真实点的高斯项，再遍历内部点 j 并额外加两次反射项； 所有项受 g∈(-3,3) 截断，省掉尾部近零计算； Y[i]=a/s 即输出概率密度。 校正后平坦序列的边界凹陷消失，成了一条直线。但图 6 的线性增长分布下，蓝色反射版仍没把边界移位吃干净——密度在边界陡升/陡降时，反射法只让曲线稍微平缓，无法全消。外汇与贵金属价格序列边界常带趋势梯度，用此法做密度分析要留这个偏误余地，属高风险品种上的近似手段。
class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Kernel density estimation with reflection of data | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void kdens(class="type">class="kw">double h) { class="type">int i,j; class="type">class="kw">double a,b,c,d,e,g,s,hh; hh=h/MathSqrt(class="num">0.5); s=sqrt(M_PI+M_PI)*N*h; c=(X[class="num">0]+X[class="num">0])/hh; d=(X[N-class="num">1]+X[N-class="num">1])/hh; for(i=class="num">0;i<Np;i++) { e=T[i]/hh; a=class="num">0; g=e-X[class="num">0]/hh; if(g>-class="num">3&&g<class="num">3)a+=MathExp(-g*g); g=e-X[N-class="num">1]/hh; if(g>-class="num">3&&g<class="num">3)a+=MathExp(-g*g); for(j=class="num">1;j<N-class="num">1;j++) { b=X[j]/hh; g=e-b; if(g>-class="num">3&&g<class="num">3)a+=MathExp(-g*g); g=d-e-b; if(g>-class="num">3&&g<class="num">3)a+=MathExp(-g*g); g=c-e-b; if(g>-class="num">3&&g<class="num">3)a+=MathExp(-g*g); } Y[i]=a/s; class=class="str">"cmt">// pdf } }
给密度估计类补上自动带宽与边界修正
在之前 CDens 的基础上，把自动选 h 区间和边界效应校正塞进同一个类，能让 KDE 跑起来更省心。核心入口是 Density(double &x[], double hh=-1)：x[] 是待分析序列，长度不能低于 8；hh 留负值就自动挑范围，给正数则强制使用该带宽，但任何情况下建议把 h 压在 0.005 以内，否则平滑过头会吃掉价格行为的细节。 方法跑完返回 0 代表正常，负数代表中途出错。成功后 T[] 存检验点、Y[] 存对应密度估计值，X[] 则是被标准化并排序过的输入——均值归零、标准差为 1，方便跨品种对比分布形态。 类里暴露的变量值得盯一下：N 是样本量，Np 是检验点数量（默认 200，可用 NTpoints(int n) 在调 Density 前改，且不能少于 10），H 是最终带宽，Pflag 标记是否用了 SJPI 法求最优 h。PluginMode(int m) 传 1 就开 SJPI，不传或传别的默认关闭；外汇和贵金属波动聚集明显，关掉自动选宽往往更容易看出极端区的真实密度，但这类品种杠杆高、回撤快，任何分布结论都只是概率倾向，别当确定性信号。
class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| CKDensity.mqh | class=class="str">"cmt">//| class="num">2012, victorg | class=class="str">"cmt">//| [MQL5官方文档] | class=class="str">"cmt">//+------------------------------------------------------------------+ class="macro">#class="kw">property copyright "class="num">2012, victorg" class="macro">#class="kw">property link "[MQL5官方文档] class="macro">#include <Object.mqh> class="macro">#include "CSJPlugin.mqh" class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Class Kernel Density Estimation | class=class="str">"cmt">//+------------------------------------------------------------------+ class CKDensity:class="kw">public CObject { class="kw">public: class="type">class="kw">double X[]; class=class="str">"cmt">// Data class="type">int N; class=class="str">"cmt">// Input data length(N >= class="num">8) class="type">class="kw">double T[]; class=class="str">"cmt">// Test points for pdf estimating class="type">class="kw">double Y[]; class=class="str">"cmt">// Estimated density(pdf) class="type">int Np; class=class="str">"cmt">// Number of test points(Npoint>=class="num">10, class="kw">default class="num">200) class="type">class="kw">double Mean; class=class="str">"cmt">// Mean(average) class="type">class="kw">double Var; class=class="str">"cmt">// Variance class="type">class="kw">double StDev; class=class="str">"cmt">// Standard deviation class="type">class="kw">double H; class=class="str">"cmt">// Bandwidth class="type">int Pflag; class=class="str">"cmt">// SJ plug-in bandwidth selection flag class="kw">public: class="type">void CKDensity(class="type">void); class="type">int Density(class="type">class="kw">double &x[],class="type">class="kw">double hh=-class="num">1);
「核密度估计的构造与默认参数落地」
这个类把核密度估计（KDE）的初始化和带宽选择拆成了几个清晰的方法。构造函数 CKDensity() 里直接把测试点数设为 200、插件模式标志 Pflag 设为 0，意味着默认不做 Sheather-Jones 插件修正，直接用基础核估计跑分布。 NTpoints(int n) 负责设定估算用的测试网格规模：传入值小于 10 会被强制抬到 10，随后用 ArrayResize 给测试点数组 T 和密度结果数组 Y 各分配 Np 长度。你在 MT5 里若想看更细的密度曲线，把 n 调大即可，但计算量会随点数线性增加。 Density() 是主入口，先取输入样本长度 N，若小于 8 直接 Print 报错并返回 -1——样本太少时核估计方差会爆炸，外汇分钟线至少攒够 30 根以上再跑才有点参考价值。它把数据拷进 X 并排序，Mean 归零留给后续均值计算。贵金属与外汇杠杆品种波动聚集明显，KDE 只是刻画历史分布形态，对跳空和极端尾风险解释力有限，实操须自担高风险。
class="type">void NTpoints(class="type">int n); class="type">void PluginMode(class="type">int m) {if(m==class="num">1)Pflag=class="num">1; else Pflag=class="num">0;} class="kw">private: class="type">void kdens(class="type">class="kw">double h); }; class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Constructor | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void CKDensity::CKDensity(class="type">void) { NTpoints(class="num">200); class=class="str">"cmt">// Default number of test points Pflag=class="num">0; class=class="str">"cmt">// Not SJ plug-in } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Setting number of test points | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void CKDensity::NTpoints(class="type">int n) { if(n<class="num">10)n=class="num">10; Np=n; class=class="str">"cmt">// Number of test points ArrayResize(T,Np); class=class="str">"cmt">// Array for test points ArrayResize(Y,Np); class=class="str">"cmt">// Array for result(pdf) } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Bandwidth selection and kernel density estimation | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">int CKDensity::Density(class="type">class="kw">double &x[],class="type">class="kw">double hh=-class="num">1) { class="type">int i; class="type">class="kw">double a,b,h; N=ArraySize(x); class=class="str">"cmt">// Input data length if(N<class="num">8) class=class="str">"cmt">// If N is too small { Print(__FUNCTION__+": Error! Not enough data length!"); class="kw">return(-class="num">1); } ArrayResize(X,N); class=class="str">"cmt">// Array for class="kw">input data ArrayCopy(X,x); class=class="str">"cmt">// Copy class="kw">input data ArraySort(X); class=class="str">"cmt">// Sort class="kw">input data Mean=class="num">0;