import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_digits
 
# 手書き数字データ（64次元）
digits = load_digits()
X, y = digits.data, digits.target
 
# スケーリング
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
 
# PCA: 全成分で固有値分析
pca_full = PCA()
pca_full.fit(X_scaled)
 
# 累積寄与率プロット
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
 
ax1.bar(range(1, len(pca_full.explained_variance_ratio_) + 1),
        pca_full.explained_variance_ratio_, alpha=0.6, label="個別寄与率")
ax1.plot(range(1, len(pca_full.explained_variance_ratio_) + 1),
         np.cumsum(pca_full.explained_variance_ratio_), "ro-", label="累積寄与率")
ax1.axhline(y=0.95, color="g", linestyle="--", label="95%ライン")
ax1.set_xlabel("主成分番号")
ax1.set_ylabel("寄与率")
ax1.set_title("PCA 寄与率分析")
ax1.legend()
ax1.grid(True, alpha=0.3)
 
# 95%を満たす次元数
n_components_95 = np.argmax(np.cumsum(pca_full.explained_variance_ratio_) >= 0.95) + 1
print(f"95%の分散を保持する次元数: {n_components_95} / {X.shape[1]}")
 
# 2Dに射影して可視化
pca_2d = PCA(n_components=2)
X_2d = pca_2d.fit_transform(X_scaled)
 
scatter = ax2.scatter(X_2d[:, 0], X_2d[:, 1], c=y, cmap="tab10", s=5, alpha=0.6)
ax2.set_xlabel(f"PC1 ({pca_2d.explained_variance_ratio_[0]:.1%})")
ax2.set_ylabel(f"PC2 ({pca_2d.explained_variance_ratio_[1]:.1%})")
ax2.set_title("PCA 2D射影（手書き数字）")
plt.colorbar(scatter, ax=ax2)
 
plt.tight_layout()
plt.savefig("reports/pca_analysis.png", dpi=150)
plt.close()

import numpy as np
 
class PCAFromScratch:
    """
    PCA をゼロから実装する。
 
    なぜ自作するのか:
    - 共分散行列 → 固有値分解の流れを理解する
    - scikit-learn の PCA が内部で何をしているかを把握する
    - SVD ベースの実装との違いを体感する
    """
 
    def __init__(self, n_components: int = 2, method: str = "eigen"):
        """
        Parameters:
            n_components: 保持する主成分数
            method: "eigen"（固有値分解）or "svd"（特異値分解）
        """
        self.n_components = n_components
        self.method = method
        self.components_ = None       # 主成分方向 (k x d)
        self.explained_variance_ = None
        self.explained_variance_ratio_ = None
        self.mean_ = None
 
    def fit(self, X: np.ndarray) -> "PCAFromScratch":
        n_samples, n_features = X.shape
 
        # Step 1: 中心化
        self.mean_ = np.mean(X, axis=0)
        X_centered = X - self.mean_
 
        if self.method == "eigen":
            # Step 2a: 共分散行列の計算
            cov_matrix = (1 / (n_samples - 1)) * X_centered.T @ X_centered
 
            # Step 3a: 固有値分解
            eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
 
            # Step 4a: 固有値の降順でソート
            sorted_idx = np.argsort(eigenvalues)[::-1]
            eigenvalues = eigenvalues[sorted_idx]
            eigenvectors = eigenvectors[:, sorted_idx]
 
            # Step 5a: 上位 k 個を選択
            self.components_ = eigenvectors[:, :self.n_components].T
            self.explained_variance_ = eigenvalues[:self.n_components]
 
        elif self.method == "svd":
            # Step 2b: SVD（共分散行列を作らない）
            U, S, Vt = np.linalg.svd(X_centered, full_matrices=False)
 
            # Step 3b: 特異値から固有値を計算
            eigenvalues = (S ** 2) / (n_samples - 1)
 
            # Step 4b: 上位 k 個を選択
            self.components_ = Vt[:self.n_components]
            self.explained_variance_ = eigenvalues[:self.n_components]
 
        # 寄与率の計算
        total_variance = np.sum(self.explained_variance_)
        # 全分散は全固有値の合計
        if self.method == "eigen":
            total_var_all = np.sum(eigenvalues)
        else:
            total_var_all = np.sum((S ** 2) / (n_samples - 1))
 
        self.explained_variance_ratio_ = self.explained_variance_ / total_var_all
 
        return self
 
    def transform(self, X: np.ndarray) -> np.ndarray:
        """データを主成分空間に射影"""
        X_centered = X - self.mean_
        return X_centered @ self.components_.T
 
    def fit_transform(self, X: np.ndarray) -> np.ndarray:
        self.fit(X)
        return self.transform(X)
 
    def inverse_transform(self, Z: np.ndarray) -> np.ndarray:
        """主成分空間から元の空間に逆射影（近似復元）"""
        return Z @ self.components_ + self.mean_
 
    def reconstruction_error(self, X: np.ndarray) -> float:
        """再構成誤差（情報損失の指標）"""
        Z = self.transform(X)
        X_reconstructed = self.inverse_transform(Z)
        return np.mean((X - X_reconstructed) ** 2)
 
