【Python実践】データに潜む共通点を見つける！因子分析（FA)の使い方ガイド（コピペで使えるサンプルコード付き）

import numpy as np
import pandas as pd
from factor_analyzer import FactorAnalyzer
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)

# 因子分析の適用（2因子を抽出）
fa = FactorAnalyzer(n_factors=2, rotation='varimax')
fa.fit(scaled_data)

# 因子負荷量を取得
loadings = fa.loadings_
print("\n因子負荷量:")
loadings_df = pd.DataFrame(loadings, index=data.columns, columns=[f"Factor{i+1}" for i in range(loadings.shape[1])])
print(loadings_df)

# 共通性（Communalities）の計算
communalities = fa.get_communalities()
print("\n共通性:")
communalities_df = pd.DataFrame(communalities, index=data.columns, columns=['共通性'])
print(communalities_df)

# 固有値（Eigenvalues）と寄与率の計算
eigenvalues, variance_explained, cumulative_explained_variance = fa.get_factor_variance()
print("\n各因子の固有値:")
print(eigenvalues)

print("\n各因子の寄与率:")
print(variance_explained)

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

# 因子得点の計算
factor_scores = fa.transform(scaled_data)
factor_scores_df = pd.DataFrame(factor_scores, columns=[f"Factor{i+1}" for i in range(fa.n_factors)])
print("\n因子得点（先頭の5点）:")
print(factor_scores_df.head())

# Factor1 vs Factor2のプロット
plt.figure(figsize=(10, 6))
plt.scatter(factor_scores[:, 0], factor_scores[:, 1], label="データ")
plt.xlabel("因子1 (Factor1)")
plt.ylabel("因子2 (Factor2)")
plt.title("因子分析: 因子1 vs 因子2")
plt.legend()
plt.grid()
plt.show()

# 累積寄与率と寄与率の棒グラフと折れ線グラフの組み合わせ
plt.figure(figsize=(10, 6))
plt.bar(range(1, len(variance_explained) + 1), variance_explained, 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(loadings.shape[1]), labels=[f"Factor{i+1}" for i in range(loadings.shape[1])])
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 factor_analyzer import FactorAnalyzer
from sklearn.preprocessing import StandardScaler

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

# 因子分析の適用（2因子を抽出）
fa = FactorAnalyzer(n_factors=2, rotation='varimax')
fa.fit(data_scaled)

# 因子得点の計算
factor_scores = fa.transform(data_scaled)
print("元のデータの次元:", data.shape)
print("削減後のデータの次元:", factor_scores.shape)
print(factor_scores)

# 必要なライブラリ
from factor_analyzer import FactorAnalyzer
from sklearn.preprocessing import StandardScaler

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

# 因子分析の適用（寄与率に基づく因子の数の選択）
# ここでは仮に、累積寄与率が0.86に達する因子数を選びます
fa = FactorAnalyzer(rotation='varimax')
fa.fit(data_scaled)
eigenvalues, variance_explained, cumulative_explained_variance = fa.get_factor_variance()

# 累積寄与率が0.86以上になる因子の数を確認
n_factors = np.argmax(cumulative_explained_variance >= 0.86) + 1

# 選ばれた因子数で再度因子分析を実行
fa = FactorAnalyzer(n_factors=n_factors, rotation='varimax')
fa.fit(data_scaled)

# 因子得点の計算
factor_scores = fa.transform(data_scaled)
print("元のデータの次元:", data.shape)
print("圧縮後のデータの次元:", factor_scores.shape)
print(factor_scores)

# 必要なライブラリ
from factor_analyzer import FactorAnalyzer
from sklearn.preprocessing import StandardScaler

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

# 因子分析の適用（寄与率に基づく因子の数の選択）
# ここでは仮に、累積寄与率が0.7に達する因子数を選びます
fa = FactorAnalyzer(rotation='varimax')
fa.fit(data_scaled)
eigenvalues, variance_explained, cumulative_explained_variance = fa.get_factor_variance()

# 累積寄与率が0.7以上になる因子の数を確認
n_factors = np.argmax(cumulative_explained_variance >= 0.7) + 1

# 選ばれた因子数で再度因子分析を実行
fa = FactorAnalyzer(n_factors=n_factors, rotation='varimax')
fa.fit(data_scaled)

