Alink漫談(二十) ：卡方檢驗原始碼解析
0x00 摘要
Alink 是阿里巴巴基於實時計算引擎 Flink 研發的新一代機器學習演算法平臺，是業界首個同時支援批式演算法、流式演算法的機器學習平臺。本文將帶領大家來分析 Alink 中 卡方檢驗 的實現。
因為Alink的公開資料太少，所以以下均為自行揣測，肯定會有疏漏錯誤，希望大家指出，我會隨時更新。
0x01 背景概念
問題：在南方的小明要去北京讀書，小明的父母都擔心起來，北方的生活能習慣嗎？是不是隻能吃麵？他們的腦子裡面提出了很多很多的假設，是不是都要驗證一下呢。
1.1 假設檢驗
統計假設是指我們對總體的猜測或判斷，比如中國廣為流行的地域劃分，北方人是不是都吃麵食？江浙是不是都喜歡吃甜食？江南的妹子是不是脾氣都很好？北方的男生是不是都大條？
為了證明，或者說絕對確定某個假設是正確或錯誤的，我們需要絕對知識，也就是說我們需要檢查所有的樣本，比如問問所有的北方人是不是都吃麵食，或者說問問江浙一帶的女生是不是都不吵架？
但是我們無法把所有樣本都統計一遍。
假設檢驗就是為了解決這個問題而誕生的，它使用隨機的樣本來判斷假設是否有理。
這來源於一個最基本的思路，那就是統計是從有限樣本推斷總體。既然我們不可能知道所有樣本，換句話說我們不可能知道總體到底是什麼統計性質，那麼我們只能拿有限的樣本做文章。
假設檢驗依據的原理是小概率事件原理。小概率事件是一個事件的發生概率，由於概率小，它在一次試驗中是幾乎不可能發生的。例如，拿我們生活的經驗來說，只買一次彩票就中大獎的機率是很小的，所以還沒看到過偶爾買彩票就中大獎的報導。
1.2 H0和H1是什麼？
小明父母的種種假設到底對不對呢？對是一種結果，不對是另一種結果。在假設檢驗中，我們也有兩部分，一部分叫做H0，一部分是H1。
H是英文單詞hypothesis的第一個字母，H0代表了原假設（null hypothesis），H1是備擇假設（alternative hypothesis）。換句話說，我們想在兩個假設中選一個，儘管這個“選”不是通常意義下的“二選一”。
- H0 原假設又叫零假設，一般來說，我們把認為想收集證據反對的假設稱為0假設，比如太陽繞著地球轉，鳥是不會飛的。小明這裡的 H0 就是 “北方不吃米飯”，是小明父母不希望，想拒絕的。
- H1又叫備選假設，一般來說，我們都希望備選假設，也就是H1為真，比如小明的父母希望北方也是吃飯的。
H0和H1並不是不能互換，但是選取的H0原則是H0必須是一個可以被拒絕的假設。對於一個假設檢驗問題，如果H1是一個不能被拒絕的假設，那麼H0和H1不能互換。
什麼樣的H1是不能被拒絕的假設？比如下面這個問題：總體是一個班級的所有同學一次的考試成績，原假設H0是全班同學的平均成績為80分，H1是全班同學的平均成績不是80分。在這個框架中，H0是一個可以被拒絕的假設，H1是一個不能被拒絕的假設。
為什麼這麼說？如果H0選做全班同學的平均成績不是80分，那麼哪怕你知道全班幾乎所有人的成績，不妨認為平均是80分，只要有一個人的成績你不知道，那麼知道其他人的成績對於你拒絕H0沒有任何幫助。只有在這個人恰好80分的情況下，原假設是不成立的，這個人是81或者79,原假設都成立。從另一個角度，也可以認為當H0是一個不能被拒絕的假設的時候，H0太過寬泛，和H1區分不了。
一般來說，假設檢驗的結果有下面兩個：
- H0錯，H1對（用專業的話來說，就是接受H1，拒絕H0，因為有足夠多的樣本支援H1，比如說小明爸媽問了5，6個到北方讀書的人，都說食堂面和飯都有。他們大可以說，北方也是吃飯的）。
- H1錯（用專業的話說，不拒絕H0，因為證據不夠。小明爸媽問了5個到北方讀書的人，1個說食堂面和飯都有，那食堂到底有沒有呢？心裡可是直打鼓了。）
有一點要注意，無法拒絕H0並不是說H0為真，而是說我們的證據不足，無法證明H1。我們經常在法庭上聽到，證據不足，無罪釋放，但是，這個人到底有沒有罪，還是要打個問號的。
1.3 P值 (P-value)
P值，也就是常見到的 P-value。P 值是一種概率，指的是在 H0 假設為真的前提下，樣本結果出現的概率。即 p值是在原假設成立的基礎上計算的。
如果 P-value 很小，則說明在原假設為真的前提下，樣本結果出現的概率很小，甚至很極端，這就反過來說明了原假設很大概率是錯誤的。
另外一個角度想，p值也是棄真錯誤的概率。也就是這個假設成立的情況下，出現這個糟糕結果的概率，也當然就是如果此時把h0拒絕，出現錯誤的概率。只要p值足夠小，我們就認為此時拒絕h0，出錯的概率很小，那就乾脆把h0拒絕好了。
通常，會設定一個顯著性水平（significance level） alpha 與 P-value 進行比較，如果 P-value < alpha ，則說明在顯著性水平 alpha 下拒絕原假設，alpha 通常情況下設定為0.05。
假如我們比較某地區男、女性的飲食口味是否存在差異，則 H0 是 "男女的飲食口味相同，不存在差異"。
最後得出 P=0.283 > 0.05，在α=0.05水平上不拒絕零假設，即不能認為該地區男女的飲食口味不同。
1.4 交叉表
在統計學中，交叉表是矩陣格式的一種表格，顯示變數的（多變數）頻率分佈。
“交叉表”物件是一個網格，用來根據指定的條件返回值。資料顯示在壓縮行和列中。這種格式易於比較資料並辨別其趨勢。它由三個元素組成：行 / 列 / 摘要欄位。
讓我們舉例說明：
| 馬軍頭領 | 步兵頭領 | |
|---|---|---|
| 二龍山 | 1 | 6 |
| 少華山 | 3 | 0 |
- “交叉表”中的行沿水平方向延伸（從一側到另一側）。在上面的示例中，”二龍山” 是一行。
- “交叉表”中的列沿垂直方向延伸（上下）。在上面的示例中，“馬軍頭領” 是一列。
- 彙總欄位位於行和列的交叉處。每個交叉處的值代表對既滿足行條件又滿足列條件的記錄的彙總（求和、計數等）。在上面的示例中，“二龍山”和“馬軍頭領”交叉處的值是1，這是在二龍山上馬軍頭領的數目。
上文中交叉表是按兩個變數交叉分類的，該列聯表稱為兩維列聯表，若按3個變數交叉分類，所得的列聯表稱為3維列聯表，依次類推。3維及以上的列聯表通常稱為“多維列聯表”或“高維列聯表”，而一維列聯表就是頻數分佈表。
1.5 卡方
交叉分類所得的表格稱為“列聯表”，統計推斷（檢驗）則要使用列聯表分析的方法------卡方檢驗。
卡方檢驗，主要用於檢驗統計樣本的實際觀測值與理論推斷值之間的偏離程度，或者是檢驗一批資料是否與某種理論分佈相符合。
Alink 文件中給出的是：卡方獨立性檢驗是檢驗兩個因素（各有兩項或以上的分類）之間是否相互影響的問題，其零假設是兩因素之間相互獨立。
1.5.1 公式
卡方值是卡方檢驗時用到的檢驗統計量，卡方值越大，說明觀測值與理論值之間的偏離就越大；反之，二者偏差越小。實際應用時，可以根據卡方值計算 P-value，從而選擇拒絕或者接受原假設。
公式如下：
1.5.2 基本思想
卡方檢驗最基本的思想就是通過觀察實際值與理論值的偏差來確定理論的正確與否。
具體做的時候常常先假設兩個變數確實是獨立的（行話就叫做“原假設”），然後觀察實際值（也可以叫做觀察值）與理論值（這個理論值是指“如果兩者確實獨立”的情況下應該有的值）的偏差程度。
- 如果偏差足夠小，我們就認為誤差是很自然的樣本誤差，是測量手段不夠精確導致或者偶然發生的，兩者確確實實是獨立的，此時就接受原假設。
- 如果偏差大到一定程度，使得這樣的誤差不太可能是偶然產生或者測量不精確所致，我們就認為兩者實際上是相關的，即否定原假設，而接受備擇假設。
1.5.3 實現過程
卡方分析的方法：
- 假設兩個變數是相互獨立，互不關聯的。這在統計上稱為原假設；
- 對於調查中得到的兩個變數的資料，用一個表格的形式來表示它們的分佈（頻數和百分數），這裡的頻數叫觀測頻數，這種表格叫列聯表；
- 如果原假設成立，在這個前提下，可以計算出上面列聯表中每個格子裡的頻數應該是多少，這叫期望頻數；
- 比較觀測頻數與期望頻數的差，如果兩者的差越大，表明實際情況與原假設相去甚遠；差越小，表明實際情況與原假設越相近。這種差值用一個卡方統計量來表示；
- 對卡方值進行檢驗，如果卡方檢驗的結果不顯著，則不能拒絕原假設，即兩變數是相互獨立、互不關聯的，如果卡方檢驗的結果顯著，則拒絕原假設，即兩變數間存在某種關聯，至於是如何關聯的，這要看列聯表中資料的分佈形態。
具體實現過程：
- 按照假設檢驗的步驟，首先我們需要確定原假設 H0(null hypothesis)：原假設是變數獨立的，實際觀測頻率和理論頻率一致。
- 其次我們根據實際觀測的聯連表，去求理論的聯連表；卡方統計值：X2，記為Statistic；
- 然後選取適合的置信度(一般為95%)同自由度一起確定臨界值Critical Value，比較卡方統計值和臨界值大小：
- If Statistic >= Critical Value: 認為變數對結果有影響，則拒絕原假設，變數不獨立
- If Statistic < Critical Value: 認為變數對結果沒有影響，接受原假設，變數獨立
1.6 自由度
自由度：取值不受限制的變數的個數。
如何理解這句簡單的話呢？給定一組資料，我們來計算不同的統計量，看看自由度的變化。這些資料分別為 1 2 4 6 8. 5個數。
先來求平均值，這幾個資料都可以任意變化成其它資料，而我們仍然可以對它們求平均值，它們的平均值也跟著變化。這時自由度為5，也就是說有幾個資料自由度就是幾。
卡方檢驗的自由度：
1）如果是獨立性檢驗，那麼自由度就等於（a-1）*（b-1），a b表示這兩個檢驗條件的對應的分類數。
2）適合性檢驗，類別數減去1。此處相當於約束條件只有一個。
卡方檢驗只有在用筆算查表時使用自由度，軟體計算不用擔心這個問題，但是最好明白自由度代表著總的變數數目減去約束條件的數目。
0x02 示例程式碼
本文示例程式碼如下，這裡需要注意的是：
- "col1","col2"是所選擇的列；
- "col4"是Label；
public class ChiSquareTestBatchOpExample {
public static void main(String[] args) throws Exception {
Row[] testArray =
new Row[]{
Row.of("a", 1.1, 1.2, 1),
Row.of("b", 0.9, 1.0, -2),
Row.of("c", -0.01, 1.0, 100),
Row.of("d", 100.9, 0.1, -99),
Row.of("a", 1.1, 1.2, 1),
Row.of("b", 0.9, 1.0, -2),
Row.of("c", -0.01, 0.2, 100),
Row.of("d", 100.9, 0.3, -99)
};
String[] colNames = new String[]{"col1", "col2", "col3", "col4"};
MemSourceBatchOp source = new MemSourceBatchOp(Arrays.asList(testArray), colNames);
ChiSquareTestBatchOp test = new ChiSquareTestBatchOp()
.setSelectedCols("col1","col2")
.setLabelCol("col4");
test.linkFrom(source).print();
}
}
輸出如下：
col|chisquare_test
---|--------------
col1|{"comment":"chi-square test","df":9.0,"p":0.004301310843500827,"value":24.0}
col2|{"comment":"chi-square test","df":9.0,"p":0.004301310843500827,"value":24.0}
轉換為圖表更好理解：
| col | chisquare_test |
|---|---|
| col1 | {"comment":"chi-square test","df":9.0,"p":0.004301310843500827,"value":24.0} |
| col2 | {"comment":"chi-square test","df":9.0,"p":0.004301310843500827,"value":24.0} |
df是自由度，p就是p-value， value就是我們前面說的卡方值，即
0x03 總體邏輯
訓練總體邏輯如下：
- 使用 flatMap 做 flatting data to triple。遍歷輸入Row，然後把Row給flat了，得到三元組<idx in row, value in row, y-label>。比如 對應輸入 Row.of("b", 0.9, 1.0, -2)，則row = {Row@9419} "b,0,9,-2"，因為col1, col2是特徵，col4是 label，則傳送兩個三元組是 <0, b, -2>, <1, 0.9, -2>；
- 使用 toTable 把前面處理的dataSet再進行轉換，生成一張表 data。{"col", "feature", "label"} 就對應著我們之前的三元組；
- 對 data 進行 計算交叉表 和 卡方校驗；
- groupBy("col,feature,label") 進行分類排序；
- select("col,feature,label,count(1) as count2")) 得出 feature 的個數作為count2；
- groupBy("col").reduceGroup 再根據col排序，歸併；
- 得到 <feature, y-label> : "count of feature" 這個map；
- Crosstab.convert(map) 利用map來做交叉表；
- map(new ChiSquareTestFromCrossTable()) 利用交叉表來構建卡方檢驗；
- test(crossTabWithId) 這裡進行計算，其中會呼叫 org.apache.commons.math3.distribution.GammaDistribution.cumulativeProbability 進行Gamma計算；
0x04 訓練
還是老套路，直奔ChiSquareTestBatchOp的linkFrom函式。
程式碼是縮減版，但原本就非常簡單，獲取“選擇的列”和“Y列”，然後用輸入資料進行訓練檢驗。
深入看下去卻很有難度。
public ChiSquareTestBatchOp linkFrom(BatchOperator<?>... inputs) {
BatchOperator<?> in = checkAndGetFirst(inputs);
String[] selectedColNames = getSelectedCols();
String labelColName = getLabelCol();
this.setOutputTable(ChiSquareTestUtil.buildResult(
ChiSquareTestUtil.test(in, selectedColNames, labelColName),
selectedColNames,
getMLEnvironmentId()));
return this;
}
最後會輾轉進入到 ChiSquareTest.test，這裡才是真章。
public static DataSet<Row> test(BatchOperator in,
String[] selectedColNames,
String labelColName) {
in = in.select(ArrayUtils.add(selectedColNames, labelColName));
return ChiSquareTest.test(in.getDataSet(), in.getMLEnvironmentId());
}
4.1 ChiSquareTest
其test函式的輸入輸出是：
- 輸入：in 的最後一列是label，其餘列是所選擇的特徵列；
- 輸出：有三列，分別是 1th is colId, 2th is pValue, 3th is chi-square value；
這裡的總體邏輯是：
- 使用 flatMap 做 flatting data to triple。遍歷輸入Row，然後把Row給flat了，得到三元組<idx in row, value in row, y-label>；
- 使用 toTable 把前面處理的dataSet再進行轉換，生成一張表 data。{"col", "feature", "label"} 就對應著我們之前的三元組；
- 對 data 進行 計算交叉表 和 卡方校驗；
具體程式碼如下：
protected static DataSet<Row> test(DataSet<Row> in, Long sessionId) {
//flatting data to triple.
//這裡就是遍歷輸入Row，然後把Row給flat了，得到三元組<idx in row, value in row, y-label>
//比如 對應輸入 Row.of("b", 0.9, 1.0, -2)，則row = {Row@9419} "b,0,9,-2"，因為col1, col2是特徵，col4是 label，則傳送兩個三元組是 <0, b, -2>, <1, 0.9, -2>
DataSet<Row> dataSet = in
.flatMap(new FlatMapFunction<Row, Row>() {
@Override
public void flatMap(Row row, Collector<Row> result) {
int n = row.getArity() - 1;
String nStr = String.valueOf(row.getField(n));
for (int i = 0; i < n; i++) {
Row out = new Row(3);
out.setField(0, i);
out.setField(1, String.valueOf(row.getField(i)));
out.setField(2, nStr);
result.collect(out);
}
}
});
// 把前面處理的dataSet再進行轉換，生成一張表。{"col", "feature", "label"} 就對應著我們之前的三元組
Table data = DataSetConversionUtil.toTable(
sessionId,
dataSet,
new String[]{"col", "feature", "label"},
new TypeInformation[]{Types.INT, Types.STRING, Types.STRING});
// 對 data 進行 計算交叉表 和 卡方校驗
//calculate cross table and chiSquare test.
return DataSetConversionUtil.fromTable(sessionId, data
.groupBy("col,feature,label") //分類排序
.select("col,feature,label,count(1) as count2")) // 為了得出 feature 的個數作為count2
.groupBy("col").reduceGroup( // 再根據col排序
new GroupReduceFunction<Row, Tuple2<Integer, Crosstab>>() {
@Override
public void reduce(Iterable<Row> iterable, Collector<Tuple2<Integer, Crosstab>> collector) {
Map<Tuple2<String, String>, Long> map = new HashMap<>();
int colIdx = -1;
for (Row row : iterable) {
// 假如有如下，row = {Row@9684} "0,a,1,2"，他對應了兩個 Row.of("a", 1.1, 1.2, 1), 就是 <col,feature,label,count(1)>, 就是 <'a'是第0列，'a'，對應 y-label是 1， 'a' 有兩個>
map.put(Tuple2.of(row.getField(1).toString(),
row.getField(2).toString()),
(long) row.getField(3));
colIdx = (Integer) row.getField(0);
}
// 得到 <feature, y-label> : "count of feature" 這個map
map = {HashMap@9676} size = 4
{Tuple2@9688} "(a,1)" -> {Long@9689} 2
{Tuple2@9690} "(b,-2)" -> {Long@9689} 2
{Tuple2@9691} "(d,-99)" -> {Long@9689} 2
{Tuple2@9692} "(c,100)" -> {Long@9689} 2
// 利用map來做交叉表
collector.collect(new Tuple2<>(colIdx, Crosstab.convert(map)));
}
})
.map(new ChiSquareTestFromCrossTable()); // 構建卡方檢驗
}
4.2 Crosstab
上面程式碼中，使用 collector.collect(new Tuple2<>(colIdx, Crosstab.convert(map))); 來構建交叉表。
Crosstab 就是 Cross Tabulations reflects the relationship between two variables。即以map key為橫軸，縱軸，value作為數值，就是feature和label之間的交叉。
public static Crosstab convert(Map<Tuple2<String, String>, Long> maps) {
Crosstab crosstab = new Crosstab();
//get row tags and col tags
Set<Tuple2<String, String>> sets = maps.keySet();
Set<String> rowTags = new HashSet<>(); // 拿到行，列
Set<String> colTags = new HashSet<>();
for (Tuple2<String, String> tuple2 : sets) {
rowTags.add(tuple2.f0);
colTags.add(tuple2.f1);
}
crosstab.rowTags = new ArrayList<>(rowTags);
crosstab.colTags = new ArrayList<>(colTags);
int rowLen = crosstab.rowTags.size();
int colLen = crosstab.colTags.size();
//compute value
crosstab.data = new long[rowLen][colLen];
for (Map.Entry<Tuple2<String, String>, Long> entry : maps.entrySet()) {
int rowIdx = crosstab.rowTags.indexOf(entry.getKey().f0);
int colIdx = crosstab.colTags.indexOf(entry.getKey().f1);
crosstab.data[rowIdx][colIdx] = entry.getValue();
}
return crosstab;
}
這裡輸入輸出如下
// 輸入如下
maps = {HashMap@9676} size = 4
{Tuple2@9688} "(a,1)" -> {Long@9689} 2
{Tuple2@9690} "(b,-2)" -> {Long@9689} 2
{Tuple2@9691} "(d,-99)" -> {Long@9689} 2
{Tuple2@9692} "(c,100)" -> {Long@9689} 2
// 交叉表如下
crosstab = {Crosstab@9703}
colTags = {ArrayList@9720} size = 4
0 = "1" 1 = "100" 2 = "-2" 3 = "-99"
rowTags = {ArrayList@9721} size = 4
0 = "a" 1 = "b" 2 = "c" 3 = "d"
data = {long[4][]@9713}
0 = {long[4]@9722} 0 = 2 1 = 0 2 = 0 3 = 0
1 = {long[4]@9723} 0 = 0 1 = 0 2 = 2 3 = 0
2 = {long[4]@9724} 0 = 0 1 = 2 2 = 0 3 = 0
3 = {long[4]@9725} 0 = 0 1 = 0 2 = 0 3 = 2
構造出來交叉表如下：
| 1 | 100 | -2 | -99 | |
|---|---|---|---|---|
| a | 2 | |||
| b | 2 | |||
| c | 2 | |||
| d | 2 |
4.3 構建卡方檢驗
4.1中，有 .map(new ChiSquareTestFromCrossTable()); ，這裡就是根據collector.collect(new Tuple2<>(colIdx, Crosstab.convert(map)));交叉表構建卡方檢驗。
/**
* calculate chi-square test value from cross table.
*/
public static class ChiSquareTestFromCrossTable implements MapFunction<Tuple2<Integer, Crosstab>, Row> {
@Override
public Row map(Tuple2<Integer, Crosstab> crossTabWithId) throws Exception {
Tuple4 tuple4 = test(crossTabWithId);
// f0 is id of cross table, f1 is pValue, f2 is chi-square Value, f3 is df
Row row = new Row(4);
row.setField(0, tuple4.f0);
row.setField(1, tuple4.f1);
row.setField(2, tuple4.f2);
row.setField(3, tuple4.f3);
return row;
}
}
test(crossTabWithId)是關鍵點，其中 distribution.cumulativeProbability 最後呼叫到 org.apache.commons.math3.distribution.GammaDistribution.cumulativeProbability。
這裡能夠看到
- df 的定義就是 (double)(rowLen - 1) * (colLen - 1)，即（行 - 1）*（列 - 1）。
- 卡方值就是嚴格按照定義來構建的。
- p-value是 呼叫到 org.apache.commons.math3.distribution.GammaDistribution.cumulativeProbability。
/**
* @param crossTabWithId: f0 is id, f1 is cross table
* @return tuple4: f0 is id which is id of cross table, f1 is pValue, f2 is chi-square Value, f3 is df
*/
protected static Tuple4<Integer, Double, Double, Double> test(Tuple2<Integer, Crosstab> crossTabWithId) {
int colIdx = crossTabWithId.f0;
Crosstab crosstab = crossTabWithId.f1;
int rowLen = crosstab.rowTags.size();
int colLen = crosstab.colTags.size();
//compute row sum and col sum 計算出列的數值和，行的數值和
double[] rowSum = crosstab.rowSum();
double[] colSum = crosstab.colSum();
double n = crosstab.sum();
//compute statistic value 計算統計值
double chiSq = 0;
for (int i = 0; i < rowLen; i++) {
for (int j = 0; j < colLen; j++) {
double nij = rowSum[i] * colSum[j] / n;
double temp = crosstab.data[i][j] - nij;
chiSq += temp * temp / nij; // 就是按照定義來構建卡方值
}
}
//set result
double p;
if (rowLen <= 1 || colLen <= 1) {
p = 1;
} else {
ChiSquaredDistribution distribution =
new ChiSquaredDistribution(null, (rowLen - 1) * (colLen - 1));
p = 1.0 - distribution.cumulativeProbability(Math.abs(chiSq));
}
// return tuple4: f0 is id which is id of cross table, f1 is pValue, f2 is chi-square Value, f3 is df
return Tuple4.of(colIdx, p, chiSq, (double)(rowLen - 1) * (colLen - 1));
}
// runtime是
tuple4 = {Tuple4@9842} "(0,0.004301310843500827,24.0,9.0)"
f0 = {Integer@9843} 0
f1 = {Double@9844} 0.004301310843500827
f2 = {Double@9847} 24.0
f3 = {Double@9848} 9.0
0xFF 參考
卡方檢驗(Chi_square_test)： 原理及python實現
Spark MLlib基本演算法【相關性分析、卡方檢驗、總結器】
spark(1.1) mllib 原始碼分析(一)-卡方檢驗