使用scipy.stats.chi2_contingency遇ValueError求无常数修复方案

问题背景 #

使用scipy.stats.chi2_contingency实现卡方决策树时，触发以下错误：

ValueError: The internally computed table of expected frequencies has a zero element at (0, 0)

目前找到的临时解决方案是给类别计数添加1e-10小常数，但该方案属于hack性质，会干扰统计结果的准确性，需要更合理的修复方式。

错误根源 #

问题出在类别计数的生成逻辑：
原代码使用np.bincount(y_left, minlength=len(classes))计算子集类别计数，但classes是从样本标签y中提取的集合（可能不是从0开始的连续整数）。np.bincount会基于标签值的索引生成数组，导致数组中包含原始数据不存在的类别（比如标签从1开始时，索引0对应的计数为0）。这些不存在的类别在左右子集的计数均为0，形成列联表中全0的列，计算期望频数时就会出现0值，触发scipy的报错。

正确修复方案 #

放弃np.bincount，直接针对classes中的每个真实类别统计子集内的样本数，确保列联表只包含实际存在的类别，避免全0列的出现。

修改代码片段 #

在best_split函数中，替换原有的left_counts和right_counts计算代码：

# 替换原bincount实现
left_counts = np.array([np.sum(y_left == c) for c in classes])
right_counts = np.array([np.sum(y_right == c) for c in classes])

修复原理 #

• 仅统计classes中存在的真实类别，不会引入不存在的类别导致全0列
• 每个类别的列和为原始数据中该类别的样本数（必然大于0），行和为左/右子集的样本数（拆分时保证子集至少有1个样本，必然大于0）
• 所有期望频数E = (行和 * 列和)/总样本数均为正数，完全满足卡方检验的前提条件

修改后的完整代码 #

import numpy as np
from scipy.stats import chi2_contingency
import matplotlib.pyplot as plt

class Node_chi2:
    def __init__(self, chi2, num_samples, num_samples_per_class, predicted_class):
        self.chi2 = chi2
        self.num_samples = num_samples
        self.num_samples_per_class = num_samples_per_class
        self.predicted_class = predicted_class
        self.feature_index = 0
        self.threshold = 0
        self.left = None
        self.right = None

class DecisionTree_chi2():

    def __init__(self, max_depth=None):
        self.max_depth = max_depth

    def fit(self, X, y, depth=0):
        classes = list(set(y))
        num_samples_per_class = [np.sum(y == i) for i in classes]
        predicted_class = classes[np.argmax(num_samples_per_class)]
        node = Node_chi2(
              chi2=0,  # chi2 is 0 at the start
              num_samples=len(y),
              num_samples_per_class=num_samples_per_class,
              predicted_class=predicted_class,
          )

        if depth < self.max_depth:
            idx, thr = self.best_split(X, y, classes)
            if idx is not None:
                indices_left = X[:, idx] < thr
                X_left, y_left = X[indices_left], y[indices_left]
                X_right, y_right = X[~indices_left], y[~indices_left]
                node.feature_index = idx
                node.threshold = thr
                node.left = self.fit(X_left, y_left, depth + 1)
                node.right = self.fit(X_right, y_right, depth + 1)

        self.node = node
        return self.node
    
    def best_split(self, X, y, classes):
        m, n = X.shape
        if m <= 1:
            return None, None

        best_chi2 = 0  # Best chi2 score
        best_idx, best_thr = None, None

        for idx in range(n):
            thresholds, classes_sorted = zip(*sorted(zip(X[:, idx], y)))
            for i in range(1, m):  # iterate through each threshold
                if thresholds[i] == thresholds[i - 1]:
                    continue
                y_left = np.array(classes_sorted[:i])
                y_right = np.array(classes_sorted[i:])
                # 修复后的类别计数逻辑
                left_counts = np.array([np.sum(y_left == c) for c in classes])
                right_counts = np.array([np.sum(y_right == c) for c in classes])
                chi2, p, _ , _ = chi2_contingency(np.array([left_counts, right_counts]))

                if chi2 > best_chi2:
                    best_chi2 = chi2
                    best_idx = idx
                    best_thr = (thresholds[i] + thresholds[i - 1]) / 2

        return best_idx, best_thr

    def predict(self, X):
        yhat = []
        for sample in X: 
            node = self.node
            while node.left:
                if sample[node.feature_index] < node.threshold:
                    node = node.left
                else:
                    node = node.right
            yhat.append(node.predicted_class)
        return np.array(yhat)

from sklearn.datasets import make_moons, make_classification

X, y = make_moons(n_samples=100, noise=0.1, random_state=42)

def plot_decision_boundaries(clf, X, y, label=''):
    # Generate a grid of points to make predictions:
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, 0.1), np.arange(x2_min, x2_max, 0.1))

    # Use the classifier to make predictions on the grid:
    Z = clf.predict(np.c_[xx1.ravel(), xx2.ravel()])

    Z = Z.reshape(xx1.shape)

    colors = {0: 'C0', 1: 'C1', 2: 'C2'}
    # colors dictionary indexed by class label

    # Create a contour plot to display the decision boundaries:
    fig = plt.figure(figsize=(4, 4))

    from matplotlib.colors import ListedColormap
    n_classes = len(np.unique(y))
    custom_cmap = ListedColormap(list(colors.values())[:n_classes])

    plt.contourf(xx1, xx2, Z, cmap=custom_cmap, alpha=0.4)

    plt.scatter(X[y==0,0], X[y==0,1], c=colors[0], marker='o', edgecolors='k', alpha=0.5)
    plt.scatter(X[y==1,0], X[y==1,1], c=colors[1], marker='o', edgecolors='k', alpha=0.5)
    plt.scatter(X[y==2,0], X[y==2,1], c=colors[2], marker='o', edgecolors='k', alpha=0.5)

    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.title('Decision Boundaries - '+label+' DT')
    plt.show()


clf_gini = DecisionTree_chi2(max_depth=5)
clf_gini.fit(X, y)

plot_decision_boundaries(clf_gini, X, y, label='Chi2')

关于临时方案的说明 #

添加1e-10的临时方案虽然能绕过报错，但会人为修改观测频数，影响卡方值的计算准确性，尤其是样本量较小时，误差会被放大，因此不推荐使用。

内容的提问来源于stack exchange，提问作者Visage