# 使用例: scikit-learn との比較
from sklearn.decomposition import PCA as SklearnPCA
from sklearn.datasets import load_digits
 
digits = load_digits()
X = digits.data
 
# ゼロから実装（SVD版）
pca_scratch = PCAFromScratch(n_components=10, method="svd")
Z_scratch = pca_scratch.fit_transform(X)
 
# scikit-learn
pca_sklearn = SklearnPCA(n_components=10)
Z_sklearn = pca_sklearn.fit_transform(X)
 
# 寄与率の比較（符号の違いを除けば一致するはず）
print("=== 寄与率の比較 ===")
for i in range(5):
    print(f"  PC{i+1}: 自作={pca_scratch.explained_variance_ratio_[i]:.6f}, "
          f"sklearn={pca_sklearn.explained_variance_ratio_[i]:.6f}")
 
# 再構成誤差
error = pca_scratch.reconstruction_error(X)
print(f"\n再構成誤差 (10成分): {error:.4f}")
print(f"累積寄与率: {sum(pca_scratch.explained_variance_ratio_):.4f}")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.datasets import load_digits
 
digits = load_digits()
X = digits.data
 
# 異なる成分数での再構成を可視化
n_components_list = [2, 5, 10, 20, 40, 64]
sample_idx = 42  # 表示するサンプル
 
fig, axes = plt.subplots(2, len(n_components_list) + 1, figsize=(20, 6))
 
# 元画像
for ax in axes[:, 0]:
    ax.imshow(X[sample_idx].reshape(8, 8), cmap="gray")
    ax.set_title("Original")
    ax.axis("off")
 
for col, n_comp in enumerate(n_components_list, 1):
    pca = PCA(n_components=n_comp)
    X_pca = pca.fit_transform(X)
    X_reconstructed = pca.inverse_transform(X_pca)
 
    # 再構成画像
    axes[0][col].imshow(X_reconstructed[sample_idx].reshape(8, 8), cmap="gray")
    axes[0][col].set_title(f"k={n_comp}")
    axes[0][col].axis("off")
 
    # 差分画像
    diff = np.abs(X[sample_idx] - X_reconstructed[sample_idx])
    axes[1][col].imshow(diff.reshape(8, 8), cmap="hot")
    axes[1][col].set_title(f"Error: {np.mean(diff**2):.2f}")
    axes[1][col].axis("off")
 
axes[1][0].set_visible(False)
 
plt.suptitle("PCA 再構成: 成分数と画質の関係", fontsize=14)
plt.tight_layout()
plt.savefig("reports/pca_reconstruction.png", dpi=150)
plt.close()
 
# 再構成誤差の推移をプロット
errors = []
cumvar = []
for k in range(1, 65):
    pca = PCA(n_components=k)
    X_pca = pca.fit_transform(X)
    X_rec = pca.inverse_transform(X_pca)
    errors.append(np.mean((X - X_rec) ** 2))
    cumvar.append(sum(pca.explained_variance_ratio_))
 
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))
ax1.plot(range(1, 65), errors, "b-")
ax1.set_xlabel("主成分数")
ax1.set_ylabel("平均二乗再構成誤差")
ax1.set_title("再構成誤差 vs 主成分数")
ax1.grid(True, alpha=0.3)
 
ax2.plot(range(1, 65), cumvar, "r-")
ax2.axhline(y=0.95, color="g", linestyle="--", label="95%")
ax2.set_xlabel("主成分数")
ax2.set_ylabel("累積寄与率")
ax2.set_title("累積寄与率 vs 主成分数")
ax2.legend()
ax2.grid(True, alpha=0.3)
 
plt.tight_layout()
plt.savefig("reports/pca_error_vs_components.png", dpi=150)
plt.close()

from sklearn.decomposition import IncrementalPCA
import numpy as np
 
# 大規模データをバッチ処理で PCA（メモリ節約）
n_samples = 100000
n_features = 500
batch_size = 5000
 
