# 必要なライブラリのインポート import numpy as np import pandas as pd import matplotlib.pyplot as plt import seaborn as sns from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA from sklearn.manifold import TSNE from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score from sklearn.ensemble import IsolationForest # 日本語フォント設定（Google Colab用） plt.rcParams[‘font.family’] = ‘DejaVu Sans’ # 乱数シード固定 np.random.seed(42)

# 顧客データの生成（実務ではCSVから読み込む） n_customers = 500 # RFM特徴量 + その他の特徴量 data = { ‘customer_id’: range(1, n_customers + 1), ‘recency’: np.random.exponential(30, n_customers), # 最終購入からの日数 ‘frequency’: np.random.poisson(5, n_customers) + 1, # 購入回数 ‘monetary’: np.random.exponential(10000, n_customers), # 合計購入金額 ‘avg_order_value’: np.random.exponential(2000, n_customers), # 平均注文金額 ‘days_since_first’: np.random.uniform(30, 365, n_customers), # 初回購入からの日数 ‘product_categories’: np.random.poisson(3, n_customers) + 1, # 購入カテゴリ数 ‘returns’: np.random.poisson(0.5, n_customers), # 返品回数 ‘customer_service_calls’: np.random.poisson(1, n_customers), # 問い合わせ回数 } df = pd.DataFrame(data) # 基本統計量 print(“=”*60) print(“顧客データの基本情報”) print(“=”*60) print(f”\nデータサイズ: {df.shape[0]}顧客 × {df.shape[1]-1}特徴量”) print(“\n基本統計量:”) print(df.describe().round(2))

# 特徴量の分布を可視化 features = [‘recency’, ‘frequency’, ‘monetary’, ‘avg_order_value’] fig, axes = plt.subplots(2, 2, figsize=(12, 10)) axes = axes.flatten() for i, col in enumerate(features): axes[i].hist(df[col], bins=30, edgecolor=’black’, alpha=0.7) axes[i].set_title(f'{col} の分布’) axes[i].set_xlabel(col) axes[i].set_ylabel(‘頻度’) plt.tight_layout() plt.show()

# 相関行列 feature_cols = [‘recency’, ‘frequency’, ‘monetary’, ‘avg_order_value’, ‘days_since_first’, ‘product_categories’, ‘returns’, ‘customer_service_calls’] plt.figure(figsize=(10, 8)) correlation = df[feature_cols].corr() sns.heatmap(correlation, annot=True, cmap=’coolwarm’, center=0, fmt=’.2f’) plt.title(‘特徴量の相関行列’) plt.tight_layout() plt.show()

# 分析用の特徴量を選択 X = df[feature_cols].values print(f”特徴量行列: {X.shape}”) # 標準化（必須！） scaler = StandardScaler() X_scaled = scaler.fit_transform(X) print(f”\n標準化後の平均: {X_scaled.mean(axis=0).round(2)}”) print(f”標準化後の標準偏差: {X_scaled.std(axis=0).round(2)}”)

# PCAで2次元に削減 pca = PCA(n_components=2) X_pca = pca.fit_transform(X_scaled) print(“PCA結果:”) print(f” 第1主成分の寄与率: {pca.explained_variance_ratio_[0]:.2%}”) print(f” 第2主成分の寄与率: {pca.explained_variance_ratio_[1]:.2%}”) print(f” 累積寄与率: {sum(pca.explained_variance_ratio_):.2%}”) # 全主成分の寄与率を確認 pca_full = PCA() pca_full.fit(X_scaled) cumsum = np.cumsum(pca_full.explained_variance_ratio_) print(f”\n累積寄与率が90%になる次元数: {np.argmax(cumsum >= 0.9) + 1}”)

# t-SNEで2次元に削減 tsne = TSNE(n_components=2, random_state=42, perplexity=30) X_tsne = tsne.fit_transform(X_scaled) # PCAとt-SNEの比較 fig, axes = plt.subplots(1, 2, figsize=(14, 6)) axes[0].scatter(X_pca[:, 0], X_pca[:, 1], s=30, alpha=0.6) axes[0].set_title(‘PCA（2次元）’) axes[0].set_xlabel(‘PC1’) axes[0].set_ylabel(‘PC2’) axes[1].scatter(X_tsne[:, 0], X_tsne[:, 1], s=30, alpha=0.6) axes[1].set_title(‘t-SNE（2次元）’) axes[1].set_xlabel(‘t-SNE 1’) axes[1].set_ylabel(‘t-SNE 2’) plt.tight_layout() plt.show()

# 最適なクラスタ数を決定（エルボー法 + シルエット分析） inertias = [] silhouettes = [] K_range = range(2, 11) for k in K_range: kmeans = KMeans(n_clusters=k, random_state=42, n_init=10) labels = kmeans.fit_predict(X_scaled) inertias.append(kmeans.inertia_) silhouettes.append(silhouette_score(X_scaled, labels)) # プロット fig, axes = plt.subplots(1, 2, figsize=(12, 4)) axes[0].plot(K_range, inertias, ‘bo-‘, linewidth=2) axes[0].set_xlabel(‘クラスタ数 (K)’) axes[0].set_ylabel(‘イナーシャ’) axes[0].set_title(‘エルボー法’) axes[0].set_xticks(K_range) axes[1].plot(K_range, silhouettes, ‘go-‘, linewidth=2) axes[1].set_xlabel(‘クラスタ数 (K)’) axes[1].set_ylabel(‘シルエットスコア’) axes[1].set_title(‘シルエット分析’) axes[1].set_xticks(K_range) plt.tight_layout() plt.show() # 結果表示 best_k = K_range[np.argmax(silhouettes)] print(f”シルエットスコアが最大のK: {best_k}”)

