【Python実践】データ分析の影の立役者！主成分分析（PCA）の使い方ガイド（コピペで使えるサンプルコード付き）

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt
from matplotlib import rcParams

rcParams['font.family'] = 'Meiryo'

# ------ダミーデータの生成 --------------------------------------
np.random.seed(42)  # シード設定
n_samples = 100  # サンプル数（時間軸を表現するサンプル）
time = np.linspace(0, 2 * np.pi, n_samples)  # 時間軸データ（0 ~ 2πの範囲）
temperature = 70 + 5 * np.sin(time) + np.random.normal(0, 0.5, n_samples)  # 温度は時間的なシーズン性を模倣（正弦波）
pressure = 50 + 3 * np.sin(2 * time) + np.random.normal(0, 0.5, n_samples)  # 圧力の周期性（速度や振動との連動性）
speed = 100 + 10 * np.sin(3 * time) + np.random.normal(0, 2, n_samples)  # 速度は工場稼働サイクルを模倣
vibration = 0.5 + 0.2 * np.sin(4 * time) + np.random.normal(0, 0.1, n_samples)  # 振動は稼働の影響を反映した正弦波
oil_flow = 5 + 1 * np.sin(time / 2) + np.random.normal(0, 0.3, n_samples)  # オイル流量はゆっくりとした変動を模倣

# DataFrameへまとめる
data = pd.DataFrame({
    '温度': temperature,
    '圧力': pressure,
    '速度': speed,
    '振動': vibration,
    'オイル流量': oil_flow,
})
# ---------------------------------------------

# 標準化
scaler = StandardScaler()
scaled_data = scaler.fit_transform(data)

# PCAの適用（必要な成分数を指定しない）
pca = PCA()
pca_data = pca.fit_transform(scaled_data)

# 累積寄与率を計算
cumulative_explained_variance = np.cumsum(pca.explained_variance_ratio_)

# 結果表示
print("各主成分の寄与率:")
print(pca.explained_variance_ratio_)

print("\n累積寄与率:")
print(cumulative_explained_variance)

print("\n主成分の負荷量:")
loadings = pd.DataFrame(pca.components_.T, index=data.columns, columns=[f"PC{i+1}" for i in range(pca.n_components_)])
print(loadings)

# PC1 vs PC2のプロット
plt.figure(figsize=(10, 6))
plt.scatter(pca_data[:, 0], pca_data[:, 1], label="データ")
plt.xlabel("主成分1 (PC1)")
plt.ylabel("主成分2 (PC2)")
plt.title("PCA: 主成分1 vs 主成分2")
plt.legend()
plt.grid()
plt.show()

# 累積寄与率と寄与率の棒グラフと折れ線グラフの組み合わせ
plt.figure(figsize=(10, 6))
plt.bar(range(1, len(pca.explained_variance_ratio_) + 1), pca.explained_variance_ratio_, alpha=0.5, align='center', label='個別の寄与率')
plt.plot(range(1, len(cumulative_explained_variance) + 1), cumulative_explained_variance, marker='o', label='累積寄与率', linestyle='--')
plt.xlabel("主成分数")
plt.ylabel("寄与率")
plt.title("寄与率と累積寄与率のグラフ")
plt.legend()
plt.grid()
plt.show()

# 主成分負荷量のヒートマップ
plt.figure(figsize=(10, 6))
plt.imshow(loadings, cmap='viridis', aspect='auto')
plt.colorbar(label='負荷量')
plt.xticks(ticks=np.arange(pca.n_components_), labels=[f"PC{i+1}" for i in range(pca.n_components_)])
plt.yticks(ticks=np.arange(len(data.columns)), labels=data.columns)
plt.xlabel("主成分")
plt.ylabel("変数")
plt.title("主成分負荷量のヒートマップ")
plt.grid()
plt.show()

# 共通部分: データ準備
import numpy as np
import pandas as pd