# データを分割して処理
ipca = IncrementalPCA(n_components=50)
 
# バッチごとに partial_fit
for i in range(0, n_samples, batch_size):
    # 実際にはファイルから読み込むなど
    X_batch = np.random.randn(min(batch_size, n_samples - i), n_features)
    ipca.partial_fit(X_batch)
 
print(f"IncrementalPCA 学習完了")
print(f"累積寄与率 (50成分): {sum(ipca.explained_variance_ratio_):.4f}")
 
# 通常の PCA との比較
# 全データをメモリに載せる場合
# from sklearn.decomposition import PCA
# pca = PCA(n_components=50)
# pca.fit(X_all)  # ← X_all が大きすぎるとメモリエラー
 
# IncrementalPCA なら OK
# データベースやファイルから chunk ごとに読み込み可能

import numpy as np
import matplotlib.pyplot as plt
from sklearn.manifold import TSNE
from sklearn.datasets import load_digits
 
digits = load_digits()
X, y = digits.data, digits.target
 
# perplexityの影響を比較
perplexities = [5, 15, 30, 50, 100]
fig, axes = plt.subplots(1, len(perplexities), figsize=(25, 5))
 
for ax, perp in zip(axes, perplexities):
    tsne = TSNE(n_components=2, perplexity=perp, random_state=42,
                n_iter=1000, learning_rate="auto", init="pca")
    X_tsne = tsne.fit_transform(X)
 
    scatter = ax.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y,
                         cmap="tab10", s=5, alpha=0.6)
    ax.set_title(f"perplexity={perp}")
    ax.set_xticks([])
    ax.set_yticks([])
 
plt.suptitle("t-SNE: perplexityの影響", fontsize=14)
plt.tight_layout()
plt.savefig("reports/tsne_perplexity.png", dpi=150)
plt.close()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.manifold import TSNE
from sklearn.datasets import load_digits
 
digits = load_digits()
X, y = digits.data, digits.target
 
# 異なる iteration 数での結果を比較
n_iters = [50, 100, 250, 500, 1000, 2000]
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
 
for ax, n_iter in zip(axes.flatten(), n_iters):
    tsne = TSNE(
        n_components=2,
        perplexity=30,
        n_iter=n_iter,
        random_state=42,
        init="pca",
        learning_rate="auto",
        method="barnes_hut",  # O(n log n) 近似
    )
    X_tsne = tsne.fit_transform(X)
 
    scatter = ax.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y,
                         cmap="tab10", s=5, alpha=0.6)
    ax.set_title(f"n_iter={n_iter}, KL={tsne.kl_divergence_:.4f}")
    ax.set_xticks([])
    ax.set_yticks([])
 
plt.suptitle("t-SNE 収束過程: イテレーション数の影響", fontsize=14)
plt.tight_layout()
plt.savefig("reports/tsne_convergence.png", dpi=150)
plt.close()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.manifold import TSNE
from sklearn.decomposition import PCA
import umap
from sklearn.datasets import load_digits
import time
 
digits = load_digits()
X, y = digits.data, digits.target
 
methods = {
    "PCA": lambda: PCA(n_components=2).fit_transform(X),
    "t-SNE": lambda: TSNE(n_components=2, random_state=42,
                           init="pca", learning_rate="auto").fit_transform(X),
    "UMAP": lambda: umap.UMAP(n_components=2, random_state=42).fit_transform(X),
}
 
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
 
for ax, (name, method) in zip(axes, methods.items()):
    start = time.time()
    X_2d = method()
    elapsed = time.time() - start
 
    scatter = ax.scatter(X_2d[:, 0], X_2d[:, 1], c=y, cmap="tab10",
                         s=5, alpha=0.6)
    ax.set_title(f"{name} ({elapsed:.2f}秒)")
    ax.set_xticks([])
    ax.set_yticks([])
 
plt.suptitle("次元削減手法の比較（手書き数字データ）", fontsize=14)
plt.colorbar(scatter, ax=axes, shrink=0.8)
plt.tight_layout()
plt.savefig("reports/dim_reduction_comparison.png", dpi=150)
plt.close()

import umap
import numpy as np
import matplotlib.pyplot as plt
 
n_neighbors_list = [5, 15, 50, 200]
min_dist_list = [0.0, 0.1, 0.5, 0.99]
 
fig, axes = plt.subplots(len(n_neighbors_list), len(min_dist_list),
                          figsize=(20, 20))
 
