MetaTrader 5 / 统计分析
English Русский Español Deutsch 日本語 Português
preview
将互信息作为渐进特征选择的准则

将互信息作为渐进特征选择的准则

MetaTrader 5 — 统计分析 |
219 0
Francis Dube
Francis Dube

概述

互信息是识别有效预测变量的有力工具，特别是在处理复杂的非线性关系时。它能够发现其他方法可能遗漏的依赖关系，因此特别适用于能够利用此类复杂关系的模型。本文探讨了互信息在特征选择中的应用，重点关注了Hanchuan Peng、Fuhui Long和Chris Ding在其研究论文《基于互信息的特征选择：最大依赖性、最大相关性和最小冗余性准则》中提出的算法。

我们将首先讨论连续变量互信息的估计方法，然后深入探讨特征选择过程本身。最后，我们将通过涉及合成数据集和真实数据集的示例来说明该算法的有效性。


连续变量互信息的估计

互信息，记为I(X;Y)，用于量化两个随机变量X和Y之间的共享信息。对于离散变量，其计算相对简单，涉及对联合概率和边缘概率求和。

离散互信息

然而，由于连续变量的可能取值范围是无限的，因此它们带来了重大挑战。解决这一问题的一种常见方法是将连续变量离散化为若干区间（bins），然后应用离散互信息公式。遗憾的是，该方法的准确性对所选区间宽度或区间数量高度敏感，这可能导致估计的互信息出现较大波动。

连续互信息

要计算连续变量的互信息，我们必须从离散求和过渡到连续积分。这需要了解联合概率密度函数（PDFs）和边缘概率密度函数，但由于数据有限，这些函数通常无法直接获得。帕尔森窗（Parzen window）技术通过直接从数据中近似概率密度函数来提供解决方案。该方法涉及在数据上放置一个窗函数，并计算窗内数据点的加权和。分配给每个数据点的权重随其与窗中心距离的增加而减小。通过在数据空间中移动此窗并计算每个点的加权和，我们可以构建概率密度函数的估计值。

帕尔森密度近似

在底层分布未知或复杂的情况下，该技术特别有用。通过调整窗的形状和大小，我们可以调整密度估计的平滑度和准确性。然而，值得注意的是，窗参数的选择会显著影响估计的质量。帕尔森密度估计的有效性取决于选择合适的加权函数，该函数理想情况下应在整个定义域上的积分等于1，并且当自变量偏离0时迅速衰减。该加权函数的一个常用选择是高斯函数，其特征是具有一个尺度参数σ。

高斯加权函数

然而，确定σ的最优值是一项颇具挑战性的任务，往往会导致估计的密度出现较大波动。帕尔森窗方法对σ的选择如此敏感，这是该方法的一大缺点，使得它在诸如评估预测变量候选等任务中不太理想，因为这些任务要求密度估计必须准确无误。

除了帕尔森窗方法，自适应分区（Adaptive Partitioning）是另一种用于估计连续变量间互信息的方法。与固定分箱的简单方法相比，自适应分区有了显著的改进。通过迭代细分信息含量高的区域，自适应分区能有效地将计算资源集中在数据空间中最相关的区域。这种有针对性的方法能够更准确地估计互信息，尤其是在处理复杂、非均匀分布时。自适应分区的一个关键优势在于其能够适应底层数据结构。递归分区过程确保具有显著依赖关系的区域得到进一步细分，而信息量极少的区域则保持不变。这种动态方法避免了固定分箱的常见陷阱，即可能对数据进行过度平滑或过度拟合。

此外，自适应分区还融入了统计检验来指导分区过程。通过评估潜在细分的统计显著性，自适应分区能够区分真实信息与随机噪声。这有助于防止过度拟合，否则可能导致互信息的估计值虚高。原始方法的主要问题在于，它过度关注双变量域中案例稀少或没有案例的区域，而可能忽略了大部分信息所在的密集区域。  

自适应分区

自适应分区从数据空间的粗粒度划分开始。对于每个分区，计算变量间的互信息。如果该值超过预设阈值，则将该分区进一步细分为更小的子分区。这一递归分区过程持续进行，直至满足停止准则。若分区过程过早停止，所得分区可能过于粗略，导致细节丢失，并使估计的互信息产生向下偏差。反之，若分区过程持续过久，分区可能变得过于精细，导致过拟合，且因噪声而使互信息估计值虚高。

为确定分区是否需进一步细分，应采用卡方检验。该检验用于评估分区内两个变量之间的独立性。将分区划分为四个子分区，形成一个2x2列联表。统计落入四个子分区中的每个分区的数据点数量。在两个变量独立的原假设下，根据边缘总数计算每个子分区中数据点的期望数量。然后，将计算得到的卡方统计量与自由度为1的卡方分布临界值进行比较。若计算值超过临界值，则拒绝独立性的原假设，表明分区内两个变量之间存在显著关系。

若卡方检验结果显著，则将分区进一步细分为更小的分区。此过程递归进行，直至分区变得足够均匀，或满足停止准则。若卡方检验结果不显著，则认为分区均匀，无需进一步细分。

为应对分区内可能存在的更复杂关系（而简单的2x2分割可能无法捕捉到这种关系），该算法还纳入了额外的检查步骤。若2x2卡方检验未能检测到显著关系，但分区仍然相对较大，则进行更精细的4x4分割。然后，对该4x4分区应用卡方检验，以评估是否存在非随机分布。若4x4检验也未能检测到显著关系，则认为分区均匀，并计算其对整体互信息的贡献。在此情况下，无需进一步细分。

若2x2或4x4卡方检验表明分区内存在非均匀分布，则将该分区细分为四个更小的子分区。这些子分区的处理方式因大小而异。若子分区过小，则视为终节点。计算其对整体互信息的贡献，且不再考虑该分区进行进一步细分。否则，若子分区足够大，则将其标记为未来细分的候选分区。

标准的2x2卡方检验统计量如下式所示。

2乘2卡方检验统计量


MQL5中的自适应分区

mutualinfo.mqh中的Capm类实现了自适应分区算法，用于估计连续变量的互信息。为了管理递归分区过程，使用了自定义结构体IntStack。该结构体存储有关分区及其所包含的任何子分区的信息。该结构体有六个成员，具体如下：

  • Xstart和Xstop：这两个成员定义了当前包围分区内一个变量的索引范围。
  • Ystart和Ystop：这两个成员定义了当前包围分区内第二个变量的索引范围。
  • DataStart和DataStop：这两个成员指定了包围分区内数据点的起始和结束索引，用于标记子分区。
//+------------------------------------------------------------------+
//| int type stack structure                                         |
//+------------------------------------------------------------------+
struct IntStack
  {
   int Xstart ;     // X value (rank) at which this rectangle starts
   int Xstop ;      // And stops
   int Ystart ;     // Ditto for Y
   int Ystop ;
   int DataStart ;  // Starting index into indices for the cases in this
   int DataStop ;   // rectangle, and the (inclusive) ending index
  };

Capm构造函数接受三个参数：

  • no_split：一个布尔标志，用于确定如何处理近似相等的值。如果设置为 true，算法会在分区过程中确保将这些值分组在一起，防止它们被分到不同的分区中。
  • dep_vals：一个向量，包含待分析的其中一个变量的值。该变量通常是用于与众多其他变量计算互信息的变量。
  • chi_critical_value：分区过程中使用的卡方检验的阈值。该值决定检验的显著性水平，并影响分区的深度。一般而言，取值在4.0到8.0之间较为合适，其中6.0是一个普遍适用且常用的选择。
