Rで線形代数:偏差行列、分散共分散行列、相関行列

# --------------------------------------------------
# 偏差行列、分散共分散行列、相関行列のシミュレーション
# --------------------------------------------------

# --- 準備: データ行列の作成 ---
# サンプル数 n=5, 変数の数 p=3
# 5人の生徒の「数学」「理科」「英語」の点数をイメージ
seed <- 20250624
set.seed(seed)
X <- matrix(round(runif(15, min = 50, max = 100)), nrow = 5, ncol = 3)
colnames(X) <- c("Math", "Science", "English")
rownames(X) <- paste0("Student", 1:5)

cat("--- 0. 元のデータ行列 X ---\n")
print(X)

# --- 1. 偏差行列 D の計算 (中心化) ---
cat("\n--- 1. 偏差行列 D の計算 ---\n")

# (a) 各列（変数）の平均を計算
col_means <- colMeans(X)
cat("各変数の平均:\n")
print(col_means)

# (b) Rのsweep()関数を使って、各要素から列平均を引く
# sweep(X, 2, STATS, FUN) は、行列Xの各行/列(2は列)にSTATSを適用する
D <- sweep(X, 2, col_means, FUN = "-")
cat("\n偏差行列 D:\n")
print(round(D, 2))

#
cat("\n(参考) Rのscale()関数も利用できます:\n")
scale(X, center = TRUE, scale = FALSE)


# --- 2. 分散共分散行列 Σ の計算 ---
cat("\n--- 2. 分散共分散行列 Σ の計算 ---\n")
n <- nrow(X) # サンプル数

# (a) D^T * D を計算（偏差平方和・積和行列）
DtD <- t(D) %*% D
cat("偏差平方和・積和行列 (D^T * D):\n")
print(round(DtD, 2))

# (b) サンプル数-1 (n-1)で割る
Sigma_manual <- (1 / (n - 1)) * DtD
cat("\n手計算による分散共分散行列 Σ:\n")
print(round(Sigma_manual, 2))

# (c) Rの組み込み関数 cov() の結果と比較
Sigma_builtin <- cov(X)
cat("\nRのcov()関数による分散共分散行列:\n")
print(round(Sigma_builtin, 2))
cat("→ 手計算の結果と一致しました。\n")


# --- 3. 相関行列 R の計算 ---
cat("\n--- 3. 相関行列 R の計算 ---\n")

# (a) 各変数の標準偏差を計算 (分散の平方根)
sds <- sqrt(diag(Sigma_manual))
cat("各変数の標準偏差:\n")
print(sds)

# (b) 標準偏差の逆数からなる対角行列を作成
# D_inv_sd = diag(1/s)
D_inv_sd <- diag(1 / sds)

# (c) R = D_inv_sd * Σ * D_inv_sd で計算
R_manual <- D_inv_sd %*% Sigma_manual %*% D_inv_sd
cat("\n手計算による相関行列 R:\n")
print(round(R_manual, 2))

# (d) Rの組み込み関数 cor() の結果と比較
R_builtin <- cor(X)
cat("\nRのcor()関数による相関行列:\n")
print(round(R_builtin, 2))
cat("→ 手計算の結果と一致しました。\n")