for i, nn in enumerate(n_neighbors_list):
    for j, md in enumerate(min_dist_list):
        reducer = umap.UMAP(n_neighbors=nn, min_dist=md,
                            n_components=2, random_state=42)
        X_umap = reducer.fit_transform(X)
 
        axes[i][j].scatter(X_umap[:, 0], X_umap[:, 1], c=y,
                           cmap="tab10", s=3, alpha=0.5)
        axes[i][j].set_title(f"nn={nn}, md={md}", fontsize=9)
        axes[i][j].set_xticks([])
        axes[i][j].set_yticks([])
 
    axes[i][0].set_ylabel(f"n_neighbors={nn}")
 
for j, md in enumerate(min_dist_list):
    axes[0][j].set_xlabel(f"min_dist={md}")
 
plt.suptitle("UMAP パラメータグリッド", fontsize=14)
plt.tight_layout()
plt.savefig("reports/umap_params.png", dpi=150)
plt.close()

import umap
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_digits
 
digits = load_digits()
X, y = digits.data, digits.target
 
# 訓練/テスト分割
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42
)
 
# UMAP を訓練データで学習
reducer = umap.UMAP(n_components=2, random_state=42)
X_train_2d = reducer.fit_transform(X_train)
 
# テストデータに transform を適用（t-SNE ではこれが不可能）
X_test_2d = reducer.transform(X_test)
 
# 可視化
import matplotlib.pyplot as plt
 
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
 
ax1.scatter(X_train_2d[:, 0], X_train_2d[:, 1], c=y_train,
            cmap="tab10", s=5, alpha=0.5)
ax1.set_title("訓練データ (fit_transform)")
ax1.set_xticks([])
ax1.set_yticks([])
 
ax2.scatter(X_test_2d[:, 0], X_test_2d[:, 1], c=y_test,
            cmap="tab10", s=5, alpha=0.5)
ax2.set_title("テストデータ (transform)")
ax2.set_xticks([])
ax2.set_yticks([])
 
plt.suptitle("UMAP: 新規データへの適用", fontsize=14)
plt.tight_layout()
plt.savefig("reports/umap_transform.png", dpi=150)
plt.close()
 
# UMAP + 分類器のパイプライン
from sklearn.pipeline import make_pipeline
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
 
pipe = make_pipeline(
    umap.UMAP(n_components=10, random_state=42),
    RandomForestClassifier(n_estimators=100, random_state=42)
)
 
scores = cross_val_score(pipe, X, y, cv=5, scoring="accuracy")
print(f"UMAP(10) + RF: Accuracy = {scores.mean():.4f} ± {scores.std():.4f}")

from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
import numpy as np
 
# 高次元データに対してPCAで前処理
from sklearn.datasets import fetch_openml
# mnist = fetch_openml("mnist_784", version=1, as_frame=False)
# X, y = mnist.data[:10000], mnist.target[:10000]
 
# サンプル: 手書き数字データ
from sklearn.datasets import load_digits
digits = load_digits()
X, y = digits.data, digits.target
 
# PCA前処理の有無で比較
configs = [
    ("RF (全次元)", make_pipeline(
        StandardScaler(),
        RandomForestClassifier(n_estimators=100, random_state=42)
    )),
    ("PCA(95%) + RF", make_pipeline(
        StandardScaler(),
        PCA(n_components=0.95),  # 95%の分散を保持
        RandomForestClassifier(n_estimators=100, random_state=42)
    )),
    ("PCA(10) + RF", make_pipeline(
        StandardScaler(),
        PCA(n_components=10),
        RandomForestClassifier(n_estimators=100, random_state=42)
    )),
]
 
import time
for name, pipe in configs:
    start = time.time()
    scores = cross_val_score(pipe, X, y, cv=5, scoring="accuracy")
    elapsed = time.time() - start
    print(f"{name:20s}  Acc={scores.mean():.4f}±{scores.std():.4f}  "
          f"時間={elapsed:.2f}秒")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA, KernelPCA
from sklearn.datasets import make_moons, make_circles
 
# 非線形データ: 半月型と同心円
datasets = {
    "半月型": make_moons(n_samples=500, noise=0.1, random_state=42),
    "同心円": make_circles(n_samples=500, noise=0.05, factor=0.3, random_state=42),
}
 
fig, axes = plt.subplots(len(datasets), 4, figsize=(20, 10))
 
