卡方选择器

  特征选择试图识别相关的特征用于模型构建。它改变特征空间的大小，它可以提高速度以及统计学习行为。ChiSqSelector实现卡方特征选择，它操作于带有类别特征的标注数据。
ChiSqSelector根据独立的卡方测试对特征进行排序，然后选择排序最高的特征。下面是一个使用的例子。

  1. import org.apache.spark.SparkContext._
  2. import org.apache.spark.mllib.linalg.Vectors
  3. import org.apache.spark.mllib.regression.LabeledPoint
  4. import org.apache.spark.mllib.util.MLUtils
  5. import org.apache.spark.mllib.feature.ChiSqSelector
  6. // 加载数据
  7. val data = MLUtils.loadLibSVMFile(sc, "data/mllib/sample_libsvm_data.txt")
  8. // 卡方分布需要类别特征，所以对特征除一个整数。虽然特征是double类型，
  9. //但是ChiSqSelector将每个唯一的值当做一个类别
  10. val discretizedData = data.map { lp =>
  11.   LabeledPoint(lp.label, Vectors.dense(lp.features.toArray.map { x => (x / 16).floor } ) )
  12. }
  13. // Create ChiSqSelector that will select top 50 of 692 features
  14. val selector = new ChiSqSelector(50)
  15. // Create ChiSqSelector model (selecting features)
  16. val transformer = selector.fit(discretizedData)
  17. // Filter the top 50 features from each feature vector
  18. val filteredData = discretizedData.map { lp => 
  19.   LabeledPoint(lp.label, transformer.transform(lp.features)) 
  20. }

  下面看看选择特征的实现，入口函数是fit。

  1.  def fit(data: RDD[LabeledPoint]): ChiSqSelectorModel = {
  2.     //计算数据卡方值
  3.     val indices = Statistics.chiSqTest(data)
  4.       .zipWithIndex.sortBy { case (res, _) => -res.statistic }
  5.       .take(numTopFeatures)
  6.       .map { case (_, indices) => indices }
  7.       .sorted
  8.     new ChiSqSelectorModel(indices)
  9.   }

  这里通过Statistics.chiSqTest计算卡方检测的值。下面需要了解卡方检测的理论基础。

1 卡方检测

