PCA と FA のモデル選択をプロットする

はじめに

この実験では、確率的主成分分析（Probabilistic PCA）と因子分析（Factor Analysis）の 2 つの確率モデルを検討し、モデル選択と共分散推定におけるそれらの有効性を比較します。同分散または異分散ノイズが含まれる低ランクデータに対して交差検証を行います。また、シュリンク共分散推定器から得られる尤度と比較して、モデルの尤度を求めます。

VM のヒント

VM の起動が完了したら、左上隅をクリックしてノートブックタブに切り替え、Jupyter Notebook を使って練習しましょう。

Jupyter Notebook の読み込みには数秒かかる場合があります。Jupyter Notebook の制限により、操作の検証は自動化できません。

学習中に問題がある場合は、Labby にお問い合わせください。セッション後にフィードバックを提供してください。すぐに問題を解決いたします。

Skills Graph

データの作成

500 個のサンプル、25 個の特徴量、ランク 5 の疑似データセットを作成します。また、同分散および異分散ノイズをデータセットに追加します。

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

モデルの適合

確率的主成分分析（Probabilistic PCA）と因子分析（Factor Analysis）のモデルをデータセットに適合させ、交差検証を使ってそれらの性能を評価します。また、シュリンク共分散推定器のスコアを計算し、結果を比較します。

import matplotlib.pyplot as plt
from sklearn.decomposition import PCA, FactorAnalysis
from sklearn.covariance import ShrunkCovariance, LedoitWolf
from sklearn.model_selection import cross_val_score, GridSearchCV

n_components = np.arange(0, n_features, 5)  ## n_components のオプション

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("PCA CV による最適な n_components = %d" % n_components_pca)
    print("FactorAnalysis CV による最適な n_components = %d" % n_components_fa)
    print("PCA MLE による最適な n_components = %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="FactorAnalysis 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" % n_components_pca_mle,
        linestyle="-.",
    )

    plt.xlabel("成分数")
    plt.ylabel("CV スコア")
    plt.legend(loc="lower right")
    plt.title(title)

plt.show()

まとめ

この実験では、確率的主成分分析（Probabilistic PCA）と因子分析（Factor Analysis）のモデルがモデル選択と共分散推定における有効性を検討しました。同分散および異分散ノイズを含む疑似データセットを作成し、交差検証を使ってモデルの性能を比較しました。また、シュリンク共分散推定器から得られる尤度と比較して、モデルの尤度を求めました。結果は、同分散ノイズの存在下では PCA と FA の両方が低ランク部分空間のサイズを回復するのに有効であることを示しています。ただし、異分散ノイズが存在する場合、PCA は失敗し、ランクを過大評価しました。適切な状況下では、低ランクモデルのための予備データは、シュリンクモデルよりも尤度が高くなります。