for row, (name, (X, y)) in enumerate(datasets.items()):
    # 元データ
    axes[row][0].scatter(X[:, 0], X[:, 1], c=y, cmap="coolwarm", s=10)
    axes[row][0].set_title(f"{name}: 元データ")
 
    # 通常の PCA
    pca = PCA(n_components=2)
    X_pca = pca.fit_transform(X)
    axes[row][1].scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap="coolwarm", s=10)
    axes[row][1].set_title("PCA")
 
    # Kernel PCA (RBF)
    kpca_rbf = KernelPCA(n_components=2, kernel="rbf", gamma=10)
    X_kpca_rbf = kpca_rbf.fit_transform(X)
    axes[row][2].scatter(X_kpca_rbf[:, 0], X_kpca_rbf[:, 1], c=y,
                          cmap="coolwarm", s=10)
    axes[row][2].set_title("Kernel PCA (RBF)")
 
    # Kernel PCA (Polynomial)
    kpca_poly = KernelPCA(n_components=2, kernel="poly", degree=3, gamma=1)
    X_kpca_poly = kpca_poly.fit_transform(X)
    axes[row][3].scatter(X_kpca_poly[:, 0], X_kpca_poly[:, 1], c=y,
                          cmap="coolwarm", s=10)
    axes[row][3].set_title("Kernel PCA (Poly)")
 
plt.tight_layout()
plt.savefig("reports/kernel_pca_comparison.png", dpi=150)
plt.close()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import KernelPCA
from sklearn.model_selection import GridSearchCV
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import make_circles
 
X, y = make_circles(n_samples=1000, noise=0.1, factor=0.3, random_state=42)
 
# gamma のグリッドサーチ
gammas = [0.01, 0.1, 1, 5, 10, 50, 100]
 
fig, axes = plt.subplots(1, len(gammas), figsize=(28, 4))
 
for ax, gamma in zip(axes, gammas):
    kpca = KernelPCA(n_components=2, kernel="rbf", gamma=gamma)
    X_kpca = kpca.fit_transform(X)
    ax.scatter(X_kpca[:, 0], X_kpca[:, 1], c=y, cmap="coolwarm", s=5)
    ax.set_title(f"γ={gamma}")
    ax.set_xticks([])
    ax.set_yticks([])
 
plt.suptitle("Kernel PCA: gamma の影響 (RBF)", fontsize=14)
plt.tight_layout()
plt.savefig("reports/kpca_gamma.png", dpi=150)
plt.close()
 
# パイプラインでの最適 gamma 探索
pipe = make_pipeline(
    KernelPCA(kernel="rbf"),
    LogisticRegression()
)
 
param_grid = {
    "kernelpca__n_components": [2, 5, 10],
    "kernelpca__gamma": [0.01, 0.1, 1, 10],
}
 
grid = GridSearchCV(pipe, param_grid, cv=5, scoring="accuracy")
grid.fit(X, y)
 
print(f"最良パラメータ: {grid.best_params_}")
print(f"最良スコア: {grid.best_score_:.4f}")

import numpy as np
import matplotlib.pyplot as plt
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.decomposition import PCA
from sklearn.datasets import load_wine
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
 
# Wine データセット（13次元、3クラス）
wine = load_wine()
X, y = wine.data, wine.target
 
X_scaled = StandardScaler().fit_transform(X)
 
# PCA vs LDA の 2D 射影を比較
pca = PCA(n_components=2)
lda = LinearDiscriminantAnalysis(n_components=2)
 
X_pca = pca.fit_transform(X_scaled)
X_lda = lda.fit_transform(X_scaled, y)  # LDA はラベルが必要
 
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
 
for cls in np.unique(y):
    mask = y == cls
    ax1.scatter(X_pca[mask, 0], X_pca[mask, 1], label=f"Class {cls}", s=30, alpha=0.7)
    ax2.scatter(X_lda[mask, 0], X_lda[mask, 1], label=f"Class {cls}", s=30, alpha=0.7)
 
ax1.set_title(f"PCA 2D (寄与率: {sum(pca.explained_variance_ratio_):.1%})")
ax1.legend()
ax1.grid(True, alpha=0.3)
 
ax2.set_title("LDA 2D (教師あり)")
ax2.legend()
ax2.grid(True, alpha=0.3)
 
plt.suptitle("PCA vs LDA: Wine データセット", fontsize=14)
plt.tight_layout()
plt.savefig("reports/pca_vs_lda.png", dpi=150)
plt.close()
 