//+------------------------------------------------------------------+
//|  constructor                                                     |
//+------------------------------------------------------------------+
Capm::Capm(bool no_split,vector &dep_vals, double chi_critical_value = 6.0)
  {
   m_no_split = no_split;
   m_size = int(dep_vals.Size()) ;
   m_chi_crit = chi_critical_value;

   if(ArrayResize(m_indices,m_size)<0 || !m_tempvals.Resize(m_size) || ArrayResize(m_ranks,m_size)<0
      || (m_no_split && ArrayResize(m_tied_ranks,m_size)<0))
     {
      Print(__FUNCTION__, " Arrayresize error ", GetLastError());
      return;
     }

   for(int i=0 ; i<m_size ; i++)
     {
      m_tempvals[i] = dep_vals[i] ;
      m_indices[i] = i ;
     }

   qsortdsi(0, m_size-1, m_tempvals, m_indices) ;

   for(int i=0 ; i<m_size ; i++)
     {
      m_ranks[m_indices[i]] = i ;
      if(! m_no_split)
         continue ;
      if(i < m_size-1  &&
         m_tempvals[i+1] - m_tempvals[i] < 1.e-12 * (1.0+fabs(m_tempvals[i])+fabs(m_tempvals[i+1])))
         m_tied_ranks[i] = 1 ;
      else
         m_tied_ranks[i] = 0 ;
     }
  }

Capm类构造函数首先会设置no_split标识（该标识决定是否考虑数据中的并列值情况）和chi_critical_value（该值设定用于识别变量间显著关系的卡方检验阈值）。接下来，构造函数会为存储数据、排序结果及并列值的数组分配内存。它将变量值复制到一个临时向量中并对其进行排序。数据点的原始索引则保存在m_indices数组中。最后，构造函数会为排序后的数据点分配排名，并识别任何并列值；如果设置了no_split标识，则会在m_tied_ranks数组中对这些并列值进行标记。

//+------------------------------------------------------------------+
//| Mutual information using the Adaptive partitioning method        |
//+------------------------------------------------------------------+
class Capm
  {

public:
                     Capm(bool no_split,vector &dep_vals, double chi_critical_value = 6.0) ;
                    ~Capm(void) ;
   double            fit(vector &raw) ;

private:
   bool               m_no_split;  // respect ties
   int m_size ;                        // Number of cases
   int m_ranks[] ;                     // 'Dependent' variable ranks
   int m_tied_ranks[] ;                // tied[i] != 0 if case with rank i == case with rank i+1
   double m_chi_crit ;                 // Chi-square test criterion
   int               m_indices[];      // indices
   vector            m_tempvals;       // temp values
  } ;

Capm类中的fit()方法旨在估计一个变量（在构造函数中提供）与作为参数传入的另一个变量之间的互信息。它接收一个向量作为输入，该向量代表第二个变量的取值。此向量的维度应与提供给构造函数的向量维度相同。通过多次调用fit()方法并传入不同的向量，我们可以高效地计算该变量与多个不同变量之间的互信息。