# 因子得点の計算
factor_scores = fa.transform(data_scaled)
print("元のデータの次元:", data.shape)
print("圧縮後のデータの次元:", factor_scores.shape)
print(factor_scores)

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

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

# 因子分析の適用（2因子を抽出）
fa = FactorAnalyzer(n_factors=2, rotation='varimax')
fa.fit(data_scaled)

# 因子得点の計算
factor_scores = fa.transform(data_scaled)

# 必要なライブラリ
from factor_analyzer import FactorAnalyzer
import pandas as pd

# 因子分析の適用（すべての因子を考慮）
fa = FactorAnalyzer(rotation='varimax')
fa.fit(data)

# 特徴量ごとの寄与率（因子負荷量の平均値として計算）
loadings = fa.loadings_
average_loadings = loadings.mean(axis=1)

print("各特徴量の寄与率（因子負荷量の平均）:")
print(pd.DataFrame({
    '特徴量': data.columns,
    '寄与率': average_loadings
}))

# 必要なライブラリ
from factor_analyzer import FactorAnalyzer
import pandas as pd
import numpy as np
import pickle
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from matplotlib import rcParams
rcParams['font.family'] = 'Meiryo'

# ------------------------------------------------
# 正常データを用いたモデルの作成
# ------------------------------------------------
# データのスケーリング
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)

# 因子分析の適用（2因子を抽出）
fa = FactorAnalyzer(n_factors=2, rotation='varimax')
fa.fit(data_scaled)

# 因子負荷量の取得
loadings = fa.loadings_

# モデルの保存
with open('fa_model.pkl', 'wb') as f:
    pickle.dump({'model': fa, 'scaler': scaler, 'loadings': loadings}, 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)
})

# データのスケーリング
new_data_scaled = scaler.transform(new_data)

# 因子分析モデルの読み込み
with open('fa_model.pkl', 'rb') as f:
    saved = pickle.load(f)
    loaded_fa = saved['model']
    loaded_scaler = saved['scaler']
    loadings = saved['loadings']

# 新しいデータを因子分析モデルで変換
transformed_new_data = loaded_fa.transform(new_data_scaled)

# 再構成誤差の代わりに因子負荷量を用いた再構成
reconstructed_data = np.dot(transformed_new_data, loadings.T)

# スケーリングを元に戻す
reconstructed_data_original_scale = loaded_scaler.inverse_transform(reconstructed_data)

# 再構成誤差の計算（※１）
reconstruction_error = np.mean((new_data.values - reconstructed_data_original_scale)**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 factor_analyzer import FactorAnalyzer
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import numpy as np
import pandas as pd

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

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

# 因子分析による次元削減（2因子を抽出）
fa = FactorAnalyzer(n_factors=2, rotation='varimax')
fa.fit(data_scaled)
data_fa = fa.transform(data_scaled)

#-------------------------------------------------------------------------
# 以下は、次元圧縮した結果を使って、ランダムフォレストのモデルを作成するサンプル
#-------------------------------------------------------------------------
# 学習データとテストデータに分割
X_train, X_test, y_train, y_test = train_test_split(data_fa, 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))

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

# ステップ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: 因子分析モデルの作成と学習
fa_util = FAUtil(df=data)
fa_util.fit(columns=['feature1', 'feature2', 'feature3', 'feature4', 'feature5'])

# 共通性の取得
communalities = fa_util.model.get_communalities()
print("共通性:")
print(pd.DataFrame(communalities, index=fa_util.columns, columns=['共通性']))

# 固有値と寄与率の取得
eigenvalues, variance_explained, cumulative_explained_variance = fa_util.model.get_factor_variance()
print("固有値:")
print(eigenvalues)
print("寄与率:")
print(variance_explained)

# 因子負荷量の取得
loadings = fa_util.loadings
print("因子負荷量:")
print(loadings)

# グラフの描画
plot_fa_results(fa_util)

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

# ステップ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: 因子分析モデルの作成と学習
fa_util = FAUtil(df=data)
fa_util.fit(columns=['feature1', 'feature2', 'feature3', 'feature4', 'feature5'])

# モデルの保存
model_path = 'fa_model.pkl'
fa_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)
})

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

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