# サンプルデータの生成
np.random.seed(42)
data = pd.DataFrame({
    'センサー1': np.random.normal(0, 1, 100),
    'センサー2': np.random.normal(1, 2, 100),
    'センサー3': np.random.normal(-1, 1, 100),
    'センサー4': np.random.normal(0.5, 1.5, 100),
    'センサー5': np.random.normal(-0.5, 1, 100),
})

# 必要なライブラリ
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# データのスケーリング
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)

# PCAの適用
pca = PCA(n_components=2)
data_pca = pca.fit_transform(data_scaled)

print("元のデータの次元:", data.shape)
print("削減後のデータの次元:", data_pca.shape)
print(data_pca)

# 必要なライブラリ
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# データのスケーリング
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)

# PCAで95%の分散を保持
pca = PCA(n_components=0.86)
data_reduced = pca.fit_transform(data_scaled)

print("元のデータの次元:", data.shape)
print("圧縮後のデータの次元:", data_reduced.shape)
print(data_reduced)

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# データのスケーリング
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)

# PCAの適用（80%の分散を保持）
pca = PCA(n_components=0.7)
data_reduced = pca.fit_transform(data_scaled)

print("元のデータの次元:", data.shape)
print("ノイズの除去の次元:", data_reduced.shape)

# 必要なライブラリ
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'Meiryo'

# PCAの適用（スケーリングなし）
pca = PCA(n_components=2)
data_pca = pca.fit_transform(data)

# 可視化
plt.figure(figsize=(8, 6))
plt.scatter(data_pca[:, 0], data_pca[:, 1], c='blue', alpha=0.6)
plt.title('PCAによる2次元可視化（スケーリングなし）')
plt.xlabel('主成分1')
plt.ylabel('主成分2')
plt.grid()
plt.show()

# 必要なライブラリ
from sklearn.decomposition import PCA

# PCAの適用（すべての成分を考慮）
pca = PCA(n_components=None)
pca.fit(data)

# 特徴量ごとの寄与率
print("各特徴量の寄与率:")
print(pd.DataFrame({
    '特徴量': data.columns,
    '寄与率': pca.explained_variance_ratio_
}))

from sklearn.decomposition import PCA
import pickle
import matplotlib.pyplot as plt
from matplotlib import rcParams
rcParams['font.family'] = 'Meiryo'

# ------------------------------------------------
# 正常データを用いたモデルの作成
# ------------------------------------------------
# PCAの適用
pca = PCA(n_components=2)
pca.fit(data)

# モデルの保存
with open('pca_model.pkl', 'wb') as f:
    pickle.dump(pca, f)

# ------------------------------------------------
# 新しいデータで異常検知を実施
# ------------------------------------------------
# 新しいデータの生成
new_data = pd.DataFrame({
    'センサー1': np.random.rand(10),  # 特徴量名を元のデータと一致させる
    'センサー2': np.random.rand(10),
    'センサー3': np.random.rand(10),
    'センサー4': np.random.rand(10),
    'センサー5': np.random.rand(10)
})

# PCAモデルの読み込み
with open('pca_model.pkl', 'rb') as f:
    loaded_pca = pickle.load(f)

# 新しいデータをPCAモデルで変換
transformed_new_data = loaded_pca.transform(new_data)

# 新しいデータの再構成
reconstructed_data = loaded_pca.inverse_transform(transformed_new_data)

# 再構成誤差の計算（※１）
reconstruction_error = np.mean((new_data.values - reconstructed_data)**2, axis=1)

# 再構成誤差の閾値設定（※２）。平均値に標準偏差の3倍（3σ）を加えたもの
threshold = reconstruction_error.mean() + 3 * reconstruction_error.std()

# 誤差（※１）が閾値（※2）を超えているかで異常を判定
is_anomalous = reconstruction_error > threshold

print("新しいデータの再構成誤差:")
print(reconstruction_error)
print("\n異常検知の結果（True = 異常, False = 正常）:")
print(is_anomalous)