//+------------------------------------------------------------------+
//|  mutual information for continuous variables                     |
//+------------------------------------------------------------------+
double Capm::fit(vector &raw)
  {
   int i, k, ix, iy, nstack, splittable ;
   int current_indices[], x[], fullXstart, fullXstop, fullYstart, fullYstop, ipos ;
   int trialXstart[4], trialXstop[4], trialYstart[4], trialYstop[4] ;
   int ipx, ipy, xcut[4], ycut[4], iSubRec, x_tied[], ioff ;
   int X_AllTied, Y_AllTied ;
   int centerX, centerY, currentDataStart, currentDataStop ;
   int actual[4], actual44[16] ;
   vector expected(16), xfrac(4), yfrac(4) ;
   double diff, testval;
   double px, py, pxy, MI ;

   IntStack stack[];

   if(ArrayResize(stack,1,256)<0)
     {
      Print(__FUNCTION__," arrayresize error ", GetLastError());
      return EMPTY_VALUE;
     }

   if(ArrayResize(current_indices,m_size)<0 || ArrayResize(x,m_size)<0
      || (m_no_split && ArrayResize(x_tied,m_size)<0))
     {
      Print(__FUNCTION__, " Arrayresize error ", GetLastError());
      return EMPTY_VALUE;
     }

   for(i=0 ; i<m_size ; i++)
     {
      m_tempvals[i] = raw[i] ;
      m_indices[i] = i ;
     }

   if(!qsortdsi(0, m_size-1, m_tempvals, m_indices))
     {
      Print(__FUNCTION__, " Arrayresize error ", GetLastError());
      return EMPTY_VALUE;
     }

   for(i=0 ; i<m_size ; i++)
     {
      x[m_indices[i]] = i ;
      if(! m_no_split)
         continue ;
      if(i < m_size-1  &&
         m_tempvals[i+1] - m_tempvals[i] < 1.e-12 * (1.0+fabs(m_tempvals[i])+fabs(m_tempvals[i+1])))
         x_tied[i] = 1 ;
      else
         x_tied[i] = 0 ;
     }

   for(i=0 ; i<m_size ; i++)
     {
      m_indices[i] = i ;
     }

   stack[0].Xstart = 0 ;
   stack[0].Xstop = m_size-1 ;
   stack[0].Ystart = 0 ;
   stack[0].Ystop = m_size-1 ;
   stack[0].DataStart = 0 ;
   stack[0].DataStop = m_size-1 ;
   nstack = 1 ;

   MI = 0.0 ;
   while(nstack > 0)
     {

      --nstack ;
      fullXstart = stack[nstack].Xstart ;
      fullXstop = stack[nstack].Xstop ;
      fullYstart = stack[nstack].Ystart ;
      fullYstop = stack[nstack].Ystop ;
      currentDataStart = stack[nstack].DataStart ;
      currentDataStop = stack[nstack].DataStop ;


      centerX = (fullXstart + fullXstop) / 2 ;
      X_AllTied = (x_tied.Size())  && (x_tied[centerX] != 0) ;
      if(X_AllTied)
        {
         for(ioff=1 ; centerX-ioff >= fullXstart ; ioff++)
           {
            if(! x_tied[centerX-ioff])
              {
               X_AllTied = 0 ;
               centerX -= ioff ;
               break ;
              }
            if(centerX + ioff == fullXstop)
               break ;
            if(! x_tied[centerX+ioff])
              {
               X_AllTied = 0 ;
               centerX += ioff ;
               break ;
              }
           }
        }

      centerY = (fullYstart + fullYstop) / 2 ;
      Y_AllTied = (m_tied_ranks.Size())  && (m_tied_ranks[centerY] != 0) ;
      if(Y_AllTied)
        {
         for(ioff=1 ; centerY-ioff >= fullYstart ; ioff++)
           {
            if(! m_tied_ranks[centerY-ioff])
              {
               Y_AllTied = 0 ;
               centerY -= ioff ;
               break ;
              }
            if(centerY + ioff == fullYstop)
               break ;
            if(! m_tied_ranks[centerY+ioff])
              {
               Y_AllTied = 0 ;
               centerY += ioff ;
               break ;
              }
           }
        }

      if(X_AllTied  ||  Y_AllTied)
         splittable = 0 ;
      else
        {
         trialXstart[0] = trialXstart[1] = fullXstart ;
         trialXstop[0] = trialXstop[1] = centerX ;
         trialXstart[2] = trialXstart[3] = centerX+1 ;
         trialXstop[2] = trialXstop[3] = fullXstop ;
         trialYstart[0] = trialYstart[2] = fullYstart ;
         trialYstop[0] = trialYstop[2] = centerY ;
         trialYstart[1] = trialYstart[3] = centerY+1 ;
         trialYstop[1] = trialYstop[3] = fullYstop ;

         for(i=0 ; i<4 ; i++)
            expected[i] = (currentDataStop - currentDataStart + 1) *
                          (trialXstop[i]-trialXstart[i]+1.0) / (fullXstop-fullXstart+1.0) *
                          (trialYstop[i]-trialYstart[i]+1.0) / (fullYstop-fullYstart+1.0) ;

         actual[0] = actual[1] = actual[2] = actual[3] = 0 ;
         for(i=currentDataStart ; i<=currentDataStop ; i++)
           {
            k = m_indices[i] ;
            if(x[k] <= centerX)
              {
               if(m_ranks[k] <= centerY)
                  ++actual[0] ;
               else
                  ++actual[1] ;
              }
            else
              {
               if(m_ranks[k] <= centerY)
                  ++actual[2] ;
               else
                  ++actual[3] ;
              }
           }

         testval = 0.0 ;
         for(i=0 ; i<4 ; i++)
           {
            diff = fabs(actual[i] - expected[i]) - 0.5 ;
            testval += diff * diff / expected[i] ;
           }

         splittable = (testval > m_chi_crit)  ?  1 : 0 ;

该方法首先初始化一个栈，用于管理分区过程。最初，整个数据集作为根分区被压入栈中。随后，算法从栈中迭代地弹出一个分区，并评估其同质性（均匀性）。通过应用卡方检验，来评估该分区内变量间潜在关系的统计显著性。如果检验结果表明存在显著关系，则将该分区拆分为四个子分区，并将每个子分区重新压入栈中，以便进一步处理。这一过程持续进行，直至所有分区均满足停止准则或变为同质（均匀）。

if(! splittable && fullXstop-fullXstart > 30 && fullYstop-fullYstart > 30)
           {
            ipx = fullXstart - 1 ;
            ipy = fullYstart - 1 ;
            for(i=0 ; i<4 ; i++)
              {
               xcut[i] = (fullXstop - fullXstart + 1) * (i+1) / 4 + fullXstart - 1 ;
               xfrac[i] = (xcut[i] - ipx) / (fullXstop - fullXstart + 1.0) ;
               ipx = xcut[i] ;
               ycut[i] = (fullYstop - fullYstart + 1) * (i+1) / 4 + fullYstart - 1 ;
               yfrac[i] = (ycut[i] - ipy) / (fullYstop - fullYstart + 1.0) ;
               ipy = ycut[i] ;
              }
            for(ix=0 ; ix<4 ; ix++)
              {
               for(iy=0 ; iy<4 ; iy++)
                 {
                  expected[ix*4+iy] = xfrac[ix] * yfrac[iy] *
                                      (currentDataStop - currentDataStart + 1) ;
                  actual44[ix*4+iy] = 0 ;
                 }
              }
            for(i=currentDataStart ; i<=currentDataStop ; i++)
              {
               k = m_indices[i] ;
               for(ix=0 ; ix<3 ; ix++)
                 {
                  if(x[k] <= xcut[ix])
                     break ;
                 }
               for(iy=0 ; iy<3 ; iy++)
                 {
                  if(m_ranks[k] <= ycut[iy])
                     break ;
                 }
               ++actual44[ix*4+iy] ;
              }
            testval = 0.0 ;
            for(ix=0 ; ix<4 ; ix++)
              {
               for(iy=0 ; iy<4 ; iy++)
                 {
                  diff = fabs(actual44[ix*4+iy] - expected[ix*4+iy]) - 0.5 ;
                  testval += diff * diff / expected[ix*4+iy] ;
                 }
              }
            splittable = (testval > 3 * m_chi_crit)  ?  1 : 0 ;
           }
        }

      if(splittable)
        {

         for(i=currentDataStart ; i<=currentDataStop ; i++)
            current_indices[i] = m_indices[i] ;

         ipos = currentDataStart ;
         for(iSubRec=0 ; iSubRec<4 ; iSubRec++)
           {

            if(actual[iSubRec] >= 3)
              {
               stack[nstack].Xstart = trialXstart[iSubRec] ;
               stack[nstack].Xstop = trialXstop[iSubRec] ;
               stack[nstack].Ystart = trialYstart[iSubRec] ;
               stack[nstack].Ystop = trialYstop[iSubRec] ;
               stack[nstack].DataStart = ipos ;
               stack[nstack].DataStop = ipos + actual[iSubRec] - 1 ;
               ++nstack ;
               if(ArrayResize(stack,nstack+1,256)<0)
                 {
                   Print(__FUNCTION__," arrayresize error ", GetLastError());
                   return EMPTY_VALUE;
                 }

               if(iSubRec == 0)
                 {
                  for(i=currentDataStart ; i<=currentDataStop ; i++)
                    {
                     k = current_indices[i] ;
                     if(x[k] <= centerX  &&  m_ranks[k] <= centerY)
                        m_indices[ipos++] = current_indices[i] ;
                    }
                 }
               else
                  if(iSubRec == 1)
                    {
                     for(i=currentDataStart ; i<=currentDataStop ; i++)
                       {
                        k = current_indices[i] ;
                        if(x[k] <= centerX  &&  m_ranks[k] > centerY)
                           m_indices[ipos++] = current_indices[i] ;
                       }
                    }
                  else
                     if(iSubRec == 2)
                       {
                        for(i=currentDataStart ; i<=currentDataStop ; i++)
                          {
                           k = current_indices[i] ;
                           if(x[k] > centerX  &&  m_ranks[k] <= centerY)
                              m_indices[ipos++] = current_indices[i] ;
                          }
                       }
                     else
                       {
                        for(i=currentDataStart ; i<=currentDataStop ; i++)
                          {
                           k = current_indices[i] ;
                           if(x[k] > centerX  &&  m_ranks[k] > centerY)
                              m_indices[ipos++] = current_indices[i] ;
                          }
                       }
              }

            else
              {
               if(actual[iSubRec] > 0)
                 {
                  px = (trialXstop[iSubRec] - trialXstart[iSubRec] + 1.0) / m_size ;
                  py = (trialYstop[iSubRec] - trialYstart[iSubRec] + 1.0) / m_size ;
                  pxy = (double) actual[iSubRec] / m_size ;
                  MI += pxy * log(pxy / (px * py)) ;
                 }
              }
           }
        }
      else
        {
         px = (fullXstop - fullXstart + 1.0) / m_size ;
         py = (fullYstop - fullYstart + 1.0) / m_size ;
         pxy = (currentDataStop - currentDataStart + 1.0) / m_size ;
         MI += pxy * log(pxy / (px * py)) ;
        }
     }

   return MI;

对于因卡方检验不显著或数据点数量过少而被判定为同质的分区，会计算其互信息。该计算过程涉及估计分区内的联合概率分布和边缘概率分布，并应用以下公式。

互信息贡献量

该过程会递归进行，持续对分区进行拆分和评估，直至无法再进行显著拆分为止。最终，通过汇总所有终端分区的贡献量得出互信息的估计值，从而基于数据结构提供全面的估计。若该方法遇到错误，则返回与MQL5内置常量EMPTY_VALUE等效的值。在接下来的章节中，我们将使用Capm类来实现渐进特征选择算法。


基于互信息的预测变量选择

当将互信息应用于特征选择任务时，我们的目标是从大量候选预测变量中识别出一个子集，使该子集与目标变量的联合依赖性最大化。这种联合依赖性是互信息的推广，其中一个量是变量集合，而非单个变量。由于维度灾难，随着预测变量子集规模的扩大，计算联合依赖性变得在计算上不可行。这是因为对数据进行所需的多维分区会导致单元格极度稀疏，尤其是当维度（预测变量）数量增加时。

此外，即使对于中等规模的特征集，可能的特征组合数量也可能非常庞大。这种组合爆炸现象使得穷举搜索所有可能的子集变得不切实际，因此需要采用高效的搜索策略。

一种常见的特征选择方法是前向逐步选择法。该方法从空模型开始，根据所选标准迭代地添加能使模型性能提升最大的特征。尽管这种方法效率较高，但其存在贪婪算法的局限性。它可能会错过最优的特征组合，例如，两个特征组合在一起时，其预测能力远优于单独使用任何一个特征的情况。这是因为该算法专注于每一步的增量改进，但可能会忽略特征之间的协同关系。

尽管还存在其他技术，如高阶方法和后向选择法，但由于前向逐步选择法具有计算效率高的优势，因此在大规模数据集和计算资源有限的情况下，它仍然是最实用的选择。此外，当应用于联合依赖性的特定场景时，前向逐步选择法展现出一种令人惊讶的特性：即使它没有明确优化联合依赖性度量，也能有效地逼近最优特征子集。

Peng、Long和Ding提出了一种基于互信息的特征选择算法，即最大相关性最小冗余（MRMR）准则。该算法无需显式计算联合依赖性（这在计算上不可行），即可有效地逼近最优特征子集。MRMR准则基于相关性和冗余性的概念。特征集S与目标变量Y的相关性定义为Y与S中每个特征的平均互信息。

相关性

其中，|S| 表示特征集 S 中的特征数量，I(Xi, Y) 表示特征Xi 与目标变量 Y 之间的互信息。

虽然最大化相关性是直观的，但由于冗余性的存在，它可能导致次优的特征子集。如果我们仅根据特征与目标变量的个体相关性来选择特征，最终可能会得到一组高度相关的特征，这些特征提供的额外信息很少。为了解决这个问题，我们不仅需要考虑特征的相关性，还需要考虑其与已选特征的冗余性。

候选预测变量 Xi 相对于已选特征集 S 的冗余性，定义为Xi与S中每个特征之间的平均互信息：

冗余性

其中，I(Xi, Xj) 是候选特征Xi与已选特征集S中的特征Xj之间的互信息。

值得注意的是，冗余性的数学表达式与相关性的表达式完全相同。唯一的区别在于解释方式。当计算候选特征与目标变量之间的互信息时，我们称之为相关性。然而，当计算候选特征与已选特征集中特征之间的互信息时，我们称之为冗余性。

本质上，这两个概念都涉及测量两个变量之间共享的信息量，但具体情境决定了它是被视为相关还是冗余。在相关性中，我们关注的是特征与目标变量之间的关系，而在冗余性中，我们关注的是候选特征与已选特征之间的重叠程度。

总结来说，MRMR算法的操作步骤如下：

  1. 选择与目标变量Y互信息最高的特征作为第一个特征。
  2. 在随后的每一步中，算法选择满足以下准则的特征：

MRMR准则

该准则平衡了特征与目标变量Y的相关性（即 I(Xi, Y)）及其与已选特征的冗余性。MRMR算法的显著之处在于，尽管直接优化联合依赖性在计算上不可行，但它仍能有效地逼近最优特征子集。这一结果虽然在直观上并不明显，但在原始论文中得到了数学证明。

为了评估所选特征的统计显著性，可以采用两种蒙特卡洛置换（MCP）检验：

  • 单个特征显著性检验：该检验通过将每个特征与目标变量的相关性，与通过置换得到的0分布进行比较，来评估每个特征的显著性。通过置换特征值，我们创建一个假设特征与目标变量无关的0分布。如果观察到的相关性显著高于置换后的值，则可以得出结论，该特征很可能具有真正的信息价值。
  • 整体模型显著性检验：该检验通过将所选特征集中各个特征相关性的累积和，与通过置换得到的0分布进行比较，来评估整个所选特征集的显著性。通过置换目标变量值，我们创建一个假设特征与目标变量无关的0分布。如果观察到的累积相关性显著高于置换后的值，则可以得出结论，所选特征集很可能能够预测目标变量。


基于互信息的渐进特征选择算法的MQL5实现

mutualinfo.mqh中的Cmrmr类实现了MRMR算法，其构造函数接受以下参数：
  • num_reps：蒙特卡洛置换检验的重复次数。其设置为小于或等于1的值将禁用这些检验。
  • max_preds：要选择的最大特征数量。
  • chisquare_thresh：用于估计互信息的自适应分区方法的卡方检验阈值。
  • m_verbose：一个布尔标志，用于控制是否在算法执行期间显示详细输出。
//+-----------------------------------------------------------------------+
//|Relevance minus redundancy for building an optimal subset of predictors|
//+-----------------------------------------------------------------------+
class Cmrmr
  {
private:
   int               m_max_preds;
   int               m_reps;
   bool              m_verbose;
   int               m_samples;
   int               m_vars;
   matrix            m_preds;
   vector            m_target;
   vector            m_crits;
   double            m_chisquare_thresh;
   ulong             m_selected_vars[];
   matrix            m_critical_values;
   vector            mutualinfo(int which_size,int &which[],vector &targs);
public:
                     Cmrmr(int num_reps,int max_preds, double chisquare_thresh,bool verbose);
                    ~Cmrmr(void);
   bool              StepWise(matrix &predictors, vector &targets);
   matrix            GetCriticalValues(void);
   bool              GetSelectedVars(ulong &output[]);
  };

要使用Cmrmr类启动特征选择流程，我们需要实例化该类的一个对象，然后调用其StepWise()方法。StepWise()方法主要接受两个输入参数：

  • 候选预测变量矩阵：该矩阵应包含从业者希望考虑的所有潜在特征的数据。
  • 目标变量向量：该向量应包含目标变量的值。

该方法返回一个布尔值，用以指示特征选择过程是否成功。Cmrmr类中的StepWise()方法会初始化必要的数据结构，然后进入蒙特卡洛置换（MCP）检验循环。如果未要求进行MCP检验（即 num_reps 小于或等于 1），则该循环仅执行一次。在循环内部，该方法会使用Capm类的一个实例，计算每个候选预测变量与目标变量之间的互信息。

//+------------------------------------------------------------------+
//|  Stepwise feature selection based on mutual information          |
//+------------------------------------------------------------------+
bool Cmrmr::StepWise(matrix &predictors,vector &targets)
  {
   if(m_selected_vars.Size())
      ArrayFree(m_selected_vars);

   m_preds = predictors;
   m_samples = int(m_preds.Rows());
   m_vars = int(m_preds.Cols());
   m_max_preds = m_max_preds>=m_vars?m_vars:m_max_preds;
   m_target = targets;
   
   
   int i, j, k, ivar, irep;
   int index[], stepwise_mcpt_count[], solo_mcpt_count[], stepwise_ivar[], which_preds[],original_stepwise_ivar[] ;

   int nkept,best_ivar;
   vector casework(m_samples), sorted(m_samples), mutual(m_vars);
   vector crit(m_vars), relevance(m_vars), original_relevance(m_vars), current_crits(m_vars), sorted_crits(m_vars);
   double best_crit, dtemp, group_pvalue,solo_pvalue;
   vector stepwise_crit(m_vars), original_stepwise_crit(m_vars);
   double sum_relevance;
   vector original_sum_relevance(m_vars), sum_redundancy(m_vars);
   vector target = m_target;

   best_crit = -DBL_MAX;
   best_ivar = -1;
   nkept = m_max_preds;

   if(ArrayResize(index,m_vars)<0 || ArrayResize(stepwise_mcpt_count,m_vars)<0 ||
      ArrayResize(solo_mcpt_count,m_vars)<0 || ArrayResize(which_preds,m_vars)<0 || ArrayResize(stepwise_ivar,m_vars)<0 ||
      ArrayResize(original_stepwise_ivar,m_vars)<0)
     {
      Print(__FUNCTION__," array resize error ", GetLastError());
      return false;
     }

   int unif_error;
   for(irep=0 ; irep<m_reps ; irep++)
     {
      if(irep)
        {
         i = m_samples ;
         while(i > 1)
           {
            j = (int)(MathRandomUniform(0.0,1.0,unif_error) * i);
            if(unif_error)
              {
               Print(__FUNCTION__," Mathrandomuniform error ", GetLastError());
               return false;
              }
            if(j >= i)
               j = i - 1 ;
            dtemp = target[--i] ;
            target[i] = target[j] ;
            target[j] = dtemp ;
           }
        }

      for(i=0 ; i<m_vars ; i++)
         which_preds[i] = i ;

      crit = mutualinfo(m_vars,which_preds,target);

      for(ivar=0 ; ivar<m_vars ; ivar++)
        {
         relevance[ivar] = crit[ivar] ;
         if(ivar == 0  ||  crit[ivar] > best_crit)
           {
            best_crit = crit[ivar] ;
            best_ivar = ivar ;
           }
        }

      stepwise_crit[0] = best_crit ;
      stepwise_ivar[0] = best_ivar ;
      sum_relevance = best_crit ;
      if(irep == 0)
        {

         original_stepwise_crit[0] = best_crit ;
         original_stepwise_ivar[0] = best_ivar ;
         original_sum_relevance[0] = sum_relevance ;
         stepwise_mcpt_count[0] = 1 ;

         for(ivar=0 ; ivar<m_vars ; ivar++)
           {
            index[ivar] = ivar ;
            original_relevance[ivar] = sorted_crits[ivar] = current_crits[ivar] = crit[ivar] ;
            solo_mcpt_count[ivar] = 1 ;
           }

         if(!qsortdsi(0, m_vars-1, sorted_crits, index))
           {
            Print(__FUNCTION__, " failed qsort ");
            return false;
           }

         if(m_verbose)
            Print("       Variable         MI") ;

         for(i=m_vars-1 ; i>=0 ; i--)
           {
            k = index[i] ;

            if(m_verbose)
               PrintFormat("%15s %12.4lf",string(k), current_crits[k]) ;
           }
        }

      else
        {
         if(sum_relevance >= original_sum_relevance[0])
            ++stepwise_mcpt_count[0] ;
         for(ivar=0 ; ivar<m_vars ; ivar++)
           {
            if(relevance[ivar] >= original_relevance[ivar])
               ++solo_mcpt_count[ivar] ;
           }
        }

      for(i=0 ; i<m_vars ; i++)
         sum_redundancy[i] = 0.0 ;

      for(nkept=1 ; nkept<m_max_preds ; nkept++)
        {

         if(irep == 0 && m_verbose)
           {
            Print("Predictors so far   Relevance   Redundancy   Criterion") ;
            for(i=0 ; i<nkept ; i++)
              {
               k = stepwise_ivar[i] ;
               PrintFormat("%15s %12.4lf %12.4lf %12.4lf",string(k), relevance[k], relevance[k] - stepwise_crit[i], stepwise_crit[i]) ;
              }
           }

         k = 0 ;
         for(i=0 ; i<m_vars ; i++)
           {
            for(j=0 ; j<nkept ; j++)
              {
               if(stepwise_ivar[j] == i)
                  break ;
              }
            if(j == nkept)
               which_preds[k++] = i ;
           }
         if(k != (m_vars - nkept))
           {
            Print(__FUNCTION__, " failed assertion ", __LINE__);
            return false;
           }

         k = stepwise_ivar[nkept-1] ;

         casework = m_preds.Col(k);

         crit = mutualinfo(m_vars-nkept,which_preds,casework);

         for(i=0 ; i<(m_vars-nkept) ; i++)
           {
            k = which_preds[i] ;
            sum_redundancy[k] += crit[i] ;
            index[i] = k ;
            sorted_crits[i] = current_crits[i] = ((relevance[k] - sum_redundancy[k]) / double(nkept)) ;
            if(i == 0  ||  current_crits[i] > best_crit)
              {
               best_crit = current_crits[i] ;
               best_ivar = k ;
              }
           }

         stepwise_crit[nkept] = best_crit ;
         stepwise_ivar[nkept] = best_ivar ;
         sum_relevance += relevance[best_ivar] ;

         if(irep == 0)
           {
            original_stepwise_crit[nkept] = best_crit ;
            original_stepwise_ivar[nkept] = best_ivar ;
            original_sum_relevance[nkept] = sum_relevance ;
            stepwise_mcpt_count[nkept] = 1 ;

            if(!qsortdsi(0, m_vars-nkept-1, sorted_crits, index))
              {
               Print(__FUNCTION__, " failed qsort ");
               return false;
              }

            if(m_verbose)
              {
               Print("Additional candidates, in order of decreasing relevance minus redundancy") ;
               Print("       Variable     Relevance   Redundancy   Criterion") ;

               for(i=m_vars-nkept-1 ; i>=0 ; i--)
                 {
                  k = index[i] ;
                  PrintFormat("%15s %12.4lf %12.4lf %12.4lf",
                              string(k), relevance[k], sum_redundancy[k] / nkept,
                              relevance[k] - sum_redundancy[k] / nkept) ;
                 }
              }
           }

         else
           {
            if(sum_relevance >= original_sum_relevance[nkept])
               ++stepwise_mcpt_count[nkept] ;
           }

        }

     }
//---
   m_critical_values = matrix::Zeros(nkept,m_reps>1?5:3);
//---
   if(ArrayResize(m_selected_vars,nkept)<0)
     {
      Print(__FUNCTION__, " failed array resize ", GetLastError());
      return false;
     }
//---
   if(m_verbose)
     {
      if(m_reps > 1)
         Print("Final predictors ||  Relevance ||  Redundancy || Criterion  ||  Solo pval || Group pval") ;
      else
         Print("Final predictors  ||  Relevance ||  Redundancy ||  Criterion") ;
     }
//---
   for(i=0 ; i<nkept ; i++)
     {
      k = original_stepwise_ivar[i] ;
      m_selected_vars[i] = ulong(k);

      m_critical_values[i][0] = original_relevance[k];
      m_critical_values[i][1] = original_relevance[k] - original_stepwise_crit[i];
      m_critical_values[i][2] = original_stepwise_crit[i];

      if(m_critical_values.Cols()>3)
        {
         group_pvalue = (double) stepwise_mcpt_count[i] / (double) m_reps;
         solo_pvalue = (double) solo_mcpt_count[k] / (double) m_reps;
         m_critical_values[i][3] = solo_pvalue;
         m_critical_values[i][4] = group_pvalue;
        }

      if(m_verbose)
        {
         if(m_reps > 1)
            PrintFormat("%15s %12.4lf %12.4lf %12.4lf    %8.3lf    %8.3lf",string(k), m_critical_values[i][0], m_critical_values[i][1], m_critical_values[i][2],solo_pvalue,group_pvalue) ;
         else
            PrintFormat("%15s %12.4lf %12.4lf %12.4lf",string(k), m_critical_values[i][0], m_critical_values[i][1], m_critical_values[i][2]);
        }
     }

   return true;

  }

mutualinfo()方法是Cmrmr类中的一个私有方法，负责使用自适应分区方法估算给定预测变量与目标变量之间的互信息。该方法返回一个互信息近似值向量，这些近似值对应于作为列索引传递给mutualinfo()方法的选定预测变量。

//+------------------------------------------------------------------+
//| continuous mutual infomation                                     |
//+------------------------------------------------------------------+
vector Cmrmr::mutualinfo(int which_size,int &which[],vector &targs)
  {
   vector res = vector::Zeros(m_vars);

   res.Fill(-DBL_MAX);

   vector vars;

   Capm mia(true,targs,m_chisquare_thresh);

   for(int i = 0; i<which_size; i++)
     {
      vars = m_preds.Col(which[i]);
      res[i] = mia.fit(vars);
     }

   return res;

  }

StepWise()方法将计算得到的互信息值存储在relevance向量中，因为在后续迭代中需要这些值来计算冗余项。算法会选择互信息值最高的特征作为第一个特征，并将其信息存储在stepwise_crit和stepwise_ivar变量中。此外，初始化sum_relevance变量以跟踪所选特征的累积相关性，该变量后续将用于蒙特卡洛置换（MCP）检验。

当处于初始的、未进行置换的复制过程中，算法会存储第一个所选特征的相关性和准则的原始值。同时，还会初始化整体模型显著性检验和单个特征显著性检验的计数器。此外，算法会输出一个表格，展示每个候选特征与目标变量之间的互信息，以提供对初始特征排名的洞察。该表格有助于可视化初始选择过程，并作为显著性检验阶段与置换值进行比较的参考。

如果当前不处于初始的、未进行置换的复制过程，算法将继续执行两个MCP检验，旨在评估由MRMR（最大相关性、最小冗余性）算法所选特征集的统计显著性。为了高效管理随着新特征加入所选集合而不断演变的冗余性，算法使用sum_redundancy数组。该数组初始化为0。该数组用于存储每个候选特征相对于已选特征的累积冗余性。

开始时，由于没有特征被选中，因此该数组中的所有值都设置为0。当新特征被添加到所选集合中时，必须更新其与每个剩余候选特征的冗余性。特征相对于不断扩大的所选集合的冗余性定义为候选特征与已选特征之间的平均互信息。计算新特征与每个剩余候选特征之间的互信息。这样就可以量化新特征与每个已选特征共享的信息量。通过将新添加特征与现有特征之间的互信息值相加，更新每个剩余特征的冗余性值。

更新冗余性值后，算法为每个剩余候选特征计算MRMR（最大相关性、最小冗余性）准则。MRMR准则通过从每个特征的相关性中减去其平均冗余性来计算。一旦为所有剩余候选特征计算出MRMR准则，选择MRMR得分最高的特征。算法迭代重复此过程，每次添加新特征时更新sum_redundancy数组并重新计算MRMR准则。随着更多特征被选中，剩余特征与所选集合的冗余性会持续更新，确保特征选择过程中考虑特征间不断演变的关系。

如果处于初始的、未进行置换的复制过程，则保留这些量的原始值，以便在后续置换检验中进行比较。否则，增加“组”置换检验和“单独”置换检验的计数器。在完成所需的迭代次数后，算法结束时显示一个表格，总结所选特征、它们的相关性以及置换检验的p值。在成功调用StepWise()方法后，可以调用GetCriticalValues()来获取详细模式下显示的数值表格矩阵。同时，GetSelectedVars()可检索与类的m_max_preds属性相对应的所选变量（列）的索引。 

//+------------------------------------------------------------------+
//| get the matrix of decision variables                             |
//| matrix rows arranged in terms of descending relevance of each    |
//| candidate predictor. Most relevant predictors critical values are|
//| in the first row, etc.                                           |
//| Column 0 is the calculated relevance                             |
//| Column 1 is the calculated redundancy                            |
//| Column 2 is the calculated Stepwise Mutual information           |
//| Column 3 and 4 will exist only if MCP test is specified          |
//| Column 3 is the Solo probability value                           |
//| Column 4 are the Group probability values                        |
//+------------------------------------------------------------------+
matrix Cmrmr::GetCriticalValues(void)
{
  return m_critical_values;
}
//+--------------------------------------------------------------------------+
//|get the indices of the most relevant predictors as determined by algorithm|
//+--------------------------------------------------------------------------+
bool Cmrmr::GetSelectedVars(ulong &output[])
{
 return (ArrayCopy(output,m_selected_vars)>0);
}

在下一节中，我们将通过合成数据集和真实数据集，了解如何应用Cmrmr类。


相关性减去冗余性的示例

为了有效展示最大相关性最小冗余性（MRMR）算法的应用，我们首先将其应用于一个合成数据集。该数据集包含100个样本和10个潜在预测变量（特征）。目标变量Y是通过将矩阵的前四列相加构造的，因此目标变量与这些列的预测变量之间存在确定性关系。然而，第4至第7列不提供关于目标变量的任何真实信息，作为无关特征存在。第8列和第9列的变量分别是第{1,3}列和第{0,2}列变量的组合。 

//---
   MathSrand(125);

   matrix rdata(100,10);

   rdata.Random(0.0,1.0);

   vector dep = rdata.Col(0) + rdata.Col(1) + rdata.Col(2) + rdata.Col(3);

   vector sum02 = rdata.Col(0) + rdata.Col(2);

   vector sum13 = rdata.Col(1) + rdata.Col(3);

   if(!rdata.Col(sum13,8) || !rdata.Col(sum02,9))
     {
      Print(" Failed column insertion ", GetLastError());
      return;
     }

目标是利用最大相关性最小冗余性（MRMR）算法识别出这些列作为最具相关性的预测变量。该实验在RelevanceMinusRedundancy.mq5脚本中实现，该脚本允许自定义各种超参数，包括蒙特卡洛置换次数和要选择的最大特征数量。脚本配置为在详细模式下运行，以便在特征选择过程中提供详细的输出信息。

//inputs
input int NumReplications = 100;
input int MaxPredictors = 10;
input double ChiSquareThreshold = 6.0;

RelevanceMinusRedundancy.mq5脚本的初始输出是每个候选预测变量与目标变量之间的互信息值。这些值按互信息降序排列，便于轻松识别最具相关性的特征。

Variable         MI
      9       0.2675
      8       0.1696
      0       0.1662
      3       0.1336
      2       0.0823
      1       0.0645
      6       0.0000
      7       0.0000
      4       0.0000
      5       0.0000

不出所料，MRMR算法选择的第一个特征是与目标变量具有最高个体互信息的特征。即列索引为9的特征。

当前已选预测变量   相关性   冗余性  准则
            9       0.2675       0.0000       0.2675

第二个被选中的特征是列索引为8的变量。从分析结果来看，该变量具有较高的相关性，且与第一个被选中的变量冗余性极小。

当前已选预测变量   相关性   冗余性  准则
            9       0.2675       0.0000       0.2675
            8       0.1696       0.0000       0.1696

第三个被选中的特征展示了冗余性如何影响选择过程。尽管列0、1、2和3的特征与目标变量直接相关，但它们与已选特征（列8和9）之间的高度冗余性降低了它们的整体MRMR准则值。因此，算法选择了与现有集合冗余性较低、但直接相关性较低的其他特征。只有在选择了几个不相关的特征后，算法才会开始优先考虑列0、1、2和3的特征，因为此时它们与已选集合的冗余性降低了。

其他候选变量（按相关性减去冗余性降序排列）
变量 相关性 冗余性 准则
    5       0.0000       0.0000       0.0000
    4       0.0000       0.0000       0.0000
    7       0.0000       0.0000       0.0000
    6       0.0000       0.0000       0.0000
    1       0.0645       0.0600       0.0044
    0       0.1662       0.1351       0.0311
    2       0.0823       0.1426      -0.0603
    3       0.1336       0.1781      -0.0445

虽然在这个合成示例中，无关变量的相关性为0，但值得注意的是，在现实场景中，无关特征仍可能与目标变量存在一定程度的关系，从而导致冗余性不为0。在这种情况下，MRMR算法能够有效平衡相关性与冗余性，优先选择能提供独特信息的特征。

最终预测变量 || 相关性 || 冗余性 || 准则值 || 独立检验p值 || 组合检验p值
9                     0.2675        0.0000       0.2675       0.010       0.010
8                     0.1696        0.0000       0.1696       0.010       0.010
4                     0.0000        0.0000       0.0000       1.000       0.010
5                     0.0000        0.0000       0.0000       1.000       0.010
6                     0.0000        0.0000       0.0000       1.000       0.010
7                     0.0000        0.0000       0.0000       1.000       0.010
1                     0.0645        0.0738      -0.0093       0.010       0.010
0                     0.1662        0.1811      -0.0149       0.010       0.010
2                     0.0823        0.1077      -0.0254       0.010       0.010
3                     0.1336        0.1584      -0.0247       0.010       0.010

通过蒙特卡洛置换检验得到的p值，为所选特征的统计显著性提供了重要的参考。p值越低，表明该特征真正具有信息性而非偶然出现的证据越强。通常采用0.05的p值阈值来判断统计显著性。

接下来，我们将展示该算法在更贴近实际的数据集上的应用。StepWiseFeatureSelectionByMutualInformation.mq5脚本展示了MRMR算法在金融数据中的实际应用。该脚本基于这样一个前提：利用一组技术指标预测特定交易品种的未来收益率。脚本针对不同回溯周期计算了四种指标：资金流量指数（MFI）、移动平均线（MA）、空头力量指标（Bears Power）和多头力量指标（Bulls Power）。这些指标将作为特征选择过程的候选特征。通过应用MRMR算法，脚本旨在识别出能够有效预测未来收益率的、最具信息量的指标组合。这涉及选择与目标变量（未来收益率）相关，同时与其他已选特征冗余度最低的特征。

//+------------------------------------------------------------------+
//|                  StepWiseFeatureSelectionByMutualInformation.mq5 |
//|                                  Copyright 2024, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2024, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#property version   "1.00"
#property script_show_inputs
#resource "\\Indicators\\LogReturns.ex5"
#include<mutualinfo.mqh>
#include<ErrorDescription.mqh>
//+------------------------------------------------------------------+
//|indicator type                                                    |
//+------------------------------------------------------------------+
enum SELECT_INDICATOR
  {
   MFI=0,//MFI
   MA,//MA
   BEARS,//BEARS
   BULLS//BULLS
  };
//--- input parameters
input int NumReplications = 100;
input int MaxPredictors = 10;
input double ChiSquareThreshold = 6.0;
input bool VerboseOutPut = false;
input uint     period_inc=2;//lookback increment
input uint     max_lookback=6;
input ENUM_MA_METHOD         AppliedMA = MODE_SMA;
input ENUM_APPLIED_PRICE     AppliedPrice = PRICE_CLOSE;
input datetime SampleStartDate=D'2019.12.31';
input datetime SampleStopDate=D'2022.12.31';
input string   SetSymbol="BTCUSD";
input ENUM_TIMEFRAMES SetTF = PERIOD_D1;
//----
string predictor_names[];                      // variable names
int size_sample,                      //training set size
    size_observations,                //size of of both training and testing sets combined
    maxperiod,                        //maximum lookback
    indicator_handle=INVALID_HANDLE;  //long moving average indicator handle
//---
vector indicator[];                   //indicator indicator values;
//---
matrix feature_matrix;          //full matrix of features;
//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+
void OnStart()
  {
//---get relative shift of sample set
   int samplestart,samplestop,num_features;
   samplestart=iBarShift(SetSymbol!=""?SetSymbol:NULL,SetTF,SampleStartDate);
   samplestop=iBarShift(SetSymbol!=""?SetSymbol:NULL,SetTF,SampleStopDate);
   num_features = int((max_lookback/period_inc)*4);
//---check for errors from ibarshift calls
   if(samplestart<0 || samplestop<0)
     {
      Print(ErrorDescription(GetLastError()));
      return;
     }
//---set the size of the sample sets
   size_observations=(samplestart - samplestop) + 1 ;
   maxperiod=int(max_lookback);
//---check for input errors
   if(size_observations<=0 || maxperiod<=0)
     {
      Print("Invalid inputs ");
      return;
     }
//---allocate memory
   if(ArrayResize(indicator,num_features+1)<0)
     {
      Print(ErrorDescription(GetLastError()));
      return;
     }
//----get the full collection of indicator values
   int period_len;
   int k=0;
//---
   for(SELECT_INDICATOR select_indicator = 0; select_indicator<4; select_indicator++)
     {
      for(int iperiod=0; iperiod<int((indicator.Size()-1)/4); iperiod++)
        {
         period_len=int((iperiod+1) * period_inc);
         int try
               =10;
         predictor_names.Push(EnumToString(select_indicator)+"_"+string(period_len));
         while(try)
           {
            switch(select_indicator)
              {
               case MFI:
                  indicator_handle=iMFI(SetSymbol!=""?SetSymbol:NULL,SetTF,period_len,VOLUME_TICK);
                  break;
               case MA:
                  indicator_handle=iMA(SetSymbol!=""?SetSymbol:NULL,SetTF,period_len,0,AppliedMA,AppliedPrice);
                  break;
               case BEARS:
                  indicator_handle=iBearsPower(SetSymbol!=""?SetSymbol:NULL,SetTF,period_len);
                  break;
               case BULLS:
                  indicator_handle=iBullsPower(SetSymbol!=""?SetSymbol:NULL,SetTF,period_len);
                  break;
              }

            if(indicator_handle==INVALID_HANDLE)
               try
                  --;
            else
               break;
           }

         if(indicator_handle==INVALID_HANDLE)
           {
            Print("Invalid indicator handle ",EnumToString(select_indicator)," ", GetLastError());
            return;
           }

         Comment("copying data to buffer for indicator ",period_len);
         try
               = 0;
         while(!indicator[k].CopyIndicatorBuffer(indicator_handle,0,samplestop,size_observations) && try
                  <10)
              {
               try
                  ++;
               Sleep(5000);
              }

         if(try
               <10)
               ++k;
         else
           {
            Print("error copying to indicator buffers ",GetLastError());
            Comment("");
            return;
           }

         if(indicator_handle!=INVALID_HANDLE && IndicatorRelease(indicator_handle))
            indicator_handle=INVALID_HANDLE;
        }
     }

//---resize matrix
   if(!feature_matrix.Resize(size_observations,indicator.Size()-1))
     {
      Print(__LINE__);
      Print(ErrorDescription(GetLastError()));
      Comment("");
      return;
     }
//---copy collected data to matrix
   for(ulong i = 0; i<feature_matrix.Cols(); i++)
      if(!feature_matrix.Col(indicator[i],i))
        {
         Print(__LINE__);
         Print(ErrorDescription(GetLastError()));
         Comment("");
         return;
        }
//---
   int try
         = 10;
   while(try
            >-1 && indicator_handle == INVALID_HANDLE)
        {
         indicator_handle=iCustom(SetSymbol!=""?SetSymbol:NULL,SetTF,"\\Indicators\\LogReturns",0,1,1);
         try
            --;
        }
//---
   if(try
         <0)
        {
         Print("Could not initialize returns indicator ");
         Print(ErrorDescription(GetLastError()));
         Comment("");
         return;
        }
   else
     {
      try
            = 10;
     }
//---
   while(try
            >-1 && !indicator[indicator.Size()-1].CopyIndicatorBuffer(indicator_handle,0,samplestop-1,size_observations))
        {
         Sleep(2000);
         try
            --;
        }
//---
   if(try
         <0)
        {
         Print("Could not collect returns indicator info ");
         Print(ErrorDescription(GetLastError()));
         Comment("");
         return;
        }
   else
     {
      IndicatorRelease(indicator_handle);
      Comment("");
     }
//--- display the dataset
   string console;
   for(uint i = 0; i<predictor_names.Size(); i++)
      console+=StringFormat(" %12s ",predictor_names[i]);
   console+=" NextBar Returns (Target) ";
   Print(console);
   for(ulong i = 0; i<feature_matrix.Rows(); i++)
     {
      console = "";
      for(ulong j = 0; j<feature_matrix.Cols(); j++)
         console += StringFormat(" %12.6lf ",feature_matrix[i][j]);
      console+=StringFormat(" %12.6lf ",indicator[indicator.Size()-1][i]);
      Print(console);
     }
//---
   Cmrmr mstep(NumReplications,MaxPredictors,ChiSquareThreshold,VerboseOutPut);
//---
   if(!mstep.StepWise(feature_matrix,indicator[indicator.Size()-1]))
      return;
//---
   Print(" Final predictor labels ");
   ulong variables[];
   if(mstep.GetSelectedVars(variables))
     {
      for(uint i = 0; i<variables.Size(); i++)
         Print(predictor_names[variables[i]]);
     }
   return;
  }
//+------------------------------------------------------------------+
以下是从MetaTrader 5终端收集的指标数据集片段，最后一列为目标变量。该数据集共包含12个候选预测变量（即技术指标）。 
MFI_2         MFI_4         MFI_6          MA_2          MA_4          MA_6       BEARS_2       BEARS_4       BEARS_6       BULLS_2       BULLS_4       BULLS_6  NextBar Returns(Target)
0.000000     53.151442     46.921608   7213.654000   7279.552750   7260.007333    -76.371267   -101.867797   -107.113146     87.265733     61.769203     56.523854     -0.003564
0.000000     26.316595     34.451328   7180.962000   7244.558000   7255.420333    -41.795089    -70.889078    -83.462247     42.344911     13.250922      0.677753     -0.032447
0.000000      0.000000     36.720977   7053.740000   7133.697000   7204.281833   -127.731696   -210.813447   -251.244462    153.948304     70.866553     30.435538      0.054562

以下是使用脚本默认参数进行预测变量分析后得到的结果。首先，我们来看互信息值表格。从中可以看到，部分指标与目标变量之间的互信息为0。其中，周期长度为6的移动平均线指标与目标变量的互信息值最高。 
PS          Variable         MI
ND                 5       0.0308
LG                 4       0.0293
MJ                 3       0.0279
GM                 6       0.0227
JP                 9       0.0182
IS                 1       0.0165
MF                 8       0.0135
JI                 7       0.0126
HL                 0       0.0000
JO                 2       0.0000
GS                10       0.0000
HD                11       0.0000

显而易见，MA_6成为算法选中的第一个特征。然已经做出了首个选择，接下来我们不妨看看是否存在其他候选预测变量与MA_6存在信息共享（即存在冗余性）。这一点可通过观察下方的冗余值来推断。 

ME   当前预测变量 相关性 冗余性 准则值
GH                 5       0.0308       0.0000       0.0308
II   其他候选变量（按相关性减去冗余性降序排列）
FE          变量   相关性   冗余性   准则
LE                 0       0.0000       0.0095      -0.0095
ED                 1       0.0165       0.0869      -0.0704
FD                 2       0.0000       0.1149      -0.1149
PE                11       0.0000       0.2768      -0.2768
ID                 8       0.0135       0.3251      -0.3116
NS                 7       0.0126       0.3492      -0.3366
CR                10       0.0000       0.3484      -0.3484
ES                 6       0.0227       0.4030      -0.3802
IS                 9       0.0182       0.4285      -0.4103
CR                 3       0.0279       2.5363      -2.5084
DR                 4       0.0293       3.0096      -2.9803
MFI_2 将成为算法的下一个选择，因为在所有剩余候选预测变量中，它的冗余性最低，尽管它与目标变量的相关性为0。此外，请注意，列索引为3和4的预测变量（分别为MA_2和MA_4）与MA_6的冗余性最高。这是很合理的，因为它们是使用更短窗口长度计算的同一指标。以下是再进行两次选择后的预测变量集，这两次选择的变量与目标变量的相关性也均为0。值得注意的是，在未被选中的候选预测变量集中，排名最靠前的预测变量。由于最后两次选择对已选预测变量集整体产生了稀释效应，算法开始考虑那些与目标变量相关的预测变量。 
PQ   当前预测变量 相关性 冗余性 准则值
FL                 5       0.0308       0.0000       0.0308
JK                 0       0.0000       0.0095      -0.0095
DK                 2       0.0000       0.1796      -0.1796
EE  其他候选变量（按相关性减去冗余性降序排列）
JH         变量   相关性   冗余性   准则
CH                6       0.0227       0.1469      -0.1241
PH                9       0.0182       0.1465      -0.1283
GI               10       0.0000       0.2024      -0.2024
PG                7       0.0126       0.2133      -0.2007
PF               11       0.0000       0.2308      -0.2308
GG                8       0.0135       0.2664      -0.2529
QG                1       0.0165       0.4679      -0.4514
KF                3       0.0279       0.8840      -0.8561
LF                4       0.0293       1.0372      -1.0079
我们跳过中间步骤，直接展示算法最终选定的预测变量集，以此结束对分析结果的讨论。 
OQ  最终选定的预测变量 || 相关性 || 冗余性 || 准则值 || 单变量p值 || 组变量p值
HN                5       0.0308       0.0000       0.0308       0.010       0.010
DJ                0       0.0000       0.0095      -0.0095       1.000       0.010
JG                2       0.0000       0.1796      -0.1796       1.000       0.010
HD                6       0.0227       0.1620      -0.1393       0.010       0.010
FQ                9       0.0182       0.2023      -0.1842       0.010       0.010
DL               11       0.0000       0.2798      -0.2798       1.000       0.010
KJ                1       0.0165       0.2932      -0.2767       0.010       0.010
OG                7       0.0126       0.3298      -0.3172       0.010       0.010
OE               10       0.0000       0.4151      -0.4151       1.000       0.010
JP                8       0.0135       0.4585      -0.4450       0.010       0.010
RN   最终预测变量标签
NO  MA_6
DE  MFI_2
NN  MFI_6
KP  BEARS_2
DJ  BULLS_2
FL  BULLS_6
PQ  MFI_4
MK  BEARS_4
FM  BULLS_4
OF  BEARS_6


结论

在本文中，我们探讨了互信息在特征选择中的应用，重点聚焦于彭汉川（Peng）、龙乐（Long）和丁晓青（Ding）提出的MRMR（最大依赖性、最大相关性、最小冗余性）算法。我们首先讨论了连续变量互信息的估计方法，强调了其在捕捉特征与目标变量之间复杂非线性依赖关系方面的能力。通过对MRMR算法的深入探究，我们展示了该算法如何有效平衡两大目标：一方面最大化与目标变量的相关性，另一方面最小化与已选特征的冗余性。

该算法通过基于MRMR分数迭代添加特征，并通过蒙特卡洛置换检验评估其统计显著性，确保所选特征集能够为模型提供最具信息量和独特性的预测变量。合成数据和真实数据示例均表明，MRMR算法在实际数据场景中具有实用价值，尤其是在无关特征和多重共线性问题常常使特征选择变得复杂的情况下。MRMR算法能够优先选择相关特征并避免冗余特征，这使其在特征选择过程中始终保持重要价值，特别是在变量间关系错综复杂且非线性的情况下。

文中提及的所有代码均附于下方。下表对本文附上的所有源代码文件进行了说明。

文件名
说明
MQL5/include/np.mqh
 一个包含矩阵与向量实用函数的头文件。
MQL5/include/mutualinfo.mqh
 该头文件包含Capm类的定义，该类实现了用于估计连续变量互信息的自适应分区方法。同时还包含Cmrmr类的定义，该类实现了基于互信息的渐进特征选择算法。
MQL5/scripts/RelevanceMinusRedundancy.mq5
 一个使用合成数据演示Cmrmr类用法的脚本。
MQL5/scripts/StepWiseFeatureSelectionByMutualInformation.mq5
 另一个脚本，同样演示了Cmrmr类在更贴近实际的数据集上的应用。

本文由MetaQuotes Ltd译自英文
原文地址： https://www.mql5.com/en/articles/16416

附加的文件 |
下载ZIP
np.mqh (74.16 KB)
mutualinfo.mqh (26.74 KB)
RelevanceMinusRedundancy.mq5 (1.69 KB)
StepWiseFeatureSelectionByMutualInformation.mq5 (7.39 KB)
Mql5.zip (19.35 KB)

注意: MetaQuotes Ltd.将保留所有关于这些材料的权利。全部或部分复制或者转载这些材料将被禁止。

本文由网站的一位用户撰写，反映了他们的个人观点。MetaQuotes Ltd 不对所提供信息的准确性负责，也不对因使用所述解决方案、策略或建议而产生的任何后果负责。

该作者的其他文章

前往讨论
您应当知道的 MQL5 向导技术（第 43 部分）：依据 SARSA 进行强化学习 您应当知道的 MQL5 向导技术（第 43 部分）：依据 SARSA 进行强化学习
SARSA 是 “State-Action-Reward-State-Action” 的缩写，是另一种能在实现强化学习时运用的算法。故此，正如我们在 Q-学习 和 DQN 中看到的那样，我们考察了如何在向导汇编的智能系统中探索和实现它，将其作为独立模型，而不仅仅是一种训练机制。
开发回放系统（第 68 部分）：取得正确的时间（一） 开发回放系统（第 68 部分）：取得正确的时间（一）
今天，我们将继续努力，让鼠标指针告诉我们在流动性较低期间，一根柱形上还剩下多少时间。尽管乍一看似乎很简单，但实际上这项任务要困难得多。这涉及一些我们必须克服的障碍。因此，为了理解以下部分，您必须很好地理解子系列第一部分的材料。
MQL5 交易管理面板开发指南（第六部分）：交易管理面板（续篇） MQL5 交易管理面板开发指南（第六部分）：交易管理面板（续篇）
在本文中，我们对多功能管理面板的“交易面板”进行升级。我们引入一个强大的辅助函数，大幅简化代码，提高可读性、可维护性与运行效率。同时演示如何无缝集成更多按钮，并优化界面，以支持更广泛的交易任务。无论是持仓管理、订单调整，还是简化交互，本文将助您打造稳健且易用的交易管理面板。
数据科学和机器学习（第 31 部分）：利用 CatBoost AI 模型进行交易 数据科学和机器学习（第 31 部分）：利用 CatBoost AI 模型进行交易
CatBoost AI 模型最近在机器学习社区中广受欢迎，因为它们的预测准确性、效率、及针对分散和困难数据集的健壮性。在本文中，我们将详细讨论如何实现这些类型的模型，进而尝试进击外汇市场。