Rで線形代数:冪等行列

# --------------------------------------------------
# 冪等行列のシミュレーション
# --------------------------------------------------

# 1. データの準備（デザイン行列 X の作成）
seed <- 20250623
set.seed(seed)

# 観測数 n=10, 説明変数の数 p=3 のランダムな行列Xを作成
X <- matrix(rnorm(10 * 3), nrow = 10, ncol = 3)

# 2. ハット行列 H の計算
# H = X %*% solve(t(X) %*% X) %*% t(X)
#   t(X): Xの転置
#   solve(): 逆行列を計算
#   %*%: 行列の積
XtX_inv <- solve(t(X) %*% X)
H <- X %*% XtX_inv %*% t(X)

cat("ハット行列 H (最初の5x5部分) を表示:\n")
print(round(H[1:5, 1:5], 4))


# 3. 冪等性の確認 (H^2 = H)
cat("\n--- 1. 冪等性の確認 (H^2 = H) ---\n")

# Hの2乗を計算
H2 <- H %*% H

# H と H2 が等しいかを確認
check_idempotent <- all.equal(H, H2)

if (isTRUE(check_idempotent)) {
  cat("確認結果: H^2 = H が成り立ちます。Hは冪等行列です。\n")
} else {
  cat("確認結果: H^2 = H は成り立ちませんでした。\n")
  print(check_idempotent) # 差分を表示
}


# 4. 性質の確認
cat("\n--- 2. 性質の確認 ---\n")

# (a) 固有値が 0 または 1 であることの確認
eigen_values <- eigen(H)$values

cat("(a) 固有値の確認:\n")
# 計算誤差を考慮して、小数点以下6桁で丸めて表示
print(round(eigen_values, 6))
cat("→ 固有値がすべて 0 または 1 になっていることがわかります。\n\n")

# (b) rank(H) = tr(H) であることの確認
# ランクの計算
rank_H <- qr(H)$rank

# トレースの計算 (対角成分の和)
trace_H <- sum(diag(H))

cat("(b) ランクとトレースの確認:\n")
cat("ランク (rank):", rank_H, "\n")
cat("トレース (trace):", trace_H, "\n")
cat("→ ランクとトレースが一致していることがわかります。\n")
# 理論上、ハット行列のランクとトレースは説明変数の数(p=3)に等しくなります。


# (c) I - H も冪等行列であることの確認
cat("\n--- 3. (おまけ) I - H の冪等性の確認 ---\n")

# Hと同じサイズの単位行列Iを作成
I <- diag(nrow(H))
ImH <- I - H

# (I-H)の2乗を計算
ImH2 <- ImH %*% ImH

check_ImH <- all.equal(ImH, ImH2)

if (isTRUE(check_ImH)) {
  cat("確認結果: (I - H)^2 = (I - H) が成り立ちます。I - H も冪等行列です。\n")
} else {
  cat("確認結果: (I - H)^2 = (I - H) は成り立ちませんでした。\n")
}