# ステップ4: 異常検知
anomalies, reconstruction_error, threshold = fa_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_reconstruction_error_hist(reconstruction_error, threshold)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
from factor_analyzer import FactorAnalyzer
rcParams['font.family'] = 'Meiryo'

# 因子分析の結果をプロットする関数
def plot_fa_results(fa_util):
    """
    因子分析の結果をプロットする関数。
    
    Parameters:
    fa_util : FAUtil
        学習済みのFAUtilモデル。
    """

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

    # 因子1 vs 因子2の散布図
    scatter = axes[0].scatter(fa_util.transformed_data[:, 0], fa_util.transformed_data[:, 1], alpha=0.7, edgecolors='k', s=100)
    axes[0].set_xlabel('因子1')
    axes[0].set_ylabel('因子2')
    axes[0].set_title('因子1と因子2の散布図')

    # 寄与率と累積寄与率の棒グラフと折れ線グラフの組み合わせ
    eigenvalues, variance_explained, cumulative_explained_variance = fa_util.model.get_factor_variance()
    axes[1].bar(range(1, len(variance_explained) + 1), variance_explained, 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(fa_util.loadings, cmap='viridis', aspect='auto')
    plt.colorbar(cax, ax=axes[2], orientation='vertical', label='負荷量')
    axes[2].set_xticks(np.arange(fa_util.n_factors))
    axes[2].set_xticklabels([f"因子{i+1}" for i in range(fa_util.n_factors)])
    axes[2].set_yticks(np.arange(len(fa_util.columns)))
    axes[2].set_yticklabels([f"Feature {i+1}" for i in range(fa_util.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 factor_analyzer import FactorAnalyzer
from sklearn.preprocessing import StandardScaler
import pickle

# FAUtilクラス
class FAUtil:
    def __init__(self, df=None, model_path=None):
        """
        FAUtilクラスの初期化メソッド。
        
        Parameters:
        df : pandas.DataFrame, optional
            分析に使用するデータフレーム。デフォルトはNone。
        model_path : str, optional
            保存済み因子分析モデルのパス。デフォルトは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 # 因子分析で変換されたデータ
        self.explained_variance_ratio = None # 寄与率
        self.loadings = None # 因子負荷量
        self.n_factors = None # 因子の数 
        self.columns = []
        self.scaler = StandardScaler()

    def fit(self, columns, n_factors=2, scale=True):
        """
        因子分析モデルの作成と学習を行うメソッド。
        
        Parameters:
        columns : list of str
            因子分析を適用するデータフレームのカラム名のリスト。
        n_factors : int, optional
            因子分析で抽出する因子の数。デフォルトは2。
        scale : bool, optional
            データをスケーリングするかどうか。デフォルトはTrue。
        
        Returns:
        FactorAnalyzer
            学習済みの因子分析モデル。
        """
        self.columns = columns
        X = self.df[columns].values
        if scale:
            X = self.scaler.fit_transform(X)
        self.model = FactorAnalyzer(n_factors=n_factors, rotation='varimax').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)
        eigenvalues, variance_explained, cumulative_explained_variance = self.model.get_factor_variance()
        self.explained_variance_ratio = variance_explained
        self.n_factors = self.model.n_factors
        self.loadings = pd.DataFrame(self.model.loadings_, 
            index=self.columns, columns=[f"因子{i+1}" for i in range(self.n_factors)])
        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):
        """
        因子分析モデルをファイルに保存するメソッド。
        
        Parameters:
        model_path : str
            保存するファイルのパス。
        """
        self.model_path = model_path
        with open(model_path, 'wb') as f:
            pickle.dump({
                'model': self.model,
                'scaler': self.scaler,
                'loadings': self.loadings
            }, f)

    def load_model(self, model_path=None):
        """
        保存済み因子分析モデルを読み込むメソッド。
        
        Parameters:
        model_path : str, optional
            読み込む因子分析モデルのファイルのパス。デフォルトはNone。
        
        Returns:
        FactorAnalyzer
            読み込まれた因子分析モデル。
        """
        model_path = model_path if model_path else self.model_path
        with open(model_path, 'rb') as f:
            saved = pickle.load(f)
            self.model = saved['model']
            self.scaler = saved['scaler']
            self.loadings = saved['loadings']
        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 = np.dot(self.transformed_data, self.loadings.T)
        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