# 分類精度の比較
for name, reducer in [("PCA(2)", PCA(n_components=2)),
                       ("PCA(5)", PCA(n_components=5)),
                       ("LDA(2)", LinearDiscriminantAnalysis(n_components=2))]:
    pipe = make_pipeline(StandardScaler(), reducer, LogisticRegression())
    scores = cross_val_score(pipe, X, y, cv=5, scoring="accuracy")
    print(f"{name:10s}  Accuracy = {scores.mean():.4f} ± {scores.std():.4f}")

import torch
import torch.nn as nn
import torch.optim as optim
from torch.utils.data import DataLoader, TensorDataset
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from sklearn.preprocessing import StandardScaler
 
# データ準備
digits = load_digits()
X = StandardScaler().fit_transform(digits.data).astype(np.float32)
y = digits.target
 
X_tensor = torch.tensor(X)
dataset = TensorDataset(X_tensor, X_tensor)  # 入力 = 出力
loader = DataLoader(dataset, batch_size=64, shuffle=True)
 
# Autoencoder モデル
class Autoencoder(nn.Module):
    def __init__(self, input_dim=64, latent_dim=2):
        super().__init__()
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, 32),
            nn.ReLU(),
            nn.Linear(32, 16),
            nn.ReLU(),
            nn.Linear(16, latent_dim),
        )
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, 16),
            nn.ReLU(),
            nn.Linear(16, 32),
            nn.ReLU(),
            nn.Linear(32, input_dim),
        )
 
    def forward(self, x):
        z = self.encoder(x)
        x_hat = self.decoder(z)
        return x_hat, z
 
# 学習
model = Autoencoder(input_dim=64, latent_dim=2)
optimizer = optim.Adam(model.parameters(), lr=0.001)
criterion = nn.MSELoss()
 
losses = []
for epoch in range(100):
    epoch_loss = 0
    for x_batch, _ in loader:
        x_hat, z = model(x_batch)
        loss = criterion(x_hat, x_batch)
 
        optimizer.zero_grad()
        loss.backward()
        optimizer.step()
 
        epoch_loss += loss.item()
 
    losses.append(epoch_loss / len(loader))
    if (epoch + 1) % 20 == 0:
        print(f"Epoch {epoch+1}: Loss = {losses[-1]:.6f}")
 
# 潜在空間の可視化
model.eval()
with torch.no_grad():
    _, Z = model(X_tensor)
    Z = Z.numpy()
 
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 5))
 
# Autoencoder 2D
scatter = ax1.scatter(Z[:, 0], Z[:, 1], c=y, cmap="tab10", s=5, alpha=0.6)
ax1.set_title("Autoencoder 2D")
ax1.set_xticks([])
ax1.set_yticks([])
plt.colorbar(scatter, ax=ax1)
 
# PCA 2D（比較用）
from sklearn.decomposition import PCA
X_pca = PCA(n_components=2).fit_transform(X)
ax2.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap="tab10", s=5, alpha=0.6)
ax2.set_title("PCA 2D")
ax2.set_xticks([])
ax2.set_yticks([])
 
# 学習曲線
ax3.plot(losses)
ax3.set_xlabel("Epoch")
ax3.set_ylabel("Reconstruction Loss")
ax3.set_title("学習曲線")
ax3.grid(True, alpha=0.3)
 
plt.tight_layout()
plt.savefig("reports/autoencoder_vs_pca.png", dpi=150)
plt.close()

import numpy as np
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import TruncatedSVD
import umap
import matplotlib.pyplot as plt
 
# テキストデータ（サンプル）
texts = [
    "machine learning is a subset of artificial intelligence",
    "deep learning uses neural networks with many layers",
    "natural language processing deals with text data",
    "computer vision focuses on image recognition",
    "reinforcement learning learns from rewards",
    "soccer is a popular sport worldwide",
    "basketball requires good physical fitness",
    "tennis is played on different surfaces",
    "swimming is both a sport and exercise",
    "baseball has nine players on a team",
    "python is a programming language",
    "javascript runs in web browsers",
    "rust focuses on memory safety",
    "go is designed for concurrency",
    "java is used in enterprise applications",
]
categories = ["AI"] * 5 + ["Sports"] * 5 + ["Programming"] * 5
 
# TF-IDF ベクトル化
vectorizer = TfidfVectorizer(stop_words="english")
X_tfidf = vectorizer.fit_transform(texts)
print(f"TF-IDF 行列の形状: {X_tfidf.shape}")
 
# 高次元スパース → 低次元密ベクトルに変換
# Step 1: TruncatedSVD で 50 次元に削減（スパース行列対応）
svd = TruncatedSVD(n_components=min(10, X_tfidf.shape[1] - 1))
X_svd = svd.fit_transform(X_tfidf)
 
