Rの関数:NLSstAsymptotic {stats}

# パッケージの読み込み
library(ggplot2)

# 乱数シードの固定
seed <- 20260217
set.seed(seed)

# NLSstAsymptotic 関数による初期値推定のシミュレーション

# 1. サンプルデータの生成
# モデル: y = b0 + b1 * (1 - exp(-exp(lrc) * x))
# 化学反応の進行度や、生物の成長曲線などを模したデータを作成します。

# 真のパラメータ設定
true_b0 <- 10 # x=0 での値（初期値）
true_b1 <- 20 # 総変化量（最終的には 10 + 20 = 30 に収束する）
true_lrc <- log(0.5) # レート定数 k = 0.5, lrc = -0.693

# xデータの作成 (0から10まで)
x_vals <- seq(0, 10, length.out = 30)

# 真のy値の計算
y_true <- true_b0 + true_b1 * (1 - exp(-exp(true_lrc) * x_vals))

# 観測データ (ノイズを付加)
y_obs <- y_true + rnorm(length(x_vals), mean = 0, sd = 1.0)

# データフレーム化
df_sim <- data.frame(x = x_vals, y = y_obs)

# 2. sortedXyData オブジェクトへの変換
# NLSstAsymptotic はこのクラスを引数に取るため、変換が必要です。
xy_data <- sortedXyData(df_sim$x, df_sim$y)

cat("--- データ概要 ---\n")
cat(sprintf("真の b0 (切片): %.2f\n", true_b0))
cat(sprintf("真の b1 (変化幅): %.2f (漸近線 = %.2f)\n", true_b1, true_b0 + true_b1))
cat(sprintf("真の lrc (対数レート): %.4f (k = 0.5)\n\n", true_lrc))
cat("\n--- データの一部を確認 ---\n")
print(str(xy_data))

# 3. NLSstAsymptotic による初期値の推定
init_vals <- NLSstAsymptotic(xy_data)

cat("--- NLSstAsymptotic による推定初期値 ---\n")
print(init_vals)

cat("\n--- 真の値との比較 ---\n")
cat(sprintf("b0 (推定): %.4f vs (真): %.4f\n", init_vals["b0"], true_b0))
cat(sprintf("b1 (推定): %.4f vs (真): %.4f\n", init_vals["b1"], true_b1))
cat(sprintf("lrc(推定): %.4f vs (真): %.4f\n\n", init_vals["lrc"], true_lrc))

# 4. 推定された初期値を用いた nls の実行
# NLSstAsymptotic が返した値を start 引数に使ってフィッティングを行います。
fit_nls <- nls(
  y ~ b0 + b1 * (1 - exp(-exp(lrc) * x)),
  data = df_sim,
  start = list(b0 = init_vals["b0"], b1 = init_vals["b1"], lrc = init_vals["lrc"])
)

cat("--- nls による最終的なモデル適合結果 ---\n")
print(summary(fit_nls))