# 可視化
plt.figure(figsize=(10, 6))
plt.hist(reconstruction_error, bins=30, alpha=0.5, label='新しいデータの再構成誤差')
plt.axvline(threshold, color='r', linestyle='--', label='異常検知の閾値')
plt.xlabel('再構成誤差')
plt.ylabel('頻度')
plt.legend()
plt.title('新しいデータの再構成誤差の分布')
plt.grid()
plt.show()

# 必要なライブラリ
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

# ラベルを仮定
np.random.seed(42)
labels = np.random.choice([0, 1], size=100)

# データのスケーリング
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)

# PCAによる次元削減
pca = PCA(n_components=2)
data_pca = pca.fit_transform(data_scaled)

#-------------------------------------------------------------------------
# 以下は、次元圧縮した結果を使って、ランダムフォレストのモデルを作成するサンプル
#-------------------------------------------------------------------------
# 学習データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(data_pca, labels, test_size=0.2, random_state=42)

# モデルの学習
model = RandomForestClassifier(random_state=42)
model.fit(X_train, y_train)

# モデルの評価
y_pred = model.predict(X_test)
print("モデルの精度:")
print(accuracy_score(y_test, y_pred))

# ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～
#　ここに後述するグラフ描画関数とPCAUtilクラスを張り付けるか、インポートしてください。
# ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

# データの生成
np.random.seed(42)
data = pd.DataFrame({
    'feature1': np.random.rand(100),
    'feature2': np.random.rand(100),
    'feature3': np.random.rand(100),
    'feature4': np.random.rand(100),
    'feature5': np.random.rand(100)
})

# PCAモデルの作成と学習
pca_util = PCAUtil(df=data)
pca_util.fit(columns=['feature1', 'feature2', 'feature3', 'feature4', 'feature5'])

print("変換されたデータ:")
print(pca_util.transformed_data)

# 寄与率の取得
print("寄与率:")
print(pca_util.explained_variance_ratio)

print("負荷量")
print(pca_util.loadings)

# グラフの描画
plot_pca_results(pca_util)

# ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～
#　ここに後述するグラフ描画関数とPCAUtilクラスを張り付けるか、インポートしてください。
# ～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

# ステップ1: データの生成
np.random.seed(42)
data = pd.DataFrame({
    'feature1': np.random.rand(100),
    'feature2': np.random.rand(100),
    'feature3': np.random.rand(100),
    'feature4': np.random.rand(100),
    'feature5': np.random.rand(100)
})

# ステップ2: PCAモデルの作成と学習
pca_util = PCAUtil(df=data)
pca_util.fit(columns=['feature1', 'feature2', 'feature3', 'feature4', 'feature5'])

# モデルの保存
model_path = 'pca_model.pkl'
pca_util.save_model(model_path)

# ステップ3: 新しいデータの生成
new_data = pd.DataFrame({
    'feature1': np.random.rand(100),
    'feature2': np.random.rand(100),
    'feature3': np.random.rand(100),
    'feature4': np.random.rand(100),
    'feature5': np.random.rand(100)
})

# 新しいデータをPCAUtilクラスに設定
pca_util_new = PCAUtil(df=new_data)

# 保存したモデルの読み込み
pca_util_new.load_model(model_path)

# ステップ4: 異常検知
anomalies, reconstruction_error, threshold = pca_util_new.detect_anomalies(
    columns=['feature1', 'feature2', 'feature3', 'feature4', 'feature5'])

# ステップ5: 結果の表示
print("再構成誤差:")
print(reconstruction_error)
print("\n異常検知の結果（True = 異常, False = 正常）:")
print(anomalies)
print("\n異常検知の閾値:")
print(threshold)

# ステップ6: グラフの描画
plot_pca_results(pca_util_new)
plot_reconstruction_error_hist(reconstruction_error, threshold)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from sklearn.decomposition import PCA
rcParams['font.family'] = 'Meiryo'