1.1 什么是卡方检测

  卡方检验是一种用途很广的计数资料的假设检验方法。它属于非参数检验的范畴，主要是比较两个及两个以上样本率( 构成比）以及两个分类变量的关联性分析。
其根本思想就是在于比较理论频数和实际频数的吻合程度或拟合优度问题。

1.2 卡方检测的基本思想

  卡方检验是以${X}^{2}$分布为基础的一种常用假设检验方法，它的无效假设H0是：观察频数与期望频数没有差别。

  该检验的基本思想是：首先假设H0成立，基于此前提计算出${X}^{2}$值，它表示观察值与理论值之间的偏离程度。根据${X}^{2}$分布及自由度可以确定在H0假设成立的情况下获得当前统计量及更极端情况的概率P。
如果P值很小，说明观察值与理论值偏离程度太大，应当拒绝无效假设，表示比较资料之间有显著差异；否则就不能拒绝无效假设，尚不能认为样本所代表的实际情况和理论假设有差别。

1.3 卡方值的计算与意义

  卡方值表示观察值与理论值之问的偏离程度。计算这种偏离程度的基本思路如下。

  • 设A代表某个类别的观察频数，E代表基于H0计算出的期望频数，A与E之差称为残差。

  • 残差可以表示某一个类别观察值和理论值的偏离程度，但如果将残差简单相加以表示各类别观察频数与期望频数的差别，则有一定的不足之处。
    因为残差有正有负，相加后会彼此抵消，总和仍然为0，为此可以将残差平方后求和。

  • 另一方面，残差大小是一个相对的概念，相对于期望频数为10时，期望频数为20的残差非常大，但相对于期望频数为1000时20的残差就很小了。
    考虑到这一点，人们又将残差平方除以期望频数再求和，以估计观察频数与期望频数的差别。

  进行上述操作之后，就得到了常用的${X}^{2}$统计量。其计算公式是：

5.1

  当n比较大时，卡方统计量近似服从k-1(计算E_i时用到的参数个数)个自由度的卡方分布。由卡方的计算公式可知，当观察频数与期望频数完全一致时，卡方值为0；观察频数与期望频数越接近，两者之间的差异越小，卡方值越小；
反之，观察频数与期望频数差别越大，两者之间的差异越大，卡方值越大。

2 卡方检测的源码实现

  在MLlib中，使用chiSquaredFeatures方法实现卡方检测。它为每个特征进行皮尔森独立检测。下面看它的代码实现。

  1. def chiSquaredFeatures(data: RDD[LabeledPoint],
  2.       methodName: String = PEARSON.name): Array[ChiSqTestResult] = {
  3.     val maxCategories = 10000
  4.     val numCols = data.first().features.size
  5.     val results = new Array[ChiSqTestResult](numCols)
  6.     var labels: Map[Double, Int] = null
  7.     // 某个时刻至少1000列
  8.     val batchSize = 1000
  9.     var batch = 0
  10.     while (batch * batchSize < numCols) {
  11.       val startCol = batch * batchSize
  12.       val endCol = startCol + math.min(batchSize, numCols - startCol)
  13.       val pairCounts = data.mapPartitions { iter =>
  14.         val distinctLabels = mutable.HashSet.empty[Double]
  15.         val allDistinctFeatures: Map[Int, mutable.HashSet[Double]] =
  16.           Map((startCol until endCol).map(col => (col, mutable.HashSet.empty[Double])): _*)
  17.         var i = 1
  18.         iter.flatMap { case LabeledPoint(label, features) =>
  19.           if (i % 1000 == 0) {
  20.             if (distinctLabels.size > maxCategories) {
  21.               throw new SparkException
  22.             }
  23.             allDistinctFeatures.foreach { case (col, distinctFeatures) =>
  24.               if (distinctFeatures.size > maxCategories) {
  25.                 throw new SparkException
  26.               }
  27.             }
  28.           }
  29.           i += 1
  30.           distinctLabels += label
  31.           features.toArray.view.zipWithIndex.slice(startCol, endCol).map { case (feature, col) =>
  32.             allDistinctFeatures(col) += feature
  33.             (col, feature, label)
  34.           }
  35.         }
  36.       }.countByValue()
  37.       if (labels == null) {
  38.         // Do this only once for the first column since labels are invariant across features.
  39.         labels =
  40.           pairCounts.keys.filter(_._1 == startCol).map(_._3).toArray.distinct.zipWithIndex.toMap
  41.       }
  42.       val numLabels = labels.size
  43.       pairCounts.keys.groupBy(_._1).map { case (col, keys) =>
  44.         val features = keys.map(_._2).toArray.distinct.zipWithIndex.toMap
  45.         val numRows = features.size
  46.         val contingency = new BDM(numRows, numLabels, new Array[Double](numRows * numLabels))
  47.         keys.foreach { case (_, feature, label) =>
  48.           val i = features(feature)
  49.           val j = labels(label)
  50.           //带有标签的特征的出现次数
  51.           contingency(i, j) += pairCounts((col, feature, label))
  52.         }
  53.         results(col) = chiSquaredMatrix(Matrices.fromBreeze(contingency), methodName)
  54.       }
  55.       batch += 1
  56.     }
  57.     results
  58.   }

  上述代码主要对数据进行处理，获取带有标签的特征的出现次数，并用这个次数计算卡方值。真正获取卡方值的函数是chiSquaredMatrix。

  1.  def chiSquaredMatrix(counts: Matrix, methodName: String = PEARSON.name): ChiSqTestResult = {
  2.     val method = methodFromString(methodName)
  3.     val numRows = counts.numRows
  4.     val numCols = counts.numCols
  5.     // get row and column sums
  6.     val colSums = new Array[Double](numCols)
  7.     val rowSums = new Array[Double](numRows)
  8.     val colMajorArr = counts.toArray
  9.     val colMajorArrLen = colMajorArr.length
  10.     var i = 0
  11.     while (i < colMajorArrLen) {
  12.       val elem = colMajorArr(i)
  13.       if (elem < 0.0) {
  14.         throw new IllegalArgumentException("Contingency table cannot contain negative entries.")
  15.       }
  16.       //每列的总数
  17.       colSums(i / numRows) += elem
  18.       //每行的总数
  19.       rowSums(i % numRows) += elem
  20.       i += 1
  21.     }
  22.     //所有元素的总和
  23.     val total = colSums.sum
  24.     // second pass to collect statistic
  25.     var statistic = 0.0
  26.     var j = 0
  27.     while (j < colMajorArrLen) {
  28.       val col = j / numRows
  29.       val colSum = colSums(col)
  30.       if (colSum == 0.0) {
  31.         throw new IllegalArgumentException("Chi-squared statistic undefined for input matrix due to"
  32.           + s"0 sum in column [$col].")
  33.       }
  34.       val row = j % numRows
  35.       val rowSum = rowSums(row)
  36.       if (rowSum == 0.0) {
  37.         throw new IllegalArgumentException("Chi-squared statistic undefined for input matrix due to"
  38.           + s"0 sum in row [$row].")
  39.       }
  40.       //期望值
  41.       val expected = colSum * rowSum / total
  42.       //PEARSON
  43.       statistic += method.chiSqFunc(colMajorArr(j), expected)
  44.       j += 1
  45.     }
  46.     //自由度
  47.     val df = (numCols - 1) * (numRows - 1)
  48.     if (df == 0) {
  49.       // 1 column or 1 row. Constant distribution is independent of anything.
  50.       // pValue = 1.0 and statistic = 0.0 in this case.
  51.       new ChiSqTestResult(1.0, 0, 0.0, methodName, NullHypothesis.independence.toString)
  52.     } else {
  53.       //计算累积概率
  54.       val pValue = 1.0 - new ChiSquaredDistribution(df).cumulativeProbability(statistic)
  55.       new ChiSqTestResult(pValue, df, statistic, methodName, NullHypothesis.independence.toString)
  56.     }
  57.   }
  58.   //上述代码中的method.chiSqFunc(colMajorArr(j), expected)，调用下面的代码
  59.   val PEARSON = new Method("pearson", (observed: Double, expected: Double) => {
  60.       val dev = observed - expected
  61.       dev * dev / expected
  62.     })

  上述代码的实现和参考文献【2】中Test of independence的描述一致。

参考文献

【1】卡方检验

【2】Pearson’s chi-squared test