Rで線形代数:一般逆行列(擬似逆行列)

# --------------------------------------------------
# 一般逆行列のシミュレーション
# --------------------------------------------------

library(MASS)

# --- シナリオ1: 正則行列の場合 ---
# 一般逆行列が、通常の逆行列と一致することを確認

cat("--- シナリオ1: 正則行列の場合 ---\n")
A <- matrix(c(2, 1, 3, 4), nrow = 2)
cat("正則行列 A:\n")
print(A)

# 通常の逆行列
A_inv <- solve(A)
cat("\n通常の逆行列 (solve(A)):\n")
print(A_inv)

# 一般逆行列
A_ginv <- ginv(A)
cat("\n一般逆行列 (ginv(A)):\n")
print(A_ginv)

# 両者が一致するか確認
cat("\n両者は一致するか？ ->", all.equal(A_inv, A_ginv), "\n\n")


# --- シナリオ2: 長方形行列（回帰分析での応用） ---
# 最小二乗解を求める例

cat("--- シナリオ2: 長方形行列（回帰分析での応用） ---\n")
seed <- 20250623
set.seed(seed)

# y = 5 + 2x + (ノイズ) という関係のデータを作成
x <- 1:20
y <- 5 + 2 * x + rnorm(20, mean = 0, sd = 3)

# 線形回帰モデル y = β0 + β1*x は、行列で Y = Xβ と書ける
# X は「デザイン行列」
X <- cbind(1, x) # 第1列は切片(β0)用、第2列はx(β1)用
cat("デザイン行列 X (最初の6行のみ):\n")
print(head(X))
cat("目的変数 y (最初の6個のみ):\n")
print(head(y))

# Xは 20x2 の長方形行列なので、通常の逆行列は持たない
# 一般逆行列を使って係数βの最小二乗解を求める
# β_hat = (X^T X)^-1 X^T y と同じ結果になる
beta_hat_ginv <- ginv(X) %*% y
cat("\nginv()で求めた回帰係数 (β0, β1):\n")
print(t(beta_hat_ginv))

# Rの標準的な回帰分析関数 lm() の結果と比較
model <- lm(y ~ x)
cat("\nlm()で求めた回帰係数:\n")
print(coef(model))

cat("\n→ 長方形行列の一般逆行列を使うことで、回帰分析の最小二乗解が求められることがわかります。\n\n")


# --- シナリオ3: 特異行列の場合 ---
# 通常の逆行列は計算できないが、一般逆行列は計算できる例

cat("--- シナリオ3: 特異行列の場合 ---\n")

# 2行目が1行目の2倍になっている特異行列
B <- matrix(c(1, 2, 2, 4), nrow = 2)
cat("特異行列 B:\n")
print(B)

# solve(B)はエラーになる
cat("\nsolve(B) を実行するとエラーになります。\n")
result <- tryCatch(
  {
    solve(B)
  },
  error = function(e) {
    return(e) # エラーオブジェクトを返す
  }
)
print(result$message)

# ginv(B)は計算可能
B_ginv <- ginv(B)
cat("\n一般逆行列 ginv(B):\n")
print(B_ginv)

# 一般逆行列の性質 B * B+ * B = B を確認
B_check <- B %*% B_ginv %*% B
cat("\n性質の確認 (B %*% ginv(B) %*% B):\n")
print(round(B_check, 10))
cat("→ 元の行列 B に戻りました。\n")