# PCAの結果をプロットする関数
def plot_pca_results(pca):
    """
    PCAの結果をプロットする関数。
    
    Parameters:
    pca : PCA
        学習済みのPCAモデル。
    """

    fig, axes = plt.subplots(1, 3, figsize=(18, 6))

    # PC1 vs PC2の散布図
    scatter = axes[0].scatter(pca.transformed_data[:, 0], pca.transformed_data[:, 1], alpha=0.7, edgecolors='k', s=100)
    axes[0].set_xlabel('主成分1 (PC1)')
    axes[0].set_ylabel('主成分2 (PC2)')
    axes[0].set_title('主成分1と主成分2の散布図')

    # 寄与率と累積寄与率の棒グラフと折れ線グラフの組み合わせ
    cumulative_explained_variance = np.cumsum(pca.explained_variance_ratio)
    axes[1].bar(range(1, len(pca.explained_variance_ratio) + 1), pca.explained_variance_ratio, alpha=0.5, align='center', label='個別の寄与率')
    axes[1].plot(range(1, len(cumulative_explained_variance) + 1), cumulative_explained_variance, marker='o', label='累積寄与率', linestyle='--')
    axes[1].set_xlabel("主成分数")
    axes[1].set_ylabel("寄与率")
    axes[1].set_title("寄与率と累積寄与率のグラフ")
    axes[1].legend()

    # 主成分負荷量のヒートマップ
    cax = axes[2].imshow(pca.loadings, cmap='viridis', aspect='auto')
    plt.colorbar(cax, ax=axes[2], orientation='vertical', label='負荷量')
    axes[2].set_xticks(np.arange(pca.n_components))
    axes[2].set_xticklabels([f"PC{i+1}" for i in range(pca.n_components)])
    axes[2].set_yticks(np.arange(len(pca.columns)))
    axes[2].set_yticklabels([f"Feature {i+1}" for i in range(pca.df.shape[1])])
    axes[2].set_xlabel("主成分")
    axes[2].set_ylabel("特徴量")
    axes[2].set_title("主成分負荷量のヒートマップ")

    plt.tight_layout()
    plt.show()

# 再構成誤差のヒストグラムをプロットする関数
def plot_reconstruction_error_hist(reconstruction_error, threshold):
    """
    再構成誤差のヒストグラムをプロットする関数。
    
    Parameters:
    reconstruction_error : array-like
        再構成誤差の配列。
    threshold : float
        異常検知の閾値。
    """
    plt.figure(figsize=(10, 6))
    plt.hist(reconstruction_error, bins=30, alpha=0.5, label='再構成誤差')
    plt.axvline(threshold, color='r', linestyle='--', label='異常検知の閾値')
    plt.xlabel('再構成誤差')
    plt.ylabel('頻度')
    plt.legend()
    plt.title('再構成誤差の分布')
    plt.grid()
    plt.show()

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import pickle

