在统计学和机器学习领域，概率模型是理解和处理数据的基础。这些模型能够描述数据的生成过程，并且可以用于模型选择和参数估计。特别是，当新数据的似然性被用于模型选择和协方差估计时，可以通过比较不同模型的性能来选择最合适的模型。本文将重点比较主成分分析（PCA）和因子分析（FA）在处理低秩数据时的表现，这些数据可能被同质噪声（所有特征的噪声方差相同）或异质噪声（不同特征的噪声方差不同）所污染。
在第一步中，通过交叉验证比较PCA和FA在低秩数据上的表现。结果表明，在同质噪声条件下，PCA和FA都能够成功恢复低秩子空间的大小，且PCA的似然性高于FA。然而，当存在异质噪声时，PCA会失败并高估秩的大小。在适当的情况下（选择合适的成分数量），低秩模型的保留数据更可能被模型所接受，而不是收缩模型。
接下来，将模型的似然性与收缩协方差估计器得到的似然性进行比较。在Thomas P. Minka的论文《Automatic Choice of Dimensionality for PCA》中，也提到了自动估计降维的方法。通过实验，观察到在同质噪声下，PCA和FA都能成功地恢复低秩子空间的大小，其中PCA的似然性高于FA。但在异质噪声存在的情况下，PCA会失败并高估秩的大小。在适当条件下（选择合适的成分数量），保留数据更可能被低秩模型所接受，而不是收缩模型。
为了创建数据，使用numpy库生成随机数，并使用scipy库的线性代数模块进行奇异值分解。设定样本数量、特征数量和秩的大小，并使用随机数生成器来创建数据矩阵。接着，向数据中添加同质噪声和异质噪声，以模拟不同的数据污染情况。
import numpy as np
from scipy import linalg
n_samples, n_features, rank = 500, 25, 5
sigma = 1.0
rng = np.random.RandomState(42)
U, _, _ = linalg.svd(rng.randn(n_features, n_features))
X = np.dot(rng.randn(n_samples, rank), U[:, :rank].T)
# 添加同质噪声
X_homo = X + sigma * rng.randn(n_samples, n_features)
# 添加异质噪声
sigmas = sigma * rng.rand(n_features) + sigma / 2.0
X_hetero = X + rng.randn(n_samples, n_features) * sigmas
接下来，使用matplotlib库进行可视化，并从sklearn库中导入所需的模型和交叉验证工具。定义了一个函数来计算不同模型的分数，并使用网格搜索交叉验证来找到最佳的收缩协方差估计器。
import matplotlib.pyplot as plt
from sklearn.covariance import LedoitWolf, ShrunkCovariance
from sklearn.decomposition import PCA, FactorAnalysis
from sklearn.model_selection import GridSearchCV, cross_val_score
n_components = np.arange(0, n_features, 5)
def compute_scores(X):
pca = PCA(svd_solver="full")
fa = FactorAnalysis()
pca_scores, fa_scores = [], []
for n in n_components:
pca.n_components = n
fa.n_components = n
pca_scores.append(np.mean(cross_val_score(pca, X)))
fa_scores.append(np.mean(cross_val_score(fa, X)))
return pca_scores, fa_scores
def shrunk_cov_score(X):
shrinkages = np.logspace(-2, 0, 30)
cv = GridSearchCV(ShrunkCovariance(), {"shrinkage": shrinkages})
return np.mean(cross_val_score(cv.fit(X).best_estimator_, X))
def lw_score(X):
return np.mean(cross_val_score(LedoitWolf(), X))
for X, title in [(X_homo, "同质噪声"), (X_hetero, "异质噪声")]:
pca_scores, fa_scores = compute_scores(X)
n_components_pca = n_components[np.argmax(pca_scores)]
n_components_fa = n_components[np.argmax(fa_scores)]
pca = PCA(svd_solver="full", n_components="mle")
pca.fit(X)
n_components_pca_mle = pca.n_components_
print("最佳n_components由PCA CV确定 = %d" % n_components_pca)
print("最佳n_components由因子分析CV确定 = %d" % n_components_fa)
print("最佳n_components由PCA MLE确定 = %d" % n_components_pca_mle)
plt.figure()
plt.plot(n_components, pca_scores, "b", label="PCA分数")
plt.plot(n_components, fa_scores, "r", label="FA分数")
plt.axvline(rank, color="g", label="真实值:%d" % rank, linestyle="-")
plt.axvline(n_components_pca, color="b", label="PCA CV:%d" % n_components_pca, linestyle="--")
plt.axvline(n_components_fa, color="r", label="因子分析CV:%d" % n_components_fa, linestyle="--")
plt.axvline(n_components_pca_mle, color="k", label="PCA MLE:%d" % n_components_pca_mle, linestyle="--")
plt.axhline(shrunk_cov_score(X), color="violet", label="收缩协方差MLE", linestyle="-.")
plt.axhline(lw_score(X), color="orange", label="LedoitWolf MLE", linestyle="-.")
plt.xlabel("成分数量")
plt.ylabel("CV分数")
plt.legend(loc="lower right")
plt.title(title)
plt.show()
通过上述代码，可以比较不同模型在不同噪声条件下的表现。在同质噪声条件下，PCA和FA都能成功恢复低秩子空间的大小，且PCA的似然性高于FA。然而，在异质噪声条件下，PCA会失败并高估秩的大小。在适当条件下（选择合适的成分数量），保留数据更可能被低秩模型所接受，而不是收缩模型。此外，还将模型的似然性与收缩协方差估计器得到的似然性进行了比较。