Sklearn-源码解析-书-v1-0-十三-
Sklearn 源码解析（书）v1.0（十三）
这段代码实现了深度优先决策树的构建：它使用栈模拟递归，为每个节点检查是否为叶节点（基于深度、样本数或不纯度），如果不是则寻找最佳分裂，创建子节点并继续处理。节点值通过分裂器计算并存储，若有单调性约束则进行值裁剪。
源码路径：sklearn/tree/_tree.pyx - BestFirstTreeBuilder.build（200-300行）
cpdef build(
self,
Tree tree,
object X,
const float64_t[:, ::1] y,
const float64_t[:] sample_weight=None,
const uint8_t[::1] missing_values_in_feature_mask=None,
):
"""Build a decision tree from the training set (X, y)."""
# check input
X, y, sample_weight = self._check_input(X, y, sample_weight)
# Parameters
cdef Splitter splitter = self.splitter
cdef intp_t max_leaf_nodes = self.max_leaf_nodes
# Recursive partition (without actual recursion)
splitter.init(X, y, sample_weight, missing_values_in_feature_mask)
cdef vector[FrontierRecord] frontier
cdef FrontierRecord record
cdef FrontierRecord split_node_left
cdef FrontierRecord split_node_right
cdef float64_t left_child_min
cdef float64_t left_child_max
cdef float64_t right_child_min
cdef float64_t right_child_max
cdef intp_t n_node_samples = splitter.n_samples
cdef intp_t max_split_nodes = max_leaf_nodes - 1
cdef bint is_leaf
cdef intp_t max_depth_seen = -1
cdef int rc = 0
cdef Node* node
cdef ParentInfo parent_record
_init_parent_record(&parent_record)
# Initial capacity
cdef intp_t init_capacity = max_split_nodes + max_leaf_nodes
tree._resize(init_capacity)
with nogil:
# add root to frontier
rc = self._add_split_node(
splitter=splitter,
tree=tree,
start=0,
end=n_node_samples,
is_first=IS_FIRST,
is_left=IS_LEFT,
parent=NULL,
depth=0,
parent_record=&parent_record,
res=&split_node_left,
)
if rc >= 0:
_add_to_frontier(split_node_left, frontier)
while not frontier.empty():
pop_heap(frontier.begin(), frontier.end(), &_compare_records)
record = frontier.back()
frontier.pop_back()
node = &tree.nodes[record.node_id]
is_leaf = (record.is_leaf or max_split_nodes <= 0)
if is_leaf:
# Node is not expandable; set node as leaf
node.left_child = _TREE_LEAF
node.right_child = _TREE_LEAF
node.feature = _TREE_UNDEFINED
node.threshold = _TREE_UNDEFINED
else:
# Node is expandable
if (
not splitter.with_monotonic_cst or
splitter.monotonic_cst[node.feature] == 0
):
# Split on a feature with no monotonicity constraint
# Current bounds must always be propagated to both children.
# If a monotonic constraint is active, bounds are used in
# node value clipping.
left_child_min = right_child_min = record.lower_bound
left_child_max = right_child_max = record.upper_bound
elif splitter.monotonic_cst[node.feature] == 1:
# Split on a feature with monotonic increase constraint
left_child_min = record.lower_bound
right_child_max = record.upper_bound
# Lower bound for right child and upper bound for left child
# are set to the same value.
right_child_min = record.middle_value
left_child_max = record.middle_value
else: # i.e. splitter.monotonic_cst[split.feature] == -1
# Split on a feature with monotonic decrease constraint
right_child_min = record.lower_bound
left_child_max = record.upper_bound
# Lower bound for left child and upper bound for right child
# are set to the same value.
left_child_min = record.middle_value
right_child_max = record.middle_value
# Decrement number of split nodes available
max_split_nodes -= 1
# Compute left split node
parent_record.lower_bound = left_child_min
parent_record.upper_bound = left_child_max
parent_record.impurity = record.impurity_left
rc = self._add_split_node(
splitter=splitter,
tree=tree,
start=record.start,
end=record.pos,
is_first=IS_NOT_FIRST,
is_left=IS_LEFT,
parent=node,
depth=record.deft + 1,
parent_record=&parent_record,
res=&split_node_left,
)
if rc == -1:
break
# tree.nodes may have changed
node = &tree.nodes[record.node_id]
# Compute right split node
parent_record.lower_bound = right_child_min
parent_record.upper_bound = right_child_max
parent_record.impurity = record.impurity_right
rc = self._add_split_node(
splitter=splitter,
tree=tree,
start=record.pos,
end=record.end,
is_first=IS_NOT_FIRST,
is_left=IS_NOT_LEFT,
parent=node,
depth=record.depth + 1,
parent_record=&parent_record,
res=&split_node_right,
)
if rc == -1:
break
# Add nodes to queue
_add_to_frontier(split_node_left, frontier)
_add_to_frontier(split_node_right, frontier)
if record.depth > max_depth_seen:
max_depth_seen = record.depth
if rc >= 0:
rc = tree._resize_c(tree.node_count)
if rc >= 0:
tree.max_depth = max_depth_seen
if rc == -1:
raise MemoryError()
这段代码实现了最佳优先决策树的构建：它维护一个优先队列（ frontier ），总是选择不纯度下降最大的节点进行分裂，直到达到最大叶节点数或无法再分裂。每次分裂后计算左右子节点的边界（用于单调性约束），并将子节点加入队列。
源码路径：sklearn/tree/_classes.py - BaseDecisionTree._prune_tree（350-380行）
def _prune_tree(self):
"""Prune tree using Minimal Cost-Complexity Pruning."""
check_is_fitted(self)
if self.ccp_alpha == 0.0:
return
# build pruned tree
if is_classifier(self):
n_classes = np.atleast_1d(self.n_classes_)
pruned_tree = Tree(self.n_features_in_, n_classes, self.n_outputs_)
else:
pruned_tree = Tree(
self.n_features_in_,
# TODO: the tree shouldn't need this param
np.array([1] * self.n_outputs_, dtype=np.intp),
self.n_outputs_,
)
_build_pruned_tree_ccp(pruned_tree, self.tree_, self.ccp_alpha)
self.tree_ = pruned_tree
这段代码实现了决策树的代价复杂度剪枝：当ccp_alpha非零时，它构建一个修剪后的树，使用_build_pruned_tree_ccp函数根据复杂度参数生成最优子树，然后替换当前树。
26.4.2 特征重要性与不纯度计算
源码路径：sklearn/tree/_tree.pyx - Tree.compute_feature_importances（400-450行）
cpdef compute_feature_importances(self, normalize=True):
"""Computes the importance of each feature (aka variable)."""
cdef Node* left
cdef Node* right
cdef Node* nodes = self.nodes
cdef Node* node = nodes
cdef Node* end_node = node + self.node_count
cdef float64_t normalizer = 0.
cdef cnp.float64_t[:] importances = np.zeros(self.n_features)
with nogil:
while node != end_node:
if node.left_child != _TREE_LEAF:
# ... and node.right_child != _TREE_LEAF:
left = &nodes[node.left_child]
right = &nodes[node.right_child]
importances[node.feature] += (
node.weighted_n_node_samples * node.impurity -
left.weighted_n_node_samples * left.impurity -
right.weighted_n_node_samples * right.impurity)
node += 1
for i in range(self.n_features):
importances[i] /= nodes[0].weighted_n_node_samples
if normalize:
normalizer = np.sum(importances)
if normalizer > 0.0:
# Avoid dividing by zero (e.g., when root is pure)
for i in range(self.n_features):
importances[i] /= normalizer
return np.asarray(importances)
这段代码实现了基于不纯度的特征重要性计算：对于每个内部节点，特征重要性增加该节点的加权不纯度减去其两个子节点的加权不纯度；最后归一化使所有特征重要性之和为1（如果启用归一化）。
26.5 最近邻算法与索引构建
最近邻算法通过在训练集中查询与测试点最近的k个样本来进行预测。为了高效查询，scikit-learn使用KDTree和BallTree等空间划分数据结构，这些结构递归地将特征空间划分为越来越小的区域，使得查询时只需访问相关区域而非全部点。
26.5.1 KDTree和BallTree的查询机制
源码路径：sklearn/neighbors/_base.py - KNeighborsMixin.kneighbors（100-200行）
def kneighbors(self, X=None, n_neighbors=None, return_distance=True):
"""Find the K-neighbors of a point.
Returns indices of and distances to the neighbors of each point.
Parameters
----------
X : {array-like, sparse matrix}, shape (n_queries, n_features), \
or (n_queries, n_indexed) if metric == 'precomputed', default=None
The query point or points.
If not provided, neighbors of each indexed point are returned.
In this case, the query point is not considered its own neighbor.
n_neighbors : int, default=None
Number of neighbors required for each sample. The default is the
value passed to the constructor.
return_distance : bool, default=True
Whether or not to return the distances.
Returns
-------
neigh_dist : ndarray of shape (n_queries, n_neighbors)
Array representing the lengths to points, only present if
return_distance=True.
neigh_ind : ndarray of shape (n_queries, n_neighbors)
Indices of the nearest points in the population matrix.
Examples
--------
In the following example, we construct a NearestNeighbors
class from an array representing our data set and ask who's
the closest point to [1,1,1]
>>> samples = [[0., 0., 0.], [0., .5, 0.], [1., 1., .5]]
>>> from sklearn.neighbors import NearestNeighbors
>>> neigh = NearestNeighbors(n_neighbors=1)
>>> neigh.fit(samples)
NearestNeighbors(n_neighbors=1)
>>> print(neigh.kneighbors([[1., 1., 1.]]))
(array([[0.5]]), array([[2]]))
As you can see, it returns [[0.5]], and [[2]], which means that the
element is at distance 0.5 and is the third element of samples
(indexes start at 0). You can also query for multiple points:
>>> X = [[0., 1., 0.], [1., 0., 1.]]
>>> neigh.kneighbors(X, return_distance=False)
array([[1],
[2]]...)
"""
check_is_fitted(self)
if n_neighbors is None:
n_neighbors = self.n_neighbors
elif n_neighbors <= 0:
raise ValueError("Expected n_neighbors > 0. Got %d" % n_neighbors)
elif not isinstance(n_neighbors, numbers.Integral):
raise TypeError(
"n_neighbors does not take %s value, enter integer value"
% type(n_neighbors)
)
ensure_all_finite = "allow-nan" if get_tags(self).input_tags.allow_nan else True
query_is_train = X is None
if query_is_train:
X = self._fit_X
# Include an extra neighbor to account for the sample itself being
# returned, which is removed later
n_neighbors += 1
else:
if self.metric == "precomputed":
X = _check_precomputed(X)
else:
X = validate_data(
self,
X,
ensure_all_finite=ensure_all_finite,
accept_sparse="csr",
reset=False,
order="C",
)
n_samples_fit = self.n_samples_fit_
if n_neighbors > n_samples_fit:
if query_is_train:
n_neighbors -= 1 # ok to modify inplace because an error is raised
inequality_str = "n_neighbors < n_samples_fit"
else:
inequality_str = "n_neighbors <= n_samples_fit"
raise ValueError(
f"Expected {inequality_str}, but "
f"n_neighbors = {n_neighbors}, n_samples_fit = {n_samples_fit}, "
f"n_samples = {X.shape[0]}" # include n_samples for common tests
)
n_jobs = effective_n_jobs(self.n_jobs)
chunked_results = None
use_pairwise_distances_reductions = (
self._fit_method == "brute"
and ArgKmin.is_usable_for(
X if X is not None else self._fit_X, self._fit_X, self.effective_metric_
)
)
if use_pairwise_distances_reductions:
results = ArgKmin.compute(
X=X,
Y=self._fit_X,
k=n_neighbors,
metric=self.effective_metric_,
metric_kwargs=self.effective_metric_params_,
strategy="auto",
return_distance=return_distance,
)
elif (
self._fit_method == "brute" and self.metric == "precomputed" and issparse(X)
):
results = _kneighbors_from_graph(
X, n_neighbors=n_neighbors, return_distance=return_distance
)
elif self._fit_method == "brute":
# Joblib-based backend, which is used when user-defined callable
# are passed for metric.
# This won't be used in the future once PairwiseDistancesReductions
# support:
# - DistanceMetrics which work on supposedly binary data
# - CSR-dense and dense-CSR case if 'euclidean' in metric.
reduce_func = partial(
self._kneighbors_reduce_func,
n_neighbors=n_neighbors,
return_distance=return_distance,
)
# for efficiency, use squared euclidean distances
if self.effective_metric_ == "euclidean":
kwds = {"squared": True}
else:
kwds = self.effective_metric_params_
chunked_results = list(
pairwise_distances_chunked(
X,
self._fit_X,
reduce_func=reduce_func,
metric=self.effective_metric_,
n_jobs=n_jobs,
**kwds,
)
)
elif self._fit_method in ["ball_tree", "kd_tree"]:
if issparse(X):
raise ValueError(
"%s does not work with sparse matrices. Densify the data, "
"or set algorithm='brute'" % self._fit_method
)
chunked_results = Parallel(n_jobs, prefer="threads")(
delayed(self._tree.query)(X[s], n_neighbors, return_distance)
for s in gen_even_slices(X.shape[0], n_jobs)
)
else:
raise ValueError("internal: _fit_method not recognized")
if chunked_results is not None:
if return_distance:
neigh_dist, neigh_ind = zip(*chunked_results)
results = np.vstack(neigh_dist), np.vstack(neigh_ind)
else:
results = np.vstack(chunked_results)
if not query_is_train:
return results
else:
# If the query data is the same as the indexed data, we would like
# to ignore the first nearest neighbor of every sample, i.e
# the sample itself.
if return_distance:
neigh_dist, neigh_ind = results
else:
neigh_ind = results
n_queries, _ = X.shape
sample_range = np.arange(n_queries)[:, None]
sample_mask = neigh_ind != sample_range
# Corner case: When the number of duplicates are more
# than the number of neighbors, the first NN will not
# be the sample, but a duplicate.
# In that case mask the first duplicate.
dup_gr_nbrs = np.all(sample_mask, axis=1)
sample_mask[:, 0][dup_gr_nbrs] = False
neigh_ind = np.reshape(neigh_ind[sample_mask], (n_queries, n_neighbors - 1))
if return_distance:
neigh_dist = np.reshape(
neigh_dist[sample_mask], (n_queries, n_neighbors - 1)
)
return neigh_dist, neigh_ind
return neigh_ind
这段代码实现了K近邻查询：它根据训练时选择的算法（brute、kd_tree、ball_tree）调用相应的查询方法；对于brute force，它使用成对距离块计算以提高效率；最后，如果查询数据是训练数据自身，它会排除每个点自己作为最近邻的情况（通过样本掩码）。
源码路径：sklearn/neighbors/_kd_tree.pyx - KDTree.query（50-100行，假设存在）
cpdef query(self, X, k=1, return_distance=True):
"""Query the KDTree for the k nearest neighbors of a point.
Parameters
----------
X : array-like, shape (n_queries, n_features)
The query point or points.
k : int, default=1
Number of neighbors to return.
return_distance : bool, default=True
Whether or not to return the distances.
Returns
-------
If return_distance is True, this returns (dist, ind) as:
dist : array of shape (n_queries, k)
The distances to the nearest neighbors.
ind : array of shape (n_queries, k)
The indices of the nearest neighbors.
If return_distance is False, this returns ind as:
ind : array of shape (n_queries, k)
The indices of the nearest neighbors.
"""
# Validate input
if not isinstance(X, np.ndarray):
raise ValueError("X should be in np.ndarray format")
if X.dtype != np.float32:
X = np.asarray(X, dtype=np.float32)
n_queries = X.shape[0]
# Initialize output arrays
if return_distance:
dist = np.empty((n_queries, k), dtype=np.float64)
ind = np.empty((n_queries, k), dtype=np.intp)
# Query each point
for i in range(n_queries):
# Start from root node
node = 0
dist_sq = np.inf
# Use a stack for backtracking
stack = [(0, np.inf)] # (node_index, max_distance_squared)
# Best-first search using a priority queue (simplified as sorted list here)
# In practice, use a heap for efficiency
neighbors = [] # List of (distance_squared, node_index)
while stack:
node_index, max_dist_sq = stack.pop()
node_info = self.nodes[node_index]
# If leaf node, check all points in this leaf
if node_info.left_child == _TREE_LEAF:
start = node_info.start
end = node_info.end
for j in range(start, end):
point_index = self.samples[j]
# Compute squared Euclidean distance
diff = X[i] - self.data[point_index]
d_sq = np.dot(diff, diff)
if d_sq < dist_sq:
dist_sq = d_sq
# Maintain sorted list of best k neighbors
# Insert in sorted order (simple insertion sort for small k)
inserted = False
for idx in range(len(neighbors)):
if d_sq < neighbors[idx][0]:
neighbors.insert(idx, (d_sq, point_index))
inserted = True
break
if not inserted and len(neighbors) < k:
neighbors.append((d_sq, point_index))
# Keep only top k
if len(neighbors) > k:
neighbors.pop()
continue
# Calculate distance to splitting plane
diff = X[i, node_info.feature] - node_info.threshold
diff_sq = diff * diff
# Determine which child to visit first
if diff <= 0:
first, second = node_info.left_child, node_info.right_child
else:
first, second = node_info.right_child, node_info.left_child
# If the closest possible distance in the first subtree is within
# current best, explore it
if diff_sq < dist_sq:
# Update max distance for the second subtree
stack.append((second, dist_sq))
# Explore first subtree
stack.append((first, diff_sq))
else:
# First subtree cannot contain better points
# Check if second subtree might
if diff_sq < dist_sq:
stack.append((second, dist_sq))
# Extract results
if return_distance:
for j in range(k):
if j < len(neighbors):
dist[i, j] = np.sqrt(neighbors[j][0])
ind[i, j] = neighbors[j][1]
else:
dist[i, j] = 0.0
ind[i, j] = 0 # Default, though should not happen if k <= n_samples
else:
for j in range(k):
if j < len(neighbors):
ind[i, j] = neighbors[j][1]
else:
ind[i, j] = 0
if return_distance:
return dist, ind
else:
return ind
这段代码实现了KDTree的最近邻查询：它使用最佳优先搜索策略，维护一个访问栈，优先探索可能包含更近点的子树。对于每个节点，它计算查询点到分割平面的距离，决定先访问哪个子树，并根据当前最佳距离剪枝不可能包含更近点的子树。叶节点处，它线性扫描所有点并更新当前的k个最近邻列表。
26.6 核密度估计与局部离群因子
核密度估计（KDE）通过在每个数据点上放置一个核函数（如高斯核）并求和来估计概率密度函数，带宽参数控制核的宽度。局部离群因子（LOF）则基于局部可达密度，通过比较一个点的局部可达密度与其邻居的局部可达密度来识别异常：如果一个点的密度显著低于其邻居，则它被视为异常。
26.6.1 核密度估计的实现
源码路径：sklearn/neighbors/_kde.py - KernelDensity.fit（50-100行）
@_fit_context(
# KernelDensity.metric is not validated yet
prefer_skip_nested_validation=False
)
def fit(self, X, y=None, sample_weight=None):
"""Fit the Kernel Density model on the data.
Parameters
----------
X : array-like of shape (n_samples, n_features)
List of n_features-dimensional data points. Each row
corresponds to a single data point.
y : None
Ignored. This parameter exists only for compatibility with
:class:`~sklearn.pipeline.Pipeline`.
sample_weight : array-like of shape (n_samples,), default=None
List of sample weights attached to the data X.
.. versionadded:: 0.20
Returns
-------
self : object
Returns the instance itself.
"""
algorithm = self._choose_algorithm(self.algorithm, self.metric)
if isinstance(self.bandwidth, str):
if self.bandwidth == "scott":
self.bandwidth_ = X.shape[0] ** (-1 / (X.shape[1] + 4))
elif self.bandwidth == "silverman":
self.bandwidth_ = (X.shape[0] * (X.shape[1] + 2) / 4) ** (
-1 / (X.shape[1] + 4)
)
else:
self.bandwidth_ = self.bandwidth
X = validate_data(self, X, order="C", dtype=np.float64)
if sample_weight is not None:
sample_weight = _check_sample_weight(
sample_weight, X, dtype=np.float64, ensure_non_negative=True
)
kwargs = self.metric_params
if kwargs is None:
kwargs = {}
self.tree_ = TREE_DICT[algorithm](
X,
metric=self.metric,
leaf_size=self.leaf_size,
sample_weight=sample_weight,
**kwargs,
)
return self
这段代码实现了KDE模型的拟合：它先选择合适的算法（KDTree或BallTree），根据带宽参数（若为"scott"或"silverman"则计算默认值）设置带宽，验证并处理样本权重，最后构建指定的树结构（KDTree或BallTree）用于后续密度查询。
源码路径：sklearn/neighbors/_kde.py - KernelDensity.score_samples（100-150行）
def score_samples(self, X):
"""Compute the log-likelihood of each sample under the model.
Parameters
----------
X : array-like of shape (n_samples, n_features)
An array of points to query. Last dimension should match dimension
of training data (n_features).
Returns
-------
density : ndarray of shape (n_samples,)
Log-likelihood of each sample in `X`. These are normalized to be
probability densities, so values will be low for high-dimensional
data.
"""
check_is_fitted(self)
# The returned density is normalized to the number of points.
# For it to be a probability, we must scale it. For this reason
# we'll also scale atol.
X = validate_data(self, X, order="C", dtype=np.float64, reset=False)
if self.tree_.sample_weight is None:
N = self.tree_.data.shape[0]
else:
N = self.tree_.sum_weight
atol_N = self.atol * N
log_density = self.tree_.kernel_density(
X,
h=self.bandwidth_,
kernel=self.kernel,
atol=atol_N,
rtol=self.rtol,
breadth_first=self.breadth_first,
return_log=True,
)
log_density -= np.log(N)
return log_density
这段代码实现了KDE的对数似然计算：它调用底层树结构的kernel_density方法计算每个查询点的密度估计（带宽、核函数等参数），然后通过除以总点数N（或权重和）将原始密度转换为对数概率密度。
26.6.2 局部离群因子的实现
源码路径：sklearn/neighbors/_lof.py - LocalOutlierFactor.fit（100-200行）
@_fit_context(
# LocalOutlierFactor.metric is not validated yet
prefer_skip_nested_validation=False
)
def fit(self, X, y=None):
"""Fit the local outlier factor detector from the training dataset.
Parameters
----------
X : {array-like, sparse matrix} of shape (n_samples, n_features) or \
(n_samples, n_samples) if metric='precomputed'
Training data.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : LocalOutlierFactor
The fitted local outlier factor detector.
"""
self._fit(X)
n_samples = self.n_samples_fit_
if self.n_neighbors > n_samples:
warnings.warn(
"n_neighbors (%s) is greater than the "
"total number of samples (%s). n_neighbors "
"will be set to (n_samples - 1) for estimation."
% (self.n_neighbors, n_samples)
)
self.n_neighbors_ = max(1, min(self.n_neighbors, n_samples - 1))
self._distances_fit_X_, _neighbors_indices_fit_X_ = self.kneighbors(
n_neighbors=self.n_neighbors_
)
if self._fit_X.dtype == np.float32:
self._distances_fit_X_ = self._distances_fit_X_.astype(
self._fit_X.dtype,
copy=False,
)
self._lrd = self._local_reachability_density(
self._distances_fit_X_, _neighbors_indices_fit_X_
)
# Compute lof score over training samples to define offset_:
lrd_ratios_array = (
self._lrd[_neighbors_indices_fit_X_] / self._lrd[:, np.newaxis]
)
self.negative_outlier_factor_ = -np.mean(lrd_ratios_array, axis=1)
if self.contamination == "auto":
# inliers score around -1 (the higher, the less abnormal).
self.offset_ = -1.5
else:
self.offset_ = np.percentile(
self.negative_outlier_factor_, 100.0 * self.contamination
)
# Verify if negative_outlier_factor_ values are within acceptable range.
# Novelty must also be false to detect outliers
if np.min(self.negative_outlier_factor_) < -1e7 and not self.novelty:
warnings.warn(
"Duplicate values are leading to incorrect results. "
"Increase the number of neighbors for more accurate results."
)
return self
这段代码实现了LOF模型的拟合：它先使用k近邻查询获取训练数据中每个点的最近邻及距离，然后计算局部可达密度（LRD），最后通过比较每个点的LRD与其邻居的LRD得到LOF得分，并根据污染率设置偏移量用于二分类判断。
源码路径：sklearn/neighbors/_lof.py - LocalOutlierFactor._local_reachability_density（50-100行）
def _local_reachability_density(self, distances_X, neighbors_indices):
"""The local reachability density (LRD)
The LRD of a sample is the inverse of the average reachability
distance of its k-nearest neighbors.
Parameters
----------
distances_X : ndarray of shape (n_queries, self.n_neighbors)
Distances to the neighbors (in the training samples `self._fit_X`)
of each query point to compute the LRD.
neighbors_indices : ndarray of shape (n_queries, self.n_neighbors)
Neighbors indices (of each query point) among training samples
self._fit_X.
Returns
-------
local_reachability_density : ndarray of shape (n_queries,)
The local reachability density of each sample.
"""
dist_k = self._distances_fit_X_[neighbors_indices, self.n_neighbors_ - 1]
reach_dist_array = np.maximum(distances_X, dist_k)
# 1e-10 to avoid `nan' when nb of duplicates > n_neighbors_:
return 1.0 / (np.mean(reach_dist_array, axis=1) + 1e-10)
这段代码实现了局部可达密度的计算：对于每个查询点，它获取其k个最近邻到训练数据的距离，计算到达距离（原始距离与邻居的k近邻距离的最大值），取平均值的倒数得到LRD。加入1e-10防止除以零。
26.7 邻域成分分析
邻域成分分析（NCA）是一种有监督的度量学习方法，它学习一个线性变换以最大化在变换后的空间中，正确分类的期望（基于随机近邻规则）。目标函数是每个点正确分类的概率之和，通过梯度下降优化（通常使用L-BFGS-B）来学习变换矩阵。
26.7.1 邻域成分分析的实现
源码路径：sklearn/neighbors/_nca.py - NeighborhoodComponentsAnalysis.fit（100-200行）
@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y):
"""Fit the model according to the given training data.
Parameters
----------
X : array-like of shape (n_samples, n_features)
The training samples.
y : array-like of shape (n_samples,)
The corresponding training labels.
Returns
-------
self : object
Fitted estimator.
"""
# Validate the inputs X and y, and converts y to numerical classes.
X, y = validate_data(self, X, y, ensure_min_samples=2)
check_classification_targets(y)
y = LabelEncoder().fit_transform(y)
# Check the preferred dimensionality of the projected space
if self.n_components is not None and self.n_components > X.shape[1]:
raise ValueError(
"The preferred dimensionality of the "
f"projected space `n_components` ({self.n_components}) cannot "
"be greater than the given data "
f"dimensionality ({X.shape[1]})!"
)
# If warm_start is enabled, check that the inputs are consistent
if (
self.warm_start
and hasattr(self, "components_")
and self.components_.shape[1] != X.shape[1]
):
raise ValueError(
f"The new inputs dimensionality ({X.shape[1]}) does not "
"match the input dimensionality of the "
f"previously learned transformation ({self.components_.shape[1]})."
)
# Check how the linear transformation should be initialized
init = self.init
if isinstance(init, np.ndarray):
init = check_array(init)
# Assert that init.shape[1] = X.shape[1]
if init.shape[1] != X.shape[1]:
raise ValueError(
f"The input dimensionality ({init.shape[1]}) of the given "
"linear transformation `init` must match the "
f"dimensionality of the given inputs `X` ({X.shape[1]})."
)
# Assert that init.shape[0] <= init.shape[1]
if init.shape[0] > init.shape[1]:
raise ValueError(
f"The output dimensionality ({init.shape[0]}) of the given "
"linear transformation `init` cannot be "
f"greater than its input dimensionality ({init.shape[1]})."
)
# Assert that self.n_components = init.shape[0]
if self.n_components is not None and self.n_components != init.shape[0]:
raise ValueError(
"The preferred dimensionality of the "
f"projected space `n_components` ({self.n_components}) does"
" not match the output dimensionality of "
"the given linear transformation "
f"`init` ({init.shape[0]})!"
)
# Initialize the random generator
self.random_state_ = check_random_state(self.random_state)
# Measure the total training time
t_train = time.time()
# Compute a mask that stays fixed during optimization:
same_class_mask = y[:, np.newaxis] == y[np.newaxis, :]
# (n_samples, n_samples)
# Initialize the transformation
transformation = np.ravel(self._initialize(X, y, init))
# Create a dictionary of parameters to be passed to the optimizer
disp = self.verbose - 2 if self.verbose > 1 else -1
optimizer_params = {
"method": "L-BFGS-B",
"fun": self._loss_grad_lbfgs,
"args": (X, same_class_mask, -1.0),
"jac": True,
"x0": transformation,
"tol": self.tol,
"options": dict(
maxiter=self.max_iter,
**_get_additional_lbfgs_options_dict("disp", disp),
),
"callback": self._callback,
}
# Call the optimizer
self.n_iter_ = 0
opt_result = minimize(**optimizer_params)
# Reshape the solution found by the optimizer
self.components_ = opt_result.x.reshape(-1, X.shape[1])
# Stop timer
t_train = time.time() - t_train
if self.verbose:
cls_name = self.__class__.__name__
# Warn the user if the algorithm did not converge
if not opt_result.success:
warn(
"[{}] NCA did not converge: {}".format(
cls_name, opt_result.message
),
ConvergenceWarning,
)
print("[{}] Training took {:8.2f}s.".format(cls_name, t_train))
return self
这段代码实现了NCA模型的拟合：它验证输入，处理标签编码，检查维度一致性，初始化变换矩阵（支持多种初始化方式如'pca'、'lda'等），计算同类别掩码，定义损失函数和梯度（基于成对距离的softmax），调用L-BFGS-B优化器学习变换，最后将结果 reshape 为变换矩阵。
源码路径：sklearn/neighbors/_nca.py - NeighborhoodComponentsAnalysis._loss_grad_lbfgs（100-150行）
def _loss_grad_lbfgs(self, transformation, X, same_class_mask, sign=1.0):
"""Compute the loss and the loss gradient w.r.t. `transformation`.
Parameters
----------
transformation : ndarray of shape (n_components * n_features,)
The raveled linear transformation on which to compute loss and
evaluate gradient.
X : ndarray of shape (n_samples, n_features)
The training samples.
same_class_mask : ndarray of shape (n_samples, n_samples)
A mask where `mask[i, j] == 1` if `X[i]` and `X[j]` belong
to the same class, and `0` otherwise.
Returns
-------
loss : float
The loss computed for the given transformation.
gradient : ndarray of shape (n_components * n_features,)
The new (flattened) gradient of the loss.
"""
if self.n_iter_ == 0:
self.n_iter_ += 1
if self.verbose:
header_fields = ["Iteration", "Objective Value", "Time(s)"]
header_fmt = "{:>10} {:>20} {:>10}"
header = header_fmt.format(*header_fields)
cls_name = self.__class__.__name__
print("[{}]".format(cls_name))
print(
"[{}] {}\n[{}] {}".format(
cls_name, header, cls_name, "-" * len(header)
)
)
t_funcall = time.time()
transformation = transformation.reshape(-1, X.shape[1])
X_embedded = np.dot(X, transformation.T) # (n_samples, n_components)
# Compute softmax distances
p_ij = pairwise_distances(X_embedded, squared=True)
np.fill_diagonal(p_ij, np.inf)
p_ij = softmax(-p_ij) # (n_samples, n_samples)
# Compute loss
masked_p_ij = p_ij * same_class_mask
p = np.sum(masked_p_ij, axis=1, keepdims=True) # (n_samples, 1)
loss = np.sum(p)
# Compute gradient of loss w.r.t. `transform`
weighted_p_ij = masked_p_ij - p_ij * p
weighted_p_ij_sym = weighted_p_ij + weighted_p_ij.T
np.fill_diagonal(weighted_p_ij_sym, -weighted_p_ij.sum(axis=0))
gradient = 2 * X_embedded.T.dot(weighted_p_ij_sym).dot(X)
# time complexity of the gradient: O(n_components x n_samples x (
# n_samples + n_features))
if self.verbose:
t_funcall = time.time() - t_funcall
values_fmt = "[{}] {:>10} {:>20.6e} {:>10.2f}"
print(
values_fmt.format(
self.__class__.__name__, self.n_iter_, loss, t_funcall
)
)
sys.stdout.flush()
return sign * loss, sign * gradient.ravel()
这段代码实现了NCA的损失函数和梯度计算：它将变换矩阵应用于数据得到嵌入空间，计算成对平方距离并应用softmax得到归一化相似度；损失是同类别点之间相似度的和；梯度通过链式法则计算，涉及嵌入空间的梯度和原始数据的乘法。
26.8 高斯混合模型 —— EM算法下的参数估计与协方差建模
高斯混合模型（Gaussian Mixture Model, GMM）通过期望最大化（EM）算法迭代估计模型参数，包含混合权重、均值向量和协方差矩阵。EM算法分为两步：E步计算每个样本属于各高斯分量的责任度（后验概率），M步根据责任度更新参数。协方差类型（full、tied、diag、spherical）直接影响参数数量和估计复杂度，AIC/BIC用于模型选择。为提升数值稳定性，似然计算通过协方差矩阵的Cholesky分解进行，避免直接求逆可能导致的不稳定性。
26.8.1 高斯混合模型的EM算法实现
源码路径：sklearn/mixture/_gaussian_mixture.py - GaussianMixture._m_step（360-400行）
def _m_step(self, X, log_resp, xp=None):
"""M step.
Parameters
----------
X : array-like of shape (n_samples, n_features)
log_resp : array-like of shape (n_samples, n_components)
Logarithm of the posterior probabilities (or responsibilities) of
the point of each sample in X.
"""
xp, _ = get_namespace(X, log_resp, xp=xp)
self.weights_, self.means_, self.covariances_ = _estimate_gaussian_parameters(
X, xp.exp(log_resp), self.reg_covar, self.covariance_type, xp=xp
)
self.weights_ /= xp.sum(self.weights_)
self.precisions_cholesky_ = _compute_precision_cholesky(
self.covariances_, self.covariance_type, xp=xp
)
这段代码实现了EM算法的M步：它接收对数责任度log_resp，先通过xp.exp还原为责任度，调用_estimate_gaussian_parameters根据责任度更新权重、均值和协方差；权重被归一化；最后根据更新的协方差计算精度矩阵的Cholesky分解，用于后续似然计算。
源码路径：sklearn/mixture/_gaussian_mixture.py - GaussianMixture._estimate_log_prob（410-440行）
def _estimate_log_prob(self, X, xp=None):
return _estimate_log_gaussian_prob(
X, self.means_, self.precisions_cholesky_, self.covariance_type, xp=xp
)
这段代码实现了E步中的似然计算：它调用_estimate_log_gaussian_prob，使用当前均值、精度矩阵的Cholesky分解和协方差类型，计算每个样本对各高斯分量的对数概率。使用Cholesky分解而非协方差矩阵直接求逆，是为了数值稳定性——特别是在协方差矩阵接近奇异时，Cholesky分解更稳健。
源码路径：sklearn/mixture/_gaussian_mixture.py - GaussianMixture._compute_lower_bound（480-520行）
def _compute_lower_bound(self, _, log_prob_norm):
return log_prob_norm
这段代码返回变分下界（在GMM中即对数似然）：对于标准高斯混合模型，变分下界正是数据的对数似然。log_prob_norm 是每个样本的平均对数似然乘以样本数，即总对数似然。EM算法通过迭代增加此下界来保证收敛。
26.8.2 协方差类型与参数数量
不同协方差类型决定了模型的表达力和参数复杂度：full 类型为每个成分分配独立协方差矩阵，参数数量为 O(KD²)；tied 类型所有成分共享一个协方差矩阵，参数数量为 O(D²)；diag 类型假设特征独立，协方差为对角矩阵，参数数量为 O(KD)；spherical 类型进一步假设各特征方差相等，参数数量仅为 O(K)。这直接影响AIC/BIC的计算，从而在模型选择中起关键作用。
26.9 贝叶斯高斯混合模型 —— 变分推断与自动模型选择
贝叶斯高斯混合模型（Bayesian Gaussian Mixture Model, BGMM）在标准GMM基础上引入先验分布，使用变分推断近似后验分布，并通过权重后验的稀疏性自动确定有效成分数。核心在于权重服从狄利克特过程或狄利克特分布先验，均值服从高斯先验，精度服从威沙特先验。变分下界（ELBO）用于监控收敛与模型选择，每次迭代后单调增加。
26.9.1 贝叶斯高斯混合的变分推断实现
源码路径：sklearn/mixture/_bayesian_mixture.py - BayesianGaussianMixture._estimate_log_weights（370-400行）
def _estimate_log_weights(self, xp=None):
if self.weight_concentration_prior_type == "dirichlet_process":
digamma_sum = digamma(
self.weight_concentration_[0] + self.weight_concentration_[1]
)
digamma_a = digamma(self.weight_concentration_[0])
digamma_b = digamma(self.weight_concentration_[1])
return (
digamma_a
- digamma_sum
+ np.hstack((0, np.cumsum(digamma_b - digamma_sum)[:-1]))
)
else:
# case Variational Gaussian mixture with dirichlet distribution
return digamma(self.weight_concentration_) - digamma(
np.sum(self.weight_concentration_)
)
这段代码计算权重参数的对数期望：当使用狄利克特过程先验时，权重后验是贝塔分布，其对数期望涉及digamma函数；当使用狄利克特分布先验时，权重后验仍是狄利克特分布。这一步对应变分推断中的M步，更新权重的后验超参数。
源码路径：sklearn/mixture/_bayesian_mixture.py - BayesianGaussianMixture._compute_lower_bound（450-490行）
def _compute_lower_bound(self, log_resp, log_prob_norm):
"""Estimate the lower bound of the model.
The lower bound on the likelihood (of the training data with respect to
the model) is used to detect the convergence and has to increase at
each iteration.
Parameters
----------
X : array-like of shape (n_samples, n_features)
log_resp : array, shape (n_samples, n_components)
Logarithm of the posterior probabilities (or responsibilities) of
the point of each sample in X.
log_prob_norm : float
Logarithm of the probability of each sample in X.
Returns
-------
lower_bound : float
"""
# Contrary to the original formula, we have done some simplification
# and removed all the constant terms.
(n_features,) = self.mean_prior_.shape
# We removed `.5 * n_features * np.log(self.degrees_of_freedom_)`
# because the precision matrix is normalized.
log_det_precisions_chol = _compute_log_det_cholesky(
self.precisions_cholesky_, self.covariance_type, n_features
) - 0.5 * n_features * np.log(self.degrees_of_freedom_)
if self.covariance_type == "tied":
log_wishart = self.n_components * np.float64(
_log_wishart_norm(
self.degrees_of_freedom_, log_det_precisions_chol, n_features
)
)
else:
log_wishart = np.sum(
_log_wishart_norm(
self.degrees_of_freedom_, log_det_precisions_chol, n_features
)
)
if self.weight_concentration_prior_type == "dirichlet_process":
log_norm_weight = -np.sum(
betaln(self.weight_concentration_[0], self.weight_concentration_[1])
)
else:
log_norm_weight = _log_dirichlet_norm(self.weight_concentration_)
return (
-np.sum(np.exp(log_resp) * log_resp)
- log_wishart
- log_norm_weight
- 0.5 * n_features * np.sum(np.log(self.mean_precision_))
)
这段代码实现了变分下界（ELBO）的计算：它由五部分组成：负责任度的熵项（-np.sum(np.exp(log_resp) * log_resp)）；威沙特先验对精度的贡献（log_wishart）；狄利克特先验对权重的贡献（log_norm_weight）；均值先验项（-0.5 * n_features * np.sum(np.log(self.mean_precision_))）；以及数据对数似然项（log_prob_norm 在调用者处已加入）。ELBO的单调增加是变分推断收敛的依据，同时其值可用于模型选择（如比较不同成分数的模型）。
26.9.2 自动模型选择机制
贝叶斯高斯混合通过权重后验的稀疏性实现自动确定成分数。狄利克特过程先验鼓励稀疏权重分布，使得多余成分的权重趋近于零；权重后验的期望直接用于计算有效成分数。这种机制无需依赖AIC/BIC等外部准则，内在地将模型复杂度与数据拟合度平衡。
26.10 协方差估计方法 —— 从经典到鲁棒的协方差建模
协方差估计是多维数据分析的基础，不同方法在鲁棒性、假设和适用场景上存在权衡。经验协方差是最大似然估计，但对离群点敏感；收缩协方差引入结构先验（如均值对角矩阵）以减少噪声；鲁棒协方差（如MCD）通过忽略极端点抵抗离群影响；最后，基于鲁棒协方差的椭圆包络可进行概率化的异常检测。
26.10.1 经验协方差的实现
源码路径：sklearn/covariance/_empirical_covariance.py - EmpiricalCovariance.fit（180-210行）
def fit(self, X, y=None):
"""Fit the maximum likelihood covariance estimator to X.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data, where `n_samples` is the number of samples and
`n_features` is the number of features.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
X = validate_data(self, X)
if self.assume_centered:
self.location_ = np.zeros(X.shape[1])
else:
self.location_ = X.mean(0)
covariance = empirical_covariance(X, assume_centered=self.assume_centered)
self._set_covariance(covariance)
return self
这段代码实现了最大似然协方差估计：它先验证数据（处理缺失值等），根据assume_centered决定是否减均值，调用empirical_covariance计算协方差（除以n而非n-1），最后通过_set_covariance存储协方差并可选计算精度矩阵。这是协方差估计的基础形式，但在高维小样本或存在离群点时会不稳定。
26.10.2 收缩协方差的实现
源码路径：sklearn/covariance/_shrunk_covariance.py - ShrunkCovariance.fit（180-210行）
def fit(self, X, y=None):
"""Fit the shrunk covariance model to X.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data, where `n_samples` is the number of samples
and `n_features` is the number of features.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
X = validate_data(self, X)
# Not calling the parent object to fit, to avoid a potential
# matrix inversion when setting the precision
if self.assume_centered:
self.location_ = np.zeros(X.shape[1])
else:
self.location_ = X.mean(0)
covariance = empirical_covariance(X, assume_centered=self.assume_centered)
covariance = shrunk_covariance(covariance, self.shrinkage)
self._set_covariance(covariance)
return self
这段代码实现了Ledoit-Wolf和OAS收缩估计的框架：它先计算经验协方差，然后应用收缩公式 (1-shrinkage) * empirical + shrinkage * structured，其中结构目标为均值对角矩阵。shrinkage 参数可手动指定（如ShrunkCovariance），或由LedoitWolf/OAS类自动估计（LedoitWolf最小化均方误差，OAS假设知道最优形式）。收缩旨在减少估计噪声，特别是在样本量小、维度高时改善条件数与稳定性。
26.10.3 鲁棒协方差的实现
源码路径：sklearn/covariance/_robust_covariance.py - MinCovDet.fit（190-240行）
def fit(self, X, y=None):
"""Fit a Minimum Covariance Determinant with the FastMCD algorithm.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Training data, where `n_samples` is the number of samples
and `n_features` is the number of features.
y : Ignored
Not used, present for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
X = validate_data(self, X, ensure_min_samples=2, estimator="MinCovDet")
random_state = check_random_state(self.random_state)
n_samples, n_features = X.shape
# check that the empirical covariance is full rank
if (linalg.svdvals(np.dot(X.T, X)) > 1e-8).sum() != n_features:
warnings.warn(
"The covariance matrix associated to your dataset is not full rank"
)
# compute and store raw estimates
raw_location, raw_covariance, raw_support, raw_dist = fast_mcd(
X,
support_fraction=self.support_fraction,
cov_computation_method=self._nonrobust_covariance,
random_state=random_state,
)
if self.assume_centered:
raw_location = np.zeros(n_features)
raw_covariance = self._nonrobust_covariance(
X[raw_support], assume_centered=True
)
# get precision matrix in an optimized way
precision = linalg.pinvh(raw_covariance)
raw_dist = np.sum(np.dot(X, precision) * X, 1)
self.raw_location_ = raw_location
self.raw_covariance_ = raw_covariance
self.raw_support_ = raw_support
self.location_ = raw_location
self.support_ = raw_support
self.dist_ = raw_dist
# obtain consistency at normal models
self.correct_covariance(X)
# re-weight estimator
self.reweight_covariance(X)
return self
这段代码实现了鲁棒协方差估计的核心流程：它使用FastMCD算法寻找使协方差行列式最小的h个点子集（默认h = floor((n + p + 1)/2)），分为寻找候选子集、汇总候选、在全数据上精炼三步。correct_covariance 应用亚马逊校正因子使估计在正态分布下一致；reweight_covariance 通过重加权步骤提高效率（利用更多点估计协方差同时保持鲁棒性）。鲁棒协方差的核心优势在于抗离群点——它通过忽略最极端的点来估计协方差，因此不易被少数异常样本主导。
26.10.4 椭圆包络异常检测的实现
源码路径：sklearn/covariance/_elliptic_envelope.py - EllipticEnvelope.decision_function（230-260行）
def decision_function(self, X):
"""Compute the decision function of the given observations.
Parameters
----------
X : array-like of shape (n_samples, n_features)
The data matrix.
Returns
-------
decision : ndarray of shape (n_samples,)
Decision function of the samples.
It is equal to the shifted Mahalanobis distances.
The threshold for being an outlier is 0, which ensures a
compatibility with other outlier detection algorithms.
"""
check_is_fitted(self)
negative_mahal_dist = self.score_samples(X)
return negative_mahal_dist - self.offset_
这段代码实现了基于鲁棒协方差的异常检测决策函数：它先调用score_samples计算负Mahalanobis距离（即 - (x - μ)^T Σ^{-1} (x - μ)），再减去偏移量offset_。偏移量由污染率决定，确保训练数据中恰好有预期比例的点被标记为异常（决策函数 < 0）。这一步将鲁棒协方差转化为概率化的异常得分：决策函数越小，点越异常；决策函数为0是异常/内点的分界线。
26.11 设计中的取舍
为什么在高斯混合模型中使用Cholesky分解而非直接求逆协方差矩阵来计算似然？
直接求逆协方差矩阵在矩阵接近奇异时数值极不稳定，可能导致非正定或极大值。Cholesky分解仅要求矩阵正定（通过正则化项reg_covar保证），且计算更稳健：对数似然中，log|Σ| 等价于 2 * sum(log(diag(L))) 其中 L 是Cholesky因子，避免了求逆步骤。在协方差估计中，这一技巧同样被用于精度矩阵（协方差逆）的稳定计算。
贝叶斯高斯混合模型中的狄利克特过程先验与狄利克特分布先验有何区别？
狄利克特分布先验是有限混合模型的共轭先验，固定成分数；狄利克特过程先验是无限混合模型的近似（ Stick-breaking 表示），允许模型自动确定有效成分数。在实现上，狄利克特过程先验使权重后验成为贝塔分布的序列，权重期望呈递减趋势，从而自然地驱动多余成分权重趋于零；而狄利克特分布先验则要求预先指定成分数，无法自动精简。因此，狄利克特过程先验是贝叶斯高斯混合实现自动模型选择的关键。
决策树中，如何平衡特征分裂的贪心策略与计算效率？
决策树的分裂采用贪心策略：在每个节点，遍历候选特征和可能的阈值，选择使不纯度下降最大的分裂。虽然这能局部最优，但计算代价高，尤其在特征众多时。为了提高效率，scikit-learn引入了max_features参数，在每个节点只随机选择一部分特征考虑，降低了从O(n_features)到O(sqrt(n_features))或O(log2(n_features))的复杂度。此外，对于稀疏数据，使用专门的分裂器避免遍历零值；对于有序特征，通过排序和桶划分进一步优化阈值搜索。这些权衡使得在保持较好分裂质量的同时，显著降低了训练时间。
最近邻算法中，KDTree与BallTree的选择取决于什么因素？
KDTree和BallTree的选择主要取决于数据的维度和稀疏性。KDTree在低维（通常小于20）且数据较为均匀分布时表现优异，因为其超矩形划分能高效地剪枝无关区域。随着维度增加，KDTree的性能因“维度灾难”而下降，此时BallTree由于其球形划分在高维下具有更好的自适应能力，能更好地处理非均匀或聚类数据。此外，BallTree支持更广泛的距离度量（包括非欧氏度量），而KDTree仅支持Minkowski族距离。在scikit-learn的自动算法选择中，当特征数小于15或邻居数小于样本数一半时优先考虑树方法；否则退回 brute force。对带权重的Minkowski距离，由于KDTree不支持，会强制使用BallTree。
核密度估计（KDE）中的带宽选择如何影响偏差-方差权衡？
KDE的带宽直接控制估计的平滑程度：小带宽导致估计过于不平滑（高方差，低偏差），对噪声敏感；大带宽导致过度平滑（低方差，高偏差），可能掩藏真实的多峰结构。scikit-learn提供了两种自动带宽选择方法：“scott”规则适用于正态分布，带宽与样本数和维度相关；“silverman”规则更保守，考虑了维度的影响。带宽过小会产生尖峰，过大会使所有特征趋向均匀。在实际应用中，带宽需根据数据特征和下游任务（如聚类或异常检测）进行调整，交叉验证是常用方法。
局部离群因子（LOF）与传统基于距离的异常检测（如基于均值和协方差的Mahalanobis距离）相比，在什么情况下优势更明显？
LOF的优势在于其局部适应性：它不依赖全局分布假设，而是基于局部密度来判断异常。在数据存在密度变化（如聚类密度不均匀）时，基于全局均值和协方差的方法会失效，因为它假设数据服从单一高斯分布；而LOF通过比较一点的密度与其邻居的密度，能够识别出在稀疏区域内的点，即使该点在全局上看起来并不极端。例如，在一个由紧密簇和孤点组成的数据集中，LOB能准确识别孤点为异常，而基于均值和协方差的方法可能因簇的影响而误判。此外，LOF对非球形和非高斯分布具有鲁棒性，因为它仅依赖密度排名而非参数假设。
邻域成分分析（NCA）相对于线性判别分析（LDA）在什么情况下更具优势？
NCA相较于LDA的优势在于：它直接优化分类性能（通过随机近邻规则的期望），而LDA优化的是类间方差与类内方差的比率（ Fisher准则），这两者不总是等价，尤其当类别分布非高斯或协方差不同时；NCA能处理多类问题而不受类别数限制（LDA受限于c-1个维度）；NCA学习的是线性变换，可降维到任意维度，而LDA的降维维度固定；此外，NCA可以通过核技巧推广到非线性变换，而基础LDA是线性的。然而，NCA的训练代价更高，因为它需要优化一个非凸目标函数，而LDA有封闭解。在样本量大、特征数高且类别分布复杂时，NCA往往能提供更好的特征变换用于后续分类。
26.12 动手练习
-
阅读
sklearn/tree/_classes.py和sklearn/tree/_splitter.pyx，理解以下实现：-
BaseDecisionTree._fit如何初始化树的构建过程，包括特征选择、样本权重处理和分裂器初始化 -
BaseDecisionTree._prune_tree如何通过代价复杂度剪枝简化树结构 -
Splitter.node_split_best和Splitter.node_split_random如何根据不同的分裂策略寻找最佳特征和阈值
回答问题：
-
决策树在构建过程中如何处理连续特征和离散特征的分裂？
-
不同的分裂策略（best vs random）在计算复杂度和分裂质量上有什么trade-off？
-
剪枝如何防止决策树的过拟合，以及代价复杂度参数如何影响结果？
-
-
阅读
sklearn/neighbors/_base.py和sklearn/neighbors/_kd_tree.pyx，理解以下实现：-
NeighborsBase._fit如何根据数据特征自动选择合适的索引算法（KDTree、BallTree或brute force） -
KNeighborsMixin.kneighbors如何处理查询数据是训练数据自身的情况，以避免将样本自身视为最近邻 -
KDTree.query如何使用最佳优先搜索和剪枝技术高效查询最近邻
回答问题：
-
KDTree和BallTree在什么情况下会被优先选择，什么情况下会退回到brute force？
-
最近邻查询中，如何处理平局情况（即多个点具有相同的距离）？
-
自动算法选择中，考虑了哪些因素来决定是否使用树结构还是暴力搜索？
-
-
阅读
sklearn/neighbors/_kde.py和sklearn/neighbors/_lof.py，理解以下实现：-
KernelDensity.fit如何根据带宽参数选择合适的树结构并进行参数初始化 -
KernelDensity.score_samples如何利用底层树结构计算对数概率密度 -
LocalOutlierFactor.fit如何计算局部可达密度并通过比较邻居密度来识别异常 -
LocalOutlierFactor._local_reachability_density如何避免在重复点情况下出现数值不稳定
回答问题：
-
KDE中的带宽参数如何影响估计的平滑程度，以及如何通过“scott”和“silverman”规则自动选择？
-
LOF的核心思想是什么？它如何区分全局异常和局部异常？
-
在LOF中，为什么需要加入一个小的常数（如1e-10）来避免除以零？
-
-
阅读
sklearn/neighbors/_nca.py，理解以下实现：-
NeighborhoodComponentsAnalysis._initialize如何根据不同的初始化策略（'pca'、'lda'、'random'等）初始化线性变换矩阵 -
NeighborhoodComponentsAnalysis._loss_grad_lbfgs如何根据软最大化的成对距离计算损失和梯度
回答问题：
-
NCA的目标函数是什么？它如何与随机近邻规则相关联？
-
NCA在优化过程中如何处理同类别掩码，以及这对梯度计算有什么影响？
-
NCA相较于传统的线性变换方法（如PCA或LDA）有什么优势和局限？
-
-
阅读
sklearn/mixture/_gaussian_mixture.py和sklearn/mixture/_bayesian_mixture.py，理解以下实现：-
GaussianMixture._estimate_gaussian_parameters如何在 M 步中估计均值和协方差 -
GaussianMixture._compute_precision_cholesky如何从协方差计算精度矩阵的 Cholesky 分解 -
GaussianMixture.bic和aic方法如何基于似然和参数数量进行模型选择
回答问题：
-
EM 算法在高斯混合模型中的收敛准则是什么？
-
不同协方差类型（full、tied、diag、spherical）在参数数量上的区别是什么？
-
为什么在计算似然时需要使用精度矩阵的 Cholesky 分解而非直接使用协方差？
-
-
阅读
sklearn/mixture/_bayesian_mixture.py，理解以下实现：-
BayesianGaussianMixture._estimate_weights如何根据先验和数据更新权重参数（狄利克特过程或狄利克特分布） -
BayesianGaussianMixture._estimate_precisions如何根据 Wishart 先验更新协方差参数 -
BayesianGaussianMixture._compute_lower_bound如何计算变分下界以监控收敛
回答问题：
-
贝叶斯高斯混合模型中的权重先验是如何区分狄利克特过程和狄利克特分布的？
-
在贝叶斯框架下，均值先验如何影响后验均值的估计？
-
下界在变分推断中的作用是什么？它如何被用于判断收敛和模型选择？
-
-
阅读
sklearn/covariance/_empirical_covariance.py、sklearn/covariance/_shrunk_covariance.py和sklearn/covariance/_robust_covariance.py，理解以下实现：-
EmpiricalCovariance.fit如何计算最大似然协方差估计 -
ledoit_wolf_shrinkage如何根据经验协方差估计最优收缩系数 -
fast_mcd如何通过随机子集迭代寻找最小协方差行列式的支持集
回答问题：
-
经验协方差估计在存在离群点时的主要限制是什么？
-
Ledoit-Wolf 收缩如何在保持无偏性的前提下改善条件数？
-
最小协方差行列式（MCD）估计器如何通过重加权步骤提高效率同时保持鲁棒性？
-
26.13 本章小结
本章我们一起学习了以下概念：
| 概念 | 解释 |
|------|------|
| BaseDecisionTree | 决策树的基类，定义了通用接口和属性，不应直接实例化 |
| DecisionTreeClassifier | 用于分类的决策树，支持gini、entropy和log_loss准则 |
| DecisionTreeRegressor | 用于回归的决策树，支持squared_error、absolute_error和poisson准则 |
| ExtraTreeClassifier | 极端随机化决策树分类器，在特征选择时引入随机性 |
| ExtraTreeRegressor | 极端随机化决策树回归器，在特征选择时引入随机性 |
| BaseDecisionTree._fit | 决策树的核心拟合方法，处理数据验证、特征选择和树构建 |
| BaseDecisionTree._prune_tree | 使用代价复杂度剪枝防止过拟合 |
| Splitter.node_split_best | 贪心地寻找使不纯度下降最大的特征和阈值 |
| Splitter.node_split_random | 随机采样特征和阈值以寻找良好分裂， trade-off 计算效率与质量 |
| Tree.compute_feature_importances | 基于不纯度降低计算特征重要性，越重要的特征对不纯度的贡献越大 |
| NeighborsBase._fit | 最近邻基类的拟合方法，根据数据特征自动选择KDTree、BallTree或brute force |
| KNeighborsMixin.kneighbors | 查询K个最近邻，处理查询数据是训练数据自身的情况以避免自环 |
| KDTree.query | 使用空间划分和最佳优先搜索高效查询最近邻，利用树结构进行剪枝 |
| KernelDensity.fit | 拟合核密度估计模型，选择合适的树结构并处理带宽参数 |
| KernelDensity.score_samples | 计算每个样本的对数概率密度，利用底层树进行高效查询 |
| LocalOutlierFactor.fit | 拟合局部离群因子模型，计算局部可达密度并基于邻居比较识别异常 |
| LocalOutlierFactor._local_reachability_density | 计算局部可达密度，通过加入小常数避免重复点导致的数值问题 |
| NeighborhoodComponentsAnalysis.fit | 拟合邻域成分分析模型，学习线性变换以最大化随机近邻分类期望 |
| NeighborhoodComponentsAnalysis._loss_grad_lbfgs | 计算NCA的损失函数和梯度，基于软最大化的成对距离和同类别掩码 |
| GaussianMixture | 通过EM算法估计权重、均值、协方差；支持full/tied/diag/spherical四种协方差类型；用AIC/BIC选择模型 |
| BayesianGaussianMixture | 变分推断估计后验；狄利克特过程先验自动确定有效成分数；权重、均值、精度各有先验 |
| EmpiricalCovariance | 最大似然协方差估计；原始样本协方差；假设数据均值为零可跳过中心化 |
| ShrunkCovariance/LedoitWolf/OAS | 收缩协方差估计；经验协方差与结构目标的凸组合；LedoitWolf自动求最优收缩系数 |
| MinCovDet | 鲁棒协方差估计；最小协方差行列式估计；FastMCD算法；离群点鲁棒性 |
| EllipticEnvelope | 基于鲁棒协方差的异常检测；决策函数为负Mahalanobis距离；由污染率决定阈值 |
感谢你读到了这里，恭喜你，你已经完成了 xxx
26.14 架构与数据流图
上述图分别展示模块依赖、调用时序、数据流和架构分层。
第 27 章 —— 特征工程与预处理 —— 数据进化的“变形金刚”
27.1 学习目标
-
难度：★★★☆☆（3/5）
-
预备知识：Python 基础、面向对象编程与 Markdown/代码阅读基础
-
掌握高斯混合模型（GMM）的 EM 算法迭代流程、协方差约束类型及 BIC/AIC 模型选择
-
理解贝叶斯高斯混合模型的变分推断、Dirichlet 过程先验与 Wishart 协方差先验的作用
-
深入稳健协方差估计（Ledoit-Wolf 收缩、图拉索、最小协方差行列式 MCD）的原理与实现
-
掌握 EllipticEnvelope 基于稳健协方差的异常检测机制
-
熟练掌握数据预处理与特征工具链的实现细节
-
理解 OneHotEncoder、OrdinalEncoder、LabelEncoder 等类别特征编码器的实现原理
-
掌握 KBinsDiscretizer、FunctionTransformer 等离散化与自定义变换工具的使用
-
理解 QuantileTransformer、PowerTransformer 等分布变换器的实现原理
-
掌握 ColumnTransformer、TransformedTargetRegressor 等组合工具的实现细节
-
理解文本特征提取工具（CountVectorizer、TfidfVectorizer、HashingVectorizer）的实现机制
27.2 生活类比
想象数据预处理和特征工程是一场“大数据炼金术”实验：原始数据就像从矿井里开出的粗矿石，经过一系列精细的处理步骤，才能成为适合机器学习模型“熔炼”的纯净金属。在这个过程中，特征缩放就像对不同金属做统一的纯度检测，确保它们在合金中按比例存在；类别编码则像是把各种金属元素的符号翻译成炼金术师统一理解的代码；离散化则像把连续的温度梯度分成若干火候档位；而特征组合工具 ColumnTransformer 就像是炼金工作台上的多工位流水线，能同时处理不同类型的原料；目标变换则是对最终产品进行后处理，使其更符合使用场景的需求；文本特征提取则是把古卷轴上的文字转化为可计算的数值向量，为后续的语义分析提供基础。整个流程就像一座智能工厂，每个预处理器都是一台专门的机床，它们组合在一起，才能把看似杂乱无章的原始数据，精密加工成机器学习算法所能“消化”的高质量特征。
27.3 源码地图
sklearn/mixture/_gaussian_mixture.py
├── GaussianMixture # 高斯混合模型
│ ├── fit() # EM 算法主循环
│ ├── _initialize() # 参数初始化 (k-means/k-means++)
│ ├── _e_step() # E 步：计算责任矩阵
│ ├── _m_step() # M 步：更新参数 (均值/协方差/权重)
│ ├── _compute_lower_bound() # 计算下界 (ELBO)
│ ├── predict() # 预测簇标签
│ ├── predict_proba() # 预测后验概率
│ ├── score() # 平均对数似然
│ ├── score_samples() # 样本对数概率密度
│ ├── sample() # 从模型采样
│ ├── bic() # 贝叶斯信息准则
│ ├── aic() # 赤池信息准则
├── BaseMixture # 混合模型基类
│ ├── _check_parameters() # 参数校验
│ ├── _check_is_fitted() # 拟合状态检查
│ ├── _print_verbose_msg_init_beg() # 打印初始化信息
sklearn/mixture/_bayesian_mixture.py
├── BayesianGaussianMixture # 贝叶斯高斯混合模型
│ ├── fit() # 变分推断主循环
│ ├── _initialize() # 参数初始化
│ ├── _estimate_log_weights() # 估计权重对数
│ ├── _estimate_log_prob() # 估计对数概率
│ ├── _e_step() # E 步
│ ├── _m_step() # M 步：更新变分参数
│ ├── _compute_lower_bound() # 计算变分下界
│ ├── _get_vb_prior() # 获取变分贝叶斯先验
sklearn/covariance/_empirical_covariance.py
├── empirical_covariance() # 经验协方差计算
sklearn/covariance/_shrunk_covariance.py
├── ledoit_wolf() # Ledoit-Wolf 收缩估计
├── ledoit_wolf_shrinkage() # 计算收缩系数
├── oas() # Oracle 收缩近似估计
sklearn/covariance/_graph_lasso.py
├── graphical_lasso() # 图拉索算法
├── GraphicalLasso # 图拉索估计器
│ ├── fit() # 坐标下降求解
├── GraphicalLassoCV # 交叉验证图拉索
│ ├── fit() # CV 选择正则化参数
sklearn/covariance/_robust_covariance.py
├── fast_mcd() # 快速 MCD 算法
├── MinCovDet # 最小协方差行列式估计器
│ ├── fit() # 拟合 MCD
│ ├── _c_step() # C-步迭代
│ ├── _get_support() # 获取内点支撑集
├── _c_step() # 核心 C-步函数
sklearn/covariance/_elliptic_envelope.py
├── EllipticEnvelope # 椭圆包络异常检测
│ ├── fit() # 拟合稳健协方差
│ ├── decision_function() # 决策函数 (马氏距离)
│ ├── score_samples() # 样本得分
│ ├── predict() # 预测异常/正常
sklearn/preprocessing/_data.py
├── StandardScaler # 标准化特征缩放器
│ ├── fit() # 拟合均值和方差
│ ├── transform() # 执行标准化转换
│ ├── inverse_transform() # 逆向标准化转换
│ ├── sklearn_tags # 返回标签信息
├── MinMaxScaler # 最小-最大特征缩放器
│ ├── fit() # 拟合最小值和最大值
│ ├── transform() # 执行缩放转换
│ ├── inverse_transform() # 逆向缩放转换
│ ├── sklearn_tags # 返回标签信息
├── RobustScaler # 鲁棒特征缩放器
│ ├── fit() # 拟合中位数和四分位距
│ ├── transform() # 执行鲁棒缩放转换
│ ├── inverse_transform() # 逆向鲁棒缩放转换
│ ├── sklearn_tags # 返回标签信息
├── MaxAbsScaler # 最大绝对值特征缩放器
│ ├── fit() # 拟合最大绝对值
│ ├── transform() # 执行缩放转换
│ ├── inverse_transform() # 逆向缩放转换
│ ├── sklearn_tags # 返回标签信息
├── Normalizer # 样本归一化器
│ ├── fit() # 参数验证（无状态转换器）
│ ├── transform() # 执行归一化转换
│ += sklearn_tags # 返回标签信息
├── QuantileTransformer # 分位数变换器
│ ├── fit() # 拟合分位数
│ ├── transform() # 执行分位数变换
│ ├── inverse_transform() # 逆向分位数变换
│ += sklearn_tags # 返回标签信息
├── PowerTransformer # 力量变换器
│ ├── fit() # 拟合lambda参数
│ ├── transform() # 执行力量变换
│ ├── inverse_transform() # 逆向力量变换
│ += sklearn_tags # 返回标签信息
├── scale # 特征标准化函数
├── minmax_scale # 最小-最大缩放函数
├── maxabs_scale # 最大绝对值缩放函数
├── robust_scale # 鲁棒缩放函函数
├── normalize # 样本归一化函数
├── binarize # 二值化函数
├── quantile_transform # 分位数变换函数
├── power_transform # 力量变换函数
sklearn/preprocessing/_encoders.py
├── OneHotEncoder # 独热编码器
│ ├── fit() # 拟合类别
│ ├── transform() # 执行独热编码
│ ├── inverse_transform() # 逆向独热编码
│ += get_feature_names_out() # 获取输出特征名
├── OrdinalEncoder # 序数编码器
│ ├── fit() # 拟合类别
│ ├── transform() # 执行序数编码
│ += inverse_transform() # 逆向序数编码
sklearn/preprocessing/_label.py
├── LabelBinarizer # 标签二值化器
│ ├── fit() # 拟合标签
│ ├── transform() # 执行二值化
│ += inverse_transform() # 逆向二值化
├── LabelEncoder # 标签编码器
│ ├── fit() # 拟合标签
│ += transform() # 执行编码
│ += inverse_transform() # 逆向编码
├── MultiLabelBinarizer # 多标签二值化器
│ ├── fit() # 拟合标签
│ += transform() # 执行多标签二值化
│ += inverse_transform() # 逆向多标签二值化
├── label_binarize() # 标签二值化函数
sklearn/preprocessing/_polynomial.py
├── PolynomialFeatures # 多项式特征生成器
│ ├── fit() # 计算输出特征数
│ += transform() # 生成多项式特征
│ += get_feature_names_out() # 获取输出特征名
├── SplineTransformer # 样条变换器
│ += fit() # 计算节点位置
│ += transform() # 生成样条基
│ += get_feature_names_out() # 获取输出特征名
sklearn/preprocessing/_discretization.py
├── KBinsDiscretizer # 离散化分箱器
│ += fit() # 计算分箱边界
│ += transform() # 执行离散化
│ += inverse_transform() # 逆向离散化
│ += get_feature_names_out() # 获取输出特征名
sklearn/preprocessing/_function_transformer.py
├── FunctionTransformer # 函数变换器
│ += fit() # 参数验证
│ += transform() # 执行函数变换
│ += inverse_transform() # 逆向函数变换
│ += get_feature_names_out() # 获取输出特征名
sklearn/preprocessing/_target_encoder.py
├── TargetEncoder # 目标编码器
│ += fit() # 拟合目标编码
│ += fit_transform() # 交叉拟合目标编码
│ += transform() # 执行目标编码
│ += inverse_transform() # 逆向目标编码
│ += get_feature_names_out() # 获取输出特征名
│ += get_metadata_routing() # 获取元数据路由
sklearn/preprocessing/_target_encoder_fast.pyx
├── _fit_encoding_fast() # 快速目标编码拟合
├── _fit_encoding_fast_auto_smooth() # 自动平滑快速目标编码拟合
sklearn/compose/_column_transformer.py
├── ColumnTransformer # 列变换器
│ += fit() # 拟合所有变换器
│ += fit_transform() # 拟合并转换
│ += transform() # 转换数据
│ += get_feature_names_out() # 获取输出特征名
│ += set_output() # 设置输出容器类型
├── make_column_selector() # 列选择器工厂
├── make_column_transformer() # 列变换器工厂
sklearn/compose/_target.py
├── TransformedTargetRegressor # 目标变换回归器
│ += fit() # 拟合目标变换回归器
│ += predict() # 预测并逆变换
│ += _fit_transformer() # 拟合目标变换器
│ += _get_regressor() # 获取回归器
│ += get_metadata_routing() # 获取元数据路由
sklearn/feature_extraction/text.py
├── CountVectorizer # 词频向量化器
│ += fit() # 构建词汇表
│ += fit_transform() # 拟合并转换
│ += transform() # 向量化文档
│ += inverse_transform() # 逆向向量化
│ += get_feature_names_out() # 获取特征名
├── TfidfVectorizer # TF-IDF向量化器
│ += fit() # 拟合词汇表和IDF
│ += fit_transform() # 拟合并转换
│ += transform() # TF-IDF向量化
├── HashingVectorizer # 哈希向量化器
│ += fit() # 参数验察
│ += transform() # 哈希变换
│ += fit_transform() # 拟合并转换
├── TfidfTransformer # TF-IDF变换器
│ += fit() # 计算IDF
│ += transform() # TF-IDF变换
├── strip_accents_unicode() # Unicode去重音
├── strip_accents_ascii() # ASCII去重音
├── strip_tags() # 去除HTML标签
sklearn/feature_extraction/_hash.py
├── FeatureHasher # 特征哈希器
│ += fit() # 参数验证
│ += transform() # 哈希变换
sklearn/feature_extraction/_hashing_fast.pyx
├── transform() # 核心哈希变换实现
sklearn/feature_extraction/_dict_vectorizer.py
├── DictVectorizer # 字典向量化器
│ += fit() # 构建特征映射
│ += fit_transform() # 拟合并转换
│ += transform() # 向量化字典
│ += inverse_transform() # 逆向向量化
│ += get_feature_names_out() # 获取特征名
sklearn/feature_extraction/image.py
├── extract_patches_2d() # 提取2D图像补丁
├── reconstruct_from_patches_2d() # 从补丁重建图像
├── img_to_graph() # 图像转图
├── grid_to_graph() # 网格转图
├── PatchExtractor # 补丁提取器
│ += fit() # 参数验证
│ += transform() # 提取补丁
sklearn/feature_extraction/_stop_words.py
├── ENGLISH_STOP_WORDS # 英文停用词集
sklearn/impute/_base.py
├── SimpleImputer # 简单插补器
│ += fit() # 计算插补统计量
│ += transform() # 执行插补
│ += inverse_transform() # 逆向插补
│ += get_feature_names_out() # 获取输出特征名
├── MissingIndicator # 缺失指示器
│ += fit() # 计算缺失模式
│ += transform() # 生成缺失指示矩阵
│ += fit_transform() # 拟合并转换
│ += get_feature_names_out() # 获取输出特征名
├── _BaseImputer # 插补器基类
sklearn/impute/_iterative.py
├── IterativeImputer # 迭代插补器
│ += fit() # 拟合迭代插补器
│ += fit_transform() # 交叉拟合插补
│ += transform() # 执行插补
│ += get_feature_names_out() # 获取输出特征名
│ += get_metadata_routing() # 获取元数据路由
sklearn/impute/_knn.py
├── KNNImputer # K近邻插补器
│ += fit() # 存储训练数据
│ += transform() # K近邻插补
│ += get_feature_names_out() # 获取输出特征名
sklearn/cross_decomposition/_pls.py
├── PLSRegression # PLS回归
│ += fit() # 拟合PLS回归
│ += transform() # 降维变换
│ += predict() # 预测目标
├── PLSCanonical # PLS典型相关分析
│ += fit() # 拟合PLS典型相关
│ += transform() # 变换X和Y
├── CCA # 典型相关分析
│ += fit() # 拟合CCA
│ += transform() # 变换X和Y
├── PLSSVD # PLS奇异值分解
│ += fit() # 拟合PLSSVD
│ += transform() # 变换X和Y
27.4 高斯混合模型 —— EM 算法的“概率聚类师”
高斯混合模型（Gaussian Mixture Model, GMM）通过期望最大化（Expectation-Maximization, EM）算法迭代拟合高斯分布的混合：在 E 步中，根据当前参数估计计算每个样本属于每个高斯分量的后验概率（称为责任）；在 M 步中，根据这些责任更新均值、协方差和混合权重。这一过程会重复进行，直到收敛。
GMM 支持四种协方差矩阵的约束类型，以适应不同的数据形状：
-
full：每个分量都有自己的完整协方差矩阵，能够捕捉任意椭球形状的簇；
-
tied：所有分量共享同一个协方差矩阵，适用于簇形状相似但位置不同的情况；
-
diag：每个分量的协方差矩阵是对角线矩阵，假设特征之间不相关；
-
spherical：每个分量的协方差矩阵是标量乘以单位矩阵，即所有特施具有相同的方差且互不相关。
为了防止协方差矩阵在迭代过程中变得奇异（特别是在样本较少时），GMM 引入了一个正则化参数 reg_covar，它会被加到每个协方差矩阵的对角线上，确保矩阵始终是正定的。此外，为了在数值上保持稳定，特别是在计算责任时，GMM 使用了 log-sum-exp 技巧来避免下溢。
为了自动选择最优的高斯分量数量，GMM 提供了基于贝叶斯信息准则（BIC）和赤池信息准则（AIC）的模型选择方法。这两个准则都在对数似然的基础上加上了一个惩罚项，以抑制过于复杂的模型（即分量过多）。BIC 的惩罚项更重，倾向于选择更简单的模型，而 AIC 则较为宽松。在实际应用中，我们通常会计算不同分量数下的 BIC 或 AIC 值，并选择使其最小的那个分量数作为最终模型的复杂度。
一旦模型训练完成，GMM 不仅能够给出硬聚类标签（通过 predict() 方法），还能输出每个样本属于各个分量的软聚类概率（通过 predict_proba() 方法），这在需要不确定性量化的场景中非常有用。此外，GMM 还支持从已训练的模型中生成新的样本（通过 sample() 方法），这使得它不仅可以用于聚类，还可以用于数据生成和密度估计。
27.4.1 核心类型定义：
源码路径：sklearn/mixture/_gaussian_mixture.py - GaussianMixture（1-58行）
// 类型即图纸，先理解数据结构
class GaussianMixture(BaseMixture):
"""Gaussian Mixture.
Representation of a Gaussian mixture model probability distribution.
This class allows to estimate the parameters of a Gaussian mixture
distribution.
.. versionadded:: 0.18
Parameters
----------
n_components : int, default=1
The number of mixture components.
covariance_type : {'full', 'tied', 'diag', 'spherical'}, default='full'
String describing the type of covariance parameters to use.
Must be one of:
- 'full': each component has its own general covariance matrix.
- 'tied': all components share the same general covariance matrix.
- 'diag': each component has its own diagonal covariance matrix.
- 'spherical': each component has its own single variance.
For an example of using `covariance_type`, refer to
:ref:`sphx_glr_auto_examples_mixture_plot_gmm_selection.py`.
tol : float, default=1e-3
The convergence threshold. EM iterations will stop when the
lower bound average gain is below this threshold.
reg_covar : float, default=1e-6
Non-negative regularization added to the diagonal of covariance.
Allows to assure that the covariance matrices are all positive.
max_iter : int, default=100
The number of EM iterations to perform.
n_init : int, default=1
The number of initializations to perform. The best results are kept.
init_params : {'kmeans', 'k-means++', 'random', 'random_from_data'}, \
default='kmeans'
The method used to initialize the weights, the means and the
precisions.
String must be one of:
- 'kmeans' : responsibilities are initialized using kmeans.
- 'k-means++' : use the k-means++ method to initialize.
- 'random' : responsibilities are initialized randomly.
- 'random_from_data' : initial means are randomly selected data points.
.. versionchanged:: v1.1
`init_params` now accepts 'random_from_data' and 'k-means++' as
initialization methods.
weights_init : array-like of shape (n_components, ), default=None
The user-provided initial weights.
If it is None, weights are initialized using the `init_params` method.
means_init : array-like of shape (n_components, n_features), default=None
The user-provided initial means,
If it is None, means are initialized using the `init_params` method.
precisions_init : array-like, default=None
The user-provided initial precisions (inverse of the covariance
matrices).
If it is None, precisions are initialized using the 'init_params'
method.
The shape depends on 'covariance_type'::
(n_components,) if 'spherical',
(n_features, n_features) if 'tied',
(n_components, n_features) if 'diag',
(n_components, n_features, n_features) if 'full'
random_state : int, RandomState instance or None, default=None
Controls the random seed given to the method chosen to initialize the
parameters (see `init_params`).
In addition, it controls the generation of random samples from the
fitted distribution (see the method `sample`).
Pass an int for reproducible output across multiple function calls.
See :term:`Glossary <random_state>`.
warm_start : bool, default=False
If 'warm_start' is True, the solution of the last fitting is used as
initialization for the next call of fit(). This can speed up
convergence when fit is called several times on similar problems.
In that case, 'n_init' is ignored and only a single initialization
occurs upon the first call.
See :term:`the Glossary <warm_start>`.
verbose : int, default=0
Enable verbose output. If 1 then it prints the current
initialization and each iteration step. If greater than 1 then
it prints also the log probability and the time needed
for each step.
verbose_interval : int, default=10
Number of iteration done before the next print.
Attributes
----------
weights_ : array-like of shape (n_components,)
The weights of each mixture components.
means_ : array-like of shape (n_components, n_features)
The mean of each mixture component.
covariances_ : array-like
The covariance of each mixture component.
The shape depends on `covariance_type`::
(n_components,) if 'spherical',
(n_features, n_features) if 'tied',
(n_components, n_features) if 'diag',
(n_components, n_features, n_features) if 'full'
For an example of using covariances, refer to
:ref:`sphx_glr_auto_examples_mixture_plot_gmm_covariances.py`.
precisions_ : array-like
The precision matrices for each component in the mixture. A precision
matrix is the inverse of a covariance matrix. A covariance matrix is
symmetric positive definite so the mixture of Gaussian can be
equivalently parameterized by the precision matrices. Storing the
precision matrices instead of the covariance matrices makes it more
efficient to compute the log-likelihood of new samples at test time.
The shape depends on `covariance_type`::
(n_components,) if 'spherical',
(n_features, n_features) if 'tied',
(n_components, n_features) if 'diag',
(n_components, n_features, n_features) if 'full'
precisions_cholesky_ : array-like
The Cholesky decomposition of the precision matrices of each mixture
component. A precision matrix is the inverse of a covariance matrix.
A covariance matrix is symmetric positive definite so the mixture of
Gaussian can be equivalently parameterized by the precision matrices.
Storing the precision matrices instead of the covariance matrices makes
it more efficient to compute the log-likelihood of new samples at test
time. The shape depends on `covariance_type`::
(n_components,) if 'spherical',
(n_features, n_features) if 'tied',
(n_components, n_features) if 'diag',
(n_components, n_features, n_features) if 'full'
converged_ : bool
True when convergence of the best fit of EM was reached, False otherwise.
n_iter_ : int
Number of step used by the best fit of EM to reach the convergence.
lower_bound_ : float
Lower bound value on the log-likelihood (of the training data with
respect to the model) of the best fit of EM.
lower_bounds_ : array-like of shape (`n_iter_`,)
The list of lower bound values on the log-likelihood from each
iteration of the best fit of EM.
n_features_in_ : int
Number of features seen during :term:`fit`.
.. versionadded:: 0.24
feature_names_in_ : ndarray of shape (`n_features_in_`,)
Names of features seen during :term:`fit`. Defined only when `X`
has feature names that are all strings.
.. versionadded:: 1.0
（这段代码定义了 GMM 的核心参数和属性，为 EM 算法的实现提供了模板。）
27.4.2 逐行解析关键函数：
源码路径：sklearn/mixture/_gaussian_mixture.py - GaussianMixture.fit()（59-110行）
// 逐行注释解释
def fit(self, X, y=None):
"""Estimate the model parameters with the EM algorithm.
The method fits the model to n_samples training data X using the
Expectation-Maximization (EM) algorithm.
Parameters
----------
X : array-like of shape (n_samples, n_features)
List of n_features-dimensional data points. Each row
corresponds to a single data point.
y : Ignored
Not used, present here for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
# Validate input data
X = validate_data(
self,
X,
order="C",
dtype=FLOAT_DTYPES,
copy=False,
ensure_all_finite="allow-nan",
)
# Check that the parameters are well defined
self._check_parameters(X)
# Initialize the parameters
self._initialize_parameters(X, self.random_state)
# Initialize the monitoring variables
self.converged_ = False
self.n_iter_ = 0
self.lower_bound_ = -np.inf
self.lower_bounds_ = np.array([])
# Store the parameters of the initialisation for warm start
weights_init, means_init, precisions_init = self._get_parameters()
# Do EM iterations
for i in range(self.max_iter):
# E step
log_prob_norm, log_resp = self._e_step(X)
# Check for convergence. We use the lower bound
# for the likelihood as convergence criteria, not the increase
# in the log likelihood.
self.lower_bound_ = log_prob_norm
self.lower_bounds_ = np.append(self.lower_bounds_, log_prob_norm)
if abs(log_prob_norm - self.lower_bound_) < self.tol:
self.converged_ = True
break
# M step
self._m_step(X, log_resp)
# Check for convergence. If we converged, we break
if self.converged_:
break
else:
# If we did not converge
if self.verbose:
print(
"Initialization %d did not converge. "
"Try different init parameters, "
"or increase max_iter, tol "
"or check for degenerate columns."
% self._init
)
# Store the parameters
self._set_parameters(self._get_parameters())
return self
（这段代码实现了 GMM 的 EM 算法主循环：先验证输入数据，再检查参数合法性，然后进行参数初始化；进入迭代循环，每次迭代先执行 E 步计算责任和对数似然，再执行 M 步更新模型参数；每次迭代后检查是否满足收敛条件（基于下界变化小于容忍度），如果满足则提前退出；如果达到最大迭代次数仍未收敛且开启了详细输出，则打印警告；最后将最终参数设置回模型中。）
源码路径：sklearn/mixture/_gaussian_mixture.py - GaussianMixture._e_step()（111-140行）
// 逐行注释解释
def _e_step(self, X):
"""E step.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Returns
-------
log_prob_norm : float
Mean of the log of the probability of each sample in X
log_resp : array-like of shape (n_samples, n_components)
Logarithm of the posterior probabilities (or responsibilities) of
the point of each sample in X.
"""
return (
self._estimate_log_prob(X),
self._estimate_log_prob(X) + self._estimate_log_weights(),
)
（这段代码实现了 E 步：调用 _estimate_log_prob 计算每个样本在每个高斯分量下的对数似然，调用 _estimate_log_weights 计算混合权重的对数，然后将两者相加得到未归一化的后验对数概率（责任）。返回值包括每个样本的平均对数似然（用于监控收敛）和责任矩阵的对数形式。）
源码路径：sklearn/mixture/_gaussian_mixture.py - GaussianMixture._m_step()（141-160行）
// 逐行注释解释
def _m_step(self, X, log_resp, xp=None):
"""M step.
Parameters
----------
X : array-like of shape (n_samples, n_features)
log_resp : array-like of shape (n_samples, n_components)
Logarithm of the posterior probabilities (or responsibilities) of
the point of each sample in X.
"""
xp, _ = get_namespace(X, log_resp, xp=xp)
self.weights_, self.means_, self.covariances_ = _estimate_gaussian_parameters(
X, xp.exp(log_resp), self.reg_covar, self.covariance_type, xp=xp
)
self.weights_ /= xp.sum(self.weights_)
self.precisions_cholesky_ = _compute_precision_cholesky(
self.covariances_, self.covariance_type, xp=xp
)
（这段代码实现了 M 步：首先将责任的对数形式通过指数函数还原为实际的责任值；调用 _estimate_gaussian_parameters 函数根据责任更新均值、协方差；然后对权重进行归一化（使其和为1）；最后根据更新后的协方差计算精度矩阵的 Cholesky 分解，用于后续的似然计算。）
27.4.3 完整数据流/流程图
27.5 贝叶斯高斯混合 —— 自动确定组件数的“变分魔法”
贝叶斯高斯混合模型（Bayesian Gaussian Mixture Model, BGM）通过变分推断（Variational Inference）近似后验分布，以自动推断有效的高斯分量数量。它采用了 Dirichlet 过程先验（对权重）和 Normal-Inverse-Wishart（NIW）先验（对均值和协方差），使得模型能够在不指定具体分量数的情况下，从数据中自动确定哪些分量是“被激活的”，哪些分量的权重趋近于零（相当于被删除）。
在变分框架中，真实的后验分布被一个易于处理的因子分布近似：q(Z, π, μ, Λ) ≈ q(Z)q(π)q(μ,Λ)，其中 Z 是潜在的分配变量，π 是混合权重，μ 和 Λ 分别是均值和精度矩阵。通过最大化证据下界（Evidence Lower Bound, ELBO），可以迭代更新这些变分分布的参数。
核心创新在于对权重的建模：使用 Dirichlet 过程先验，它可以看作是一个无限混合模型的有限截断（stick-breaking 构造）。该先验由一个浓度参数 weight_concentration_prior 控制：当该值小时，倾向于产生稀疏的权重分布（即只有少数分量有显著权重）；当该值大时，则倾向于使用更多的分量。通过在变分更新中估计后验的 Dirichlet 参数（即 weight_concentration_），我们可以推断出有多少个分量具有非 negligible 的权重。
此外，模型还提供了对均值和协方差的不确定性量化：通过保存均值分布和精度分布的变分近似（分别对应于高斯和 Wishart 分布），我们可以得到每个分量的均值和协方差的后验分布，而不仅仅是点估计。这使得 BGM 不仅能够进行聚类，还能提供贝叶斯意义上的不确定性估计。
27.5.1 核心类型定义：
源码路径：sklearn/mixture/_bayesian_mixture.py - BayesianGaussianMixture（1-60行）
// 类型即图纸，先理解数据结构
class BayesianGaussianMixture(BaseMixture):
"""Variational Bayesian estimation of a Gaussian mixture.
This class allows to infer an approximate posterior distribution over the
parameters of a Gaussian mixture distribution. The effective number of
components can be inferred from the data.
This class implements two types of prior for the weights distribution: a
finite mixture model with Dirichlet distribution and an infinite mixture
model with the Dirichlet Process. In practice Dirichlet Process inference
algorithm is approximated and uses a truncated distribution with a fixed
maximum number of components (called the Stick-breaking representation).
The number of components actually used almost always depends on the data.
.. versionadded:: 0.18
Read more in the :ref:`User Guide <bgmm>`.
Parameters
----------
n_components : int, default=1
The number of mixture components. Depending on the data and the value
of the `weight_concentration_prior` the model can decide to not use
all the components by setting some component `weights_` to values very
close to zero. The number of effective components is therefore smaller
than n_components.
covariance_type : {'full', 'tied', 'diag', 'spherical'}, default='full'
String describing the type of covariance parameters to use.
Must be one of:
- 'full' (each component has its own general covariance matrix),
- 'tied' (all components share the same general covariance matrix),
- 'diag' (each component has its own diagonal covariance matrix),
- 'spherical' (each component has its own single variance).
tol : float, default=1e-3
The convergence threshold. EM iterations will stop when the
lower bound average gain on the likelihood (of the training data with
respect to the model) is below this threshold.
reg_covar : float, default=1e-6
Non-negative regularization added to the diagonal of covariance.
Allows to assure that the covariance matrices are all positive.
max_iter : int, default=100
The number of EM iterations to perform.
n_init : int, default=1
The number of initializations to perform. The result with the highest
lower bound value on the likelihood is kept.
init_params : {'kmeans', 'k-means++', 'random', 'random_from_data'}, \
default='kmeans'
The method used to initialize the weights, the means and the
covariances. String must be one of:
- 'kmeans': responsibilities are initialized using kmeans.
- 'k-means++': use the k-means++ method to initialize.
- 'random': responsibilities are initialized randomly.
- 'random_from_data': initial means are randomly selected data points.
.. versionchanged:: v1.1
`init_params` now accepts 'random_from_data' and 'k-means++' as
initialization methods.
weight_concentration_prior_type : {'dirichlet_process', 'dirichlet_distribution'}, \
default='dirichlet_process'
String describing the type of the weight concentration prior.
weight_concentration_prior : float or None, default=None
The dirichlet concentration of each component on the weight
distribution (Dirichlet). This is commonly called gamma in the
literature. The higher concentration puts more mass in
the center and will lead to more components being active, while a lower
concentration parameter will lead to more mass at the edge of the
mixture weights simplex. The value of the parameter must be greater
than 0. If it is None, it's set to ``1. / n_components``.
mean_precision_prior : float or None, default=None
The precision prior on the mean distribution (Gaussian).
Controls the extent of where means can be placed. Larger
values concentrate the cluster means around `mean_prior`.
The value of the parameter must be greater than 0.
If it is None, it is set to 1.
mean_prior : array-like, shape (n_features,), default=None
The prior on the mean distribution (Gaussian).
If it is None, it is set to the mean of X.
degrees_of_freedom_prior : float or None, default=None
The prior of the number of degrees of freedom on the covariance
distributions (Wishart). If it is None, it's set to `n_features`.
covariance_prior : float or array-like, default=None
The prior on the covariance distribution (Wishart).
If it is None, the emiprical covariance prior is initialized using the
covariance of X. The shape depends on `covariance_type`::
(n_features, n_features) if 'full',
(n_features, n_features) if 'tied',
(n_features) if 'diag',
float if 'spherical'
random_state : int, RandomState instance or None, default=None
Controls the random seed given to the method chosen to initialize the
parameters (see `init_params`).
In addition, it controls the generation of random samples from the
fitted distribution (see the method `sample`).
Pass an int for reproducible output across multiple function calls.
See :term:`Glossary <random_state>`.
warm_start : bool, default=False
If 'warm_start' is True, the solution of the last fitting is used as
initialization for the next call of fit(). This can speed up
convergence when fit is called several times on similar problems.
See :term:`the Glossary <warm_start>`.
verbose : int, default=0
Enable verbose output. If 1 then it prints the current
initialization and each iteration step. If greater than 1 then
it prints also the log probability and the time needed
for each step.
verbose_interval : int, default=10
Number of iteration done before the next print.
Attributes
----------
weights_ : array-like of shape (n_components,)
The weights of each mixture components.
means_ : array-like of shape (n_components, n_features)
The mean of each mixture component.
covariances_ : array-like
The covariance of each mixture component.
The shape depends on `covariance_type`::
(n_components,) if 'spherical',
(n_features, n_features) if 'tied',
(n_components, n_features) if 'diag',
(n_components, n_features, n_features) if 'full'
precisions_ : array-like
The precision matrices for each component in the mixture. A precision
matrix is the inverse of a covariance matrix. A covariance matrix is
symmetric positive definite so the mixture of Gaussian can be
equivalently parameterized by the precision matrices. Storing the
precision matrices instead of the covariance matrices makes it more
efficient to compute the log-likelihood of new samples at test time.
The shape depends on ``covariance_type``::
(n_components,) if 'spherical',
(n_features, n_features) if 'tied',
(n_components, n_features) if 'diag',
(n_components, n_features, n_features) if 'full'
precisions_cholesky_ : array-like
The Cholesky decomposition of the precision matrices of each mixture
component. A precision matrix is the inverse of a covariance matrix.
A covariance matrix is symmetric positive definite so the mixture of
Gaussian can be equivalently parameterized by the precision matrices.
Storing the precision matrices instead of the covariance matrices makes
it more efficient to compute the log-likelihood of new samples at test
time. The shape depends on ``covariance_type``::
(n_components,) if 'spherical',
(n_features, n_features) if 'tied',
(n_components, n_features) if 'diag',
(n_components, n_features, n_features) if 'full'
converged_ : bool
True when convergence of the best fit of inference was reached, False otherwise.
n_iter_ : int
Number of step used by the best fit of inference to reach the
convergence.
lower_bound_ : float
Lower bound value on the model evidence (of the training data) of the
best fit of inference.
lower_bounds_ : array-like of shape (`n_iter_`,)
The list of lower bound values on the model evidence from each iteration
of the best fit of inference.
weight_concentration_prior_ : tuple or float
The dirichlet concentration of each component on the weight
distribution (Dirichlet). The type depends on
``weight_concentration_prior_type``::
(float, float) if 'dirichlet_process' (Beta parameters),
float if 'dirichlet_distribution' (Dirichlet parameters).
The higher concentration puts more mass in
the center and will lead to more components being active, while a lower
concentration parameter will lead to more mass at the edge of the
simplex.
weight_concentration_ : array-like of shape (n_components,)
The dirichlet concentration of each component on the weight
distribution (Dirichlet).
mean_precision_prior_ : float
The precision prior on the mean distribution (Gaussian).
Controls the extent of where means can be placed.
Larger values concentrate the cluster means around `mean_prior`.
If mean_precision_prior is set to None, `mean_precision_prior_` is set
to 1.
mean_precision_ : array-like of shape (n_components,)
The precision of each components on the mean distribution (Gaussian).
mean_prior_ : array-like of shape (n_features,)
The prior on the mean distribution (Gaussian).
degrees_of_freedom_prior_ : float
The prior of the number of degrees of freedom on the covariance
distributions (Wishart).
degrees_of_freedom_ : array-like of shape (n_components,)
The number of degrees of freedom of each components in the model.
covariance_prior_ : float or array-like
The prior on the covariance distribution (Wishart).
The shape depends on `covariance_type`::
(n_features, n_features) if 'full',
(n_features, n_features) if 'tied',
(n_features) if 'diag',
float if 'spherical'
n_features_in_ : int
Number of features seen during :term:`fit`.
.. versionadded:: 0.24
feature_names_in_ : ndarray of shape (`n_features_in_`,)
Names of features seen during :term:`fit`. Defined only when `X`
has feature names that are all strings.
.. versionadded:: 1.0
（这段代码定义了贝叶斯高斯混合模型的核心参数和属性，展示了它如何通过先验分布来建模不确定性，并如何通过变分推断自动确定组件数。）
27.5.2 逐行解析关键函数：
源码路径：sklearn/mixture/_bayesian_mixture.py - BayesianGaussianMixture.fit()（61-120行）
// 逐行注释解释
def fit(self, X, y=None):
"""Estimate the model parameters with the variational inference algorithm.
The method fits the model to n_samples training data X using the
variational inference algorithm.
Parameters
----------
X : array-like of shape (n_samples, n_features)
List of n_features-dimensional data points. Each row
corresponds to a single data point.
y : Ignored
Not used, present here for API consistency by convention.
Returns
-------
self : object
Returns the instance itself.
"""
# Validate input data
X = validate_data(
self,
X,
order="C",
dtype=FLOAT_DTYPES,
copy=False,
ensure_all_finite="allow-nan",
)
# Check that the parameters are well defined
self._check_parameters(X)
# Initialize the parameters
self._initialize_parameters(X, self.random_state)
# Initialize the monitoring variables
self.converged_ = False
self.n_iter_ = 0
self.lower_bound_ = -np.inf
self.lower_bounds_ = np.array([])
# Store the parameters of the initialisation for warm start
weights_init, means_init, precisions_init = self._get_parameters()
# Do variational inference iterations
for i in range(self.max_iter):
# E step
log_prob_norm, log_resp = self._e_step(X)
# Check for convergence. We use the lower bound
# for the likelihood as convergence criteria, not the increase
# in the log likelihood.
self.lower_bound_ = log_prob_norm
self.lower_bounds_ = np.append(self.lower_bounds_, log_prob_norm)
if abs(log_prob_norm - self.lower_bound_) < self.tol:
self.converged_ = True
break
# M step
self._m_step(X, log_resp)
# Check for convergence. If we converged, we break
if self.converged_:
break
else:
# If we did not converge
if self.verbose:
print(
"Initialization %d did not converge. "
"Try different init parameters, "
"or increase max_iter, tol "
"or check for degenerate columns."
% self._init
)
# Store the parameters
self._set_parameters(self._get_parameters())
return self
（这段代码实现了贝叶斯 GMM 的变分推断主循环，结构与标准 GMM 类似，但 E 步和 M 步中调用的函数不同，因为它们需要处理变分参数而不是标准参数。）
源码路径：sklearn/mixture/_bayesian_mixture.py - BayesianGaussianMixture._m_step()（150-200行）
// 逐行注释解释
def _m_step(self, X, log_resp, xp=None):
"""M step.
Parameters
----------
X : array-like of shape (n_samples, n_features)
log_resp : array-like of shape (n_samples, n_components)
Logarithm of the posterior probabilities (or responsibilities) of
the point of each sample in X.
"""
n_samples, _ = X.shape
nk, xk, sk = _estimate_gaussian_parameters(
X, np.exp(log_resp), self.reg_covar, self.covariance_type
)
self._estimate_weights(nk)
self._estimate_means(nk, xk)
self._estimate_precisions(nk, xk, sk)
（这段代码实现了贝叶斯 GMM 的 M 步：首先根据责任（通过指数还原 log_resp）估计高斯参数（均值、第二阶矩阵和样本数），然后依次更新权重的 Dirichlet 参数、均值的高斯参数和协方差的 Wishart 参数。每个更新都封装在各自的辅助函数中，使得主循环保持简洁。）
源码路径：sklearn/mixture/_bayesian_mixture.py - BayesianGaussianMixture._compute_lower_bound()（201-250行）
// 逐行注释解释
def _compute_lower_bound(self, log_resp, log_prob_norm):
"""Estimate the lower bound of the model.
The lower bound on the likelihood (of the training data with respect to
the model) is used to detect the convergence and has to increase at
each iteration.
Parameters
----------
log_resp : array, shape (n_samples, n_components)
Logarithm of the posterior probabilities (or responsibilities) of
the point of each sample in X.
log_prob_norm : float
Logarithm of the probability of each sample in X.
Returns
-------
lower_bound : float
"""
# Contrary to the original formula, we have done some simplification
# and removed all the constant terms.
(n_features,) = self.mean_prior_.shape
# We removed `.5 * n_features * np.log(self.degrees_of_freedom_)`
# because the precision matrix is normalized.
log_det_precisions_chol = _compute_log_det_cholesky(
self.precisions_cholesky_, self.covariance_type, n_features
) - 0.5 * n_features * np.log(self.degrees_of_freedom_)
if self.covariance_type == "tied":
log_wishart = self.n_components * np.float64(
_log_wishart_norm(
self.degrees_of_freedom_, log_det_precisions_chol, n_features
)
)
else:
log_wishart = np.sum(
_log_wishart_norm(
self.degrees_of_freedom_, log_det_precisions_chol, n_features
)
)
if self.weight_concentration_prior_type == "dirichlet_process":
log_norm_weight = -np.sum(
betaln(self.weight_concentration_[0], self.weight_concentration_[1])
)
else:
log_norm_weight = _log_dirichlet_norm(self.weight_concentration_)
return (
-np.sum(np.exp(log_resp) * log_resp)
- log_wishart
- log_norm_weight
- 0.5 * n_features * np.sum(np.log(self.mean_precision_))
)
（这段代码实现了变分下界（ELBO）的计算，它是变分推断的核心目标函数。通过在每次迭代后计算这个下界并检查其变化是否小于容忍度，我们可以判断算法是否收敛。下界由几部分组成：负的交叉熵项（-∑ exp(log_resp) * log_resp）、Wishart 先验的对数归一化项、Dirichlet 先验的对数归一化项，以及均值先验的对数精度项。通过最大化这个下界，我们实际上是在最小化 KL 散度，从而使变分分布逼近真实后验。）
27.5.3 完整数据流/流程图
27.6 协方差估计与异常检测 —— 多维数据的“散布测量”
在多维数据分析中，协方差矩阵描述了特征之间的线性关系，是许多算法（如 PCA、LDA、高斯混合模型）的基础。然而，经验协方差矩阵（即样本协方差）在样本量小时会不稳健，并且对离群点极其敏感。为了解决这一问题，scikit-learn 提供了一系列稳健和收缩式的协方差估计方法。
经验协方差是最大似然估计，但在高维或小样本情况下，它容易过拟合噪声。Ledoit-Wolf 收缩方法通过将经验协方差向单位矩阵（按特征方差均值缩放）进行线性收缩来改善这一点，这种收缩在均方误差意义下是最优的（在某些前提下）。Oracle 收缩近似（OAS）进一步假设数据服从高斯分布，从而在高斯前提下实现了理论上的最优收缩。
图拉索（Graphical Lasso）则采用了另一种思路：它不直接估计协方差，而是估计其逆矩阵（即精度矩阵），并在精度矩阵上加入 L1 正则项，以促进稀疏性。这使得我们不仅能得到一个稳健的协方差估计，还能通过精度矩阵中的零元素推断出哪些特征是条件独立的，从而学习一个稀疏的图结构——这在基因网络、金融系统等领域尤为重要。
最小协方差行列式（Minimum Covariance Determinant, MCD）是一种高破坏点（高达 50%）的稳健估计器：它通过寻找 h 个（h ≈ n/2）点的子集，使得这个子集的协方差行列式最小，从而有效地剔除离群点的影响。快速 MCD（FastMCD）算法通过在随机子集上进行初步搜索，然后通过 C-步迭代（不断用当前估计重新加权样本并重新估计）快速收敛到近似解。
基于稳健协方差（通常默认使用 MCD 满足 50% 破坏点），EllipticEnvelope 构建了一个异常检测器：它假设服从高斯分布的内生数据落在以稳健均值为中心、以稳健协方差为形状的椭圆内部，而外点则被视为异常。通过马氏距离（该距离考虑了特征间的相关性），我们可以将每个样本映射到一个标量得分；根据假设的污染比例（如 10%），我们设定一个阈值，使得恰好有该比例的样本得分超过阈值而被标记为异常。
27.6.1 核心类型定义：
源码路径：sklearn/covariance/_shrunk_covariance.py - ledoit_wolf（50-150行）
// 类型即图纸，先理解数据结构
def ledoit_wolf(X, *, assume_centered=False, block_size=1000):
"""Estimate the shrunk Ledoit-Wolf covariance matrix.
Read more in the :ref:`User Guide <shrunk_covariance>`.
Parameters
----------
X : array-like of shape (n_samples, n_features)
Data from which to compute the covariance estimate.
assume_centered : bool, default=False
If True, data will not be centered before computation.
Useful to work with data whose mean is significantly equal to
zero but is not exactly zero.
If False, data will be centered before computation.
block_size : int, default=1000
Size of blocks into which the covariance matrix will be split.
This is purely a memory optimization and does not affect results.
Returns
-------
shrunk_cov : ndarray of shape (n_features, n_features)
Shrunk covariance.
shrinkage : float
Coefficient in the convex combination used for the computation
of the shrunk estimate.
Notes
-----
The regularized (shrunk) covariance is:
(1 - shrinkage) * cov + shrinkage * mu * np.identity(n_features)
where mu = trace(cov) / n_features
Examples
--------
>>> import numpy as np
>>> from sklearn.covariance import empirical_covariance, ledoit_wolf
>>> real_cov = np.array([[.4, .2], [.2, .8]])
>>> rng = np.random.RandomState(0)
>>> X = rng.multivariate_normal(mean=[0, 0], cov=real_cov, size=50)
>>> covariance, shrinkage = ledoit_wolf(X)
>>> covariance
array([[0.44, 0.16],
[0.16, 0.80]])
>>> shrinkage
np.float64(0.23)
"""
estimator = LedoitWolf(
assume_centered=assume_centered,
block_size=block_size,
store_precision=False,
).fit(X)
return estimator.covariance_, estimator.shrinkage_
浙公网安备 33010602011771号