Rで線形代数:レーベンバーグ･マーカート法

# --------------------------------------------------
# レーベンバーグ･マーカート法のシミュレーション
# --------------------------------------------------

library(minpack.lm)

# --- 1. データの生成 ---

# 真のパラメータを設定
A_true <- 100
k_true <- 0.5

# モデル関数を定義
model_func <- function(x, A, k) {
  A * exp(-k * x)
}

# xデータを生成
x_data <- seq(0, 5, length.out = 50)

# 真のモデルからyデータを生成し、ノイズを加える
seed <- 20250625
set.seed(seed)
y_noise <- rnorm(length(x_data), mean = 0, sd = 5)
y_data <- model_func(x_data, A_true, k_true) + y_noise

# 生成したデータをプロットして確認
# plot(x_data, y_data, main = "Generated Data with Noise", xlab = "Time (x)", ylab = "Value (y)", pch = 16)
# curve(model_func(x, A_true, k_true), from=0, to=5, col="red", lwd=2, add=TRUE)
# legend("topright", legend=c("Noisy Data", "True Model"), col=c("black", "red"), pch=c(16, NA), lty=c(NA, 1))


# --- 2. LM法による非線形最小二乗フィッティング ---

# (a) 残差関数を定義
#   nls.lm()は、パラメータベクトル par を受け取り、
#   残差ベクトル（実測値 - 予測値）を返す関数を要求する
residual_func <- function(par, x, y_obs) {
  A <- par[1]
  k <- par[2]

  y_pred <- model_func(x, A, k)

  # 残差を返す
  residuals <- y_obs - y_pred
  return(residuals)
}

# (b) パラメータの初期値を設定
# 意図的に離れた初期値から開始
par_start <- c(A = 80, k = 0.3)
cat("初期パラメータ GUESS: A =", par_start[1], ", k =", par_start[2], "\n")

# (c) nls.lm() を実行してLM法で最適化
#    controlオプションで詳細なトレース情報を表示できる
fit_lm <- nls.lm(
  par = par_start,
  fn = residual_func,
  x = x_data,
  y_obs = y_data,
  control = nls.lm.control(nprint = 1)
) # nprint=1で各ステップの情報を表示

# --- 3. 結果の確認 ---

# (a) 最適化結果のサマリーを表示
cat("\n--- Optimization Summary ---\n")
summary(fit_lm)

# (b) 推定されたパラメータを取得
par_estimated <- coef(fit_lm)
cat("\n推定されたパラメータ (A, k):\n")
print(par_estimated)
cat("\n真のパラメータ (A, k):\n")
print(c(A_true, k_true))

# (c) 結果をプロットして比較
plot(x_data, y_data,
  main = "Levenberg-Marquardt Fitting Result",
  pch = 16, col = "grey"
)
# 推定されたパラメータで曲線を描画
curve(model_func(x, par_estimated[1], par_estimated[2]),
  from = 0, to = 5, col = "blue", lwd = 2, add = TRUE
)
# 真のモデルも点線で描画
curve(model_func(x, A_true, k_true),
  from = 0, to = 5, col = "red", lty = 2, lwd = 2, add = TRUE
)

legend("topright",
  legend = c("Noisy Data", "Fitted Model (LM)", "True Model"),
  col = c("grey", "blue", "red"),
  pch = c(16, NA, NA),
  lty = c(NA, 1, 2),
  lwd = 2
)