# PCAUtilクラス
class PCAUtil:
    def __init__(self, df=None, model_path=None):
        """
        PCAUtilクラスの初期化メソッド。
        
        Parameters:
        df : pandas.DataFrame, optional
            分析に使用するデータフレーム。デフォルトはNone。
        model_path : str, optional
            保存済みPCAモデルのパス。デフォルトはNone。
        """
        self.df = None if df is None else df.copy()
        self.model = None if model_path is None else self.load_model(model_path)
        self.transformed_data = None # PCAで変換されたデータ
        self.explained_variance_ratio = None # PCAの寄与率
        self.loadings = None # 負荷量
        self.n_components = None # 主成分の数 
        self.columns = []
        self.scaler = StandardScaler()

    def fit(self, columns, n_components=2, scale=True):
        """
        PCAモデルの作成と学習を行うメソッド。
        
        Parameters:
        columns : list of str
            PCAを適用するデータフレームのカラム名のリスト。
        n_components : int, optional
            PCAで抽出する主成分の数。デフォルトは2。
        scale : bool, optional
            データをスケーリングするかどうか。デフォルトはTrue。
        
        Returns:
        PCA
            学習済みのPCAモデル。
        """
        self.columns = columns
        X = self.df[columns].values
        if scale:
            X = self.scaler.fit_transform(X)
        self.model = PCA(n_components=n_components).fit(X)
        self.transform(columns, scale=scale)
        return self.model

    def transform(self, columns=None, scale=True):
        """
        学習済みモデルを使用してデータを変換するメソッド。
        
        Parameters:
        columns : list of str
            変換するデータフレームのカラム名のリスト。
            省略するとモデル作成時に指定したカラム名のリストが使用される。
        scale : bool, optional
            データをスケーリングするかどうか。デフォルトはTrue。
        
        Returns:
        numpy.ndarray
            変換されたデータ。
        """
        if columns is None:
            columns = self.columns
        else:
            self.columns = columns

        X = self.df[columns].values
        if scale:
            X = self.scaler.fit_transform(X)
        self.transformed_data = self.model.transform(X)
        self.explained_variance_ratio = self.model.explained_variance_ratio_
        self.n_components = self.model.n_components_
        self.loadings = pd.DataFrame(self.model.components_.T, 
            index=self.columns, columns=[f"PC{i+1}" for i in range(self.n_components)])
        return self.transformed_data
     
    def read_csv(self, file_name, encoding="shift-jis"):
        """
        CSVファイルを読み込むメソッド。
        
        Parameters:
        file_name : str
            読み込むCSVファイルのパス。
        encoding : str, optional
            ファイルのエンコーディング。デフォルトは"shift-jis"。
        """
        self.df = pd.read_csv(file_name, encoding=encoding)

    def save_model(self, model_path):
        """
        PCAモデルをファイルに保存するメソッド。
        
        Parameters:
        model_path : str
            保存するファイルのパス。
        """
        self.model_path = model_path
        with open(model_path, 'wb') as f:
            pickle.dump(self.model, f)

    def load_model(self, model_path=None):
        """
        保存済みPCAモデルを読み込むメソッド。
        
        Parameters:
        model_path : str, optional
            読み込むPCAモデルのファイルのパス。デフォルトはNone。
        
        Returns:
        PCA
            読み込まれたPCAモデル。
        """
        model_path = model_path if model_path else self.model_path
        with open(model_path, 'rb') as f:
            self.model = pickle.load(f)
        return self.model

    def detect_anomalies(self, columns=None, threshold=None, scale=True):
        """
        データの異常検知を行うメソッド。
        
        Parameters:
        columns : list of str
            異常検知を行うデータフレームのカラム名のリスト。
            省略すると、fit 又は transform で指定したカラム名のリストが使用される。
        threshold : float, optional
            異常とみなす再構成誤差の閾値。デフォルトはNoneで、
            再構成誤差の平均 + 3 * 標準偏差が設定されます。
        scale : bool, optional
            データをスケーリングするかどうか。デフォルトはTrue。
        
        Returns:
        tuple
            異常フラグの配列、再構成誤差の配列、閾値。
        """
        # データを変換
        columns = self.columns if columns is None else columns
        self.transform(columns, scale=scale)
        
        # 再構成
        reconstructed_data = self.model.inverse_transform(self.transformed_data)
        if scale:
            reconstructed_data = self.scaler.inverse_transform(reconstructed_data)
        
        # 再構成誤差を計算
        reconstruction_error = np.mean((self.df[columns].values - reconstructed_data)**2, axis=1)
        
        # 閾値を設定（デフォルトは再構成誤差の平均 + 3 * 標準偏差）
        if threshold is None:
            threshold = reconstruction_error.mean() + 3 * reconstruction_error.std()
        
        # 異常検知
        anomalies = reconstruction_error > threshold
        
        return anomalies, reconstruction_error, threshold