# K=4でクラスタリング実行 n_clusters = 4 kmeans = KMeans(n_clusters=n_clusters, random_state=42, n_init=10) df[‘cluster’] = kmeans.fit_predict(X_scaled) print(“クラスタごとの顧客数:”) print(df[‘cluster’].value_counts().sort_index()) # クラスタリング結果の可視化（t-SNE上） plt.figure(figsize=(10, 8)) scatter = plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=df[‘cluster’], cmap=’viridis’, s=50, alpha=0.7) plt.colorbar(scatter, label=’Cluster’) plt.title(‘顧客セグメント（t-SNE + K-means）’) plt.xlabel(‘t-SNE 1’) plt.ylabel(‘t-SNE 2’) plt.show()

# クラスタごとの特徴量の平均値 cluster_summary = df.groupby(‘cluster’)[feature_cols].mean().round(2) print(“=”*70) print(“クラスタごとの特徴量平均”) print(“=”*70) print(cluster_summary)

# 各クラスタの特徴を可視化（レーダーチャート風） # 標準化した値でクラスタ比較 cluster_means_scaled = pd.DataFrame( scaler.transform(cluster_summary), columns=feature_cols, index=cluster_summary.index ) # ヒートマップで可視化 plt.figure(figsize=(12, 6)) sns.heatmap(cluster_means_scaled.T, annot=True, cmap=’RdYlGn’, center=0, fmt=’.2f’) plt.title(‘クラスタ特性の比較（標準化後、0=平均）’) plt.xlabel(‘クラスタ’) plt.ylabel(‘特徴量’) plt.tight_layout() plt.show()

# クラスタの命名と解釈 cluster_names = { 0: ‘優良顧客（高頻度・高額）’, 1: ‘離脱リスク顧客（低頻度・低額・休眠）’, 2: ‘ロイヤル顧客（長期・高額）’, 3: ‘新規顧客（短期・低額）’ } print(“=”*60) print(“クラスタの解釈”) print(“=”*60) for cluster_id in range(n_clusters): cluster_data = df[df[‘cluster’] == cluster_id] print(f”\n【クラスタ {cluster_id}: {cluster_names[cluster_id]}】”) print(f” 顧客数: {len(cluster_data)}人 ({len(cluster_data)/len(df)*100:.1f}%)”) print(f” 平均購入頻度: {cluster_data[‘frequency’].mean():.1f}回”) print(f” 平均購入金額: ¥{cluster_data[‘monetary’].mean():,.0f}”) print(f” 最終購入からの日数: {cluster_data[‘recency’].mean():.1f}日”)

# Isolation Forestで異常検知 iso_forest = IsolationForest(contamination=0.05, random_state=42) df[‘anomaly’] = iso_forest.fit_predict(X_scaled) # 異常スコア df[‘anomaly_score’] = iso_forest.decision_function(X_scaled) # 結果 n_anomaly = sum(df[‘anomaly’] == -1) print(f”異常と判定された顧客: {n_anomaly}人 ({n_anomaly/len(df)*100:.1f}%)”) # 異常顧客の特徴 print(“\n異常顧客の特徴量平均:”) anomaly_summary = df[df[‘anomaly’] == -1][feature_cols].mean() normal_summary = df[df[‘anomaly’] == 1][feature_cols].mean() comparison = pd.DataFrame({ ‘異常顧客’: anomaly_summary, ‘正常顧客’: normal_summary, ‘差（異常-正常）’: anomaly_summary – normal_summary }).round(2) print(comparison)

# 異常顧客の可視化 plt.figure(figsize=(10, 8)) colors = [‘red’ if a == -1 else ‘blue’ for a in df[‘anomaly’]] plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=colors, s=50, alpha=0.6) plt.title(‘異常顧客の検出（赤=異常、青=正常）’) plt.xlabel(‘t-SNE 1’) plt.ylabel(‘t-SNE 2’) plt.show() # 異常顧客のリスト print(“\n異常顧客TOP5（異常スコアが低い順）:”) anomaly_customers = df[df[‘anomaly’] == -1].nsmallest(5, ‘anomaly_score’) print(anomaly_customers[[‘customer_id’, ‘monetary’, ‘frequency’, ‘returns’, ‘anomaly_score’]])

# 最終レポートの作成 print(“=”*70) print(“顧客セグメンテーション分析レポート”) print(“=”*70) print(“\n【1. データ概要】”) print(f” 分析対象: {len(df)}顧客”) print(f” 特徴量数: {len(feature_cols)}項目”) print(“\n【2. クラスタリング結果】”) print(f” クラスタ数: {n_clusters}”) print(f” シルエットスコア: {silhouette_score(X_scaled, df[‘cluster’]):.3f}”) for c in range(n_clusters): count = sum(df[‘cluster’] == c) print(f” クラスタ{c} ({cluster_names[c]}): {count}人 ({count/len(df)*100:.1f}%)”) print(“\n【3. 異常検知結果】”) print(f” 異常顧客数: {n_anomaly}人 ({n_anomaly/len(df)*100:.1f}%)”) print(“\n【4. 推奨アクション】”) print(” クラスタ0: ロイヤリティプログラムの強化”) print(” クラスタ1: リテンションキャンペーンの実施”) print(” クラスタ2: クロスセル・アップセルの提案”) print(” クラスタ3: オンボーディングの改善”) print(” 異常顧客: 個別確認と対応”)