# Step 2: UMAP で 2D に可視化
reducer = umap.UMAP(n_components=2, n_neighbors=5, random_state=42)
X_2d = reducer.fit_transform(X_svd)
 
# 可視化
color_map = {"AI": "red", "Sports": "blue", "Programming": "green"}
colors = [color_map[c] for c in categories]
 
plt.figure(figsize=(10, 8))
for cat in color_map:
    mask = [c == cat for c in categories]
    plt.scatter(X_2d[np.array(mask), 0], X_2d[np.array(mask), 1],
                c=color_map[cat], label=cat, s=100)
 
for i, text in enumerate(texts):
    short_text = text[:25] + "..." if len(text) > 25 else text
    plt.annotate(short_text, (X_2d[i, 0], X_2d[i, 1]),
                 fontsize=7, alpha=0.7)
 
plt.title("テキストデータの次元削減可視化")
plt.legend()
plt.tight_layout()
plt.savefig("reports/text_dim_reduction.png", dpi=150)
plt.close()

import torch
import torchvision.models as models
import torchvision.transforms as transforms
from torch.utils.data import DataLoader
from torchvision.datasets import CIFAR10
import numpy as np
import umap
import matplotlib.pyplot as plt
 
# 事前学習済み CNN で特徴量を抽出
model = models.resnet18(pretrained=True)
model = torch.nn.Sequential(*list(model.children())[:-1])  # 最終FC層を除去
model.eval()
 
transform = transforms.Compose([
    transforms.Resize(224),
    transforms.ToTensor(),
    transforms.Normalize(mean=[0.485, 0.456, 0.406],
                         std=[0.229, 0.224, 0.225]),
])
 
# CIFAR-10 のサブセット
dataset = CIFAR10(root="./data", train=False, transform=transform, download=True)
loader = DataLoader(dataset, batch_size=128, shuffle=False)
 
# 特徴量抽出
features = []
labels = []
with torch.no_grad():
    for images, targets in loader:
        feat = model(images).squeeze()
        features.append(feat.numpy())
        labels.append(targets.numpy())
        if len(features) * 128 >= 2000:
            break
 
features = np.vstack(features)[:2000]
labels = np.concatenate(labels)[:2000]
print(f"特徴量行列: {features.shape}")  # (2000, 512)
 
# UMAP で 2D に射影
reducer = umap.UMAP(n_components=2, n_neighbors=15, random_state=42)
X_2d = reducer.fit_transform(features)
 
# 可視化
class_names = ['airplane', 'automobile', 'bird', 'cat', 'deer',
               'dog', 'frog', 'horse', 'ship', 'truck']
 
plt.figure(figsize=(12, 10))
for i, name in enumerate(class_names):
    mask = labels == i
    plt.scatter(X_2d[mask, 0], X_2d[mask, 1], label=name, s=5, alpha=0.6)
 
plt.legend(markerscale=5, loc="best")
plt.title("CIFAR-10: CNN特徴量の UMAP 可視化")
plt.xticks([])
plt.yticks([])
plt.tight_layout()
plt.savefig("reports/cnn_features_umap.png", dpi=150)
plt.close()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
 
# 正常データ + 異常データの生成
np.random.seed(42)
n_normal = 500
n_anomaly = 20
n_features = 20
 
# 正常データ: 低ランク構造（実際には5次元の構造）
true_dim = 5
W = np.random.randn(true_dim, n_features)
Z_normal = np.random.randn(n_normal, true_dim)
X_normal = Z_normal @ W + np.random.randn(n_normal, n_features) * 0.3
 
# 異常データ: 正常データの構造から外れたパターン
X_anomaly = np.random.randn(n_anomaly, n_features) * 3
 
X = np.vstack([X_normal, X_anomaly])
labels = np.array([0] * n_normal + [1] * n_anomaly)  # 0=正常, 1=異常
 
# PCA で再構成誤差を計算
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
 
pca = PCA(n_components=5)  # 正常データの真の次元数
X_pca = pca.fit_transform(X_scaled)
X_reconstructed = pca.inverse_transform(X_pca)
 
# 再構成誤差 = 異常スコア
reconstruction_error = np.mean((X_scaled - X_reconstructed) ** 2, axis=1)
 
# 閾値: 正常データの 99 パーセンタイル
threshold = np.percentile(reconstruction_error[:n_normal], 99)
predictions = (reconstruction_error > threshold).astype(int)
 
# 可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
 
# 再構成誤差の分布
ax1.hist(reconstruction_error[labels == 0], bins=30, alpha=0.5, label="正常", color="blue")
ax1.hist(reconstruction_error[labels == 1], bins=10, alpha=0.5, label="異常", color="red")
ax1.axvline(threshold, color="green", linestyle="--", label=f"閾値={threshold:.2f}")
ax1.set_xlabel("再構成誤差")
ax1.set_ylabel("頻度")
ax1.set_title("PCA 再構成誤差による異常検知")
ax1.legend()
 
# PCA 2D 可視化
X_2d = PCA(n_components=2).fit_transform(X_scaled)
ax2.scatter(X_2d[labels == 0, 0], X_2d[labels == 0, 1],
            c="blue", s=10, alpha=0.5, label="正常")
ax2.scatter(X_2d[labels == 1, 0], X_2d[labels == 1, 1],
            c="red", s=50, alpha=0.8, label="異常", marker="x")
ax2.set_title("PCA 2D 可視化")
ax2.legend()
 
plt.tight_layout()
plt.savefig("reports/pca_anomaly_detection.png", dpi=150)
plt.close()
 
# 検出精度
from sklearn.metrics import precision_score, recall_score, f1_score
print(f"Precision: {precision_score(labels, predictions):.4f}")
print(f"Recall: {recall_score(labels, predictions):.4f}")
print(f"F1: {f1_score(labels, predictions):.4f}")

import numpy as np
import time
from sklearn.decomposition import PCA, IncrementalPCA
from sklearn.manifold import TSNE
import umap
 
# ベンチマーク: データサイズ別の実行時間
data_sizes = [1000, 5000, 10000, 50000]
results = []
 
for n in data_sizes:
    X = np.random.randn(n, 100)
 
    # PCA
    start = time.time()
    PCA(n_components=2).fit_transform(X)
    pca_time = time.time() - start
 
    # t-SNE (5000以下のみ)
    if n <= 5000:
        start = time.time()
        TSNE(n_components=2, random_state=42).fit_transform(X)
        tsne_time = time.time() - start
    else:
        tsne_time = float("nan")
 
    # UMAP
    start = time.time()
    umap.UMAP(n_components=2, random_state=42).fit_transform(X)
    umap_time = time.time() - start
 
    results.append({
        "n": n,
        "PCA": f"{pca_time:.2f}s",
        "t-SNE": f"{tsne_time:.2f}s" if not np.isnan(tsne_time) else "N/A",
        "UMAP": f"{umap_time:.2f}s",
    })
    print(f"n={n:>6d}  PCA={pca_time:.2f}s  t-SNE={tsne_time:.2f}s  "
          f"UMAP={umap_time:.2f}s")

# BAD: PCA(n_components=2) で寄与率30%しかないのに結論を出す
pca = PCA(n_components=2)
X_2d = pca.fit_transform(X)
print(f"寄与率: {sum(pca.explained_variance_ratio_):.1%}")
# → 30%の情報しかない2D図で「クラスタが分離していない」と結論
 
# GOOD: まず寄与率を確認し、適切な手法を選択
pca_full = PCA()
pca_full.fit(X)
cumvar = np.cumsum(pca_full.explained_variance_ratio_)
print(f"2D寄与率: {cumvar[1]:.1%}")
 
if cumvar[1] < 0.5:
    print("PCA 2Dでは情報が不足。t-SNE or UMAPで可視化推奨")
    # t-SNE or UMAP で可視化

# BAD: 単位の異なる特徴量をそのまま PCA
# 面積(m²: 10-200) と 部屋数(1-5) → 面積の分散が圧倒的に大きい
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_raw)
# → PC1 ≈ 面積の方向（単位の影響）
 
# GOOD: StandardScaler で正規化してから PCA
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
 
pipe = make_pipeline(StandardScaler(), PCA(n_components=2))
X_pca = pipe.fit_transform(X_raw)
# → 全特徴量が等しくスケーリングされた状態で PCA

# BAD: 全データで次元削減 → 分割 → モデル評価
# → テストデータの情報が次元削減に漏洩（データリーク）
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
 
pca = PCA(n_components=10)
X_reduced = pca.fit_transform(X)  # 全データで fit（リーク!）
X_train, X_test = train_test_split(X_reduced)
 
# GOOD: Pipeline 内で次元削減 → CV 内で fit/transform が分離される
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import cross_val_score
 
pipe = make_pipeline(PCA(n_components=10), LogisticRegression())
scores = cross_val_score(pipe, X, y, cv=5)
# → 各 fold で PCA の fit は訓練データのみに適用される