［R］重み付き非線形最小二乗法を行う
nls関数を使う。weightsオプションに重みを指定すると、その重みを考慮した計算を行う。
以下の例では、与えられた8組の測定値を使用して、非線形最小二乗法により回帰モデルy = b1 + b2 * exp(b3 * x)（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。
> # 計算の用意
> x <- c(1, 2, 3, 5, 6, 7, 8, 9)
> y <- c(1, 4, 10, 25, 32, 49, 64, 75) # 測定値（図の黒丸）
> w <- rep(1, length(x))
> w[c(3, 5)] <- 100
> dtf <- data.frame(x, y, w)
> print(dtf)
x y w
1 1 1 1
2 2 4 1
3 3 10 100
4 5 25 1
5 6 32 100
6 7 49 1
7 8 64 1
8 9 75 1
> f <- function(x, b) b[1] + b[2] * exp(b[3] * x)
> b0 <- c(1, 1, 1)
>
> # 重みを考慮しない計算（図の青実線）
> r1 <- nls(y ~ f(x, b), start = list(b = b0))
> sm <- summary(r1)
> print(r1) # nls関数による計算結果の概要
Nonlinear regression model
model: y ~ f(x, b)
data: parent.frame()
b1 b2 b3
-25.4828 21.1533 0.1756
residual sum-of-squares: 32.98
Number of iterations to convergence: 21
Achieved convergence tolerance: 3.185e-06
> print(sm$coef) # 推定した回帰係数,その標準誤差,t値,p値
Estimate Std. Error t value Pr(>|t|)
b1 -25.4827638 10.22340189 -2.492591 0.054986347
b2 21.1532994 8.08590380 2.616071 0.047321602
b3 0.1755871 0.03299953 5.320898 0.003137863
> print(sum(sm$resi ^ 2)) # 残差平方和
[1] 32.98248
> print(sum((y - fitted(r1)) ^ 2))
[1] 32.98248
> print(sm$df) # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2]) # 残差分散
[1] 6.596497
> print(sm$sigma) # 回帰値の標準誤差
[1] 2.568365
> print(sm$sigma ^ 2 * sm$cov.unscaled) # 推定した回帰係数の分散共分散行列
b1 b2 b3
b1 104.5179462 -81.8524151 0.329318514
b2 -81.8524151 65.3818402 -0.265714306
b3 0.3293185 -0.2657143 0.001088969
>
> # 重みを考慮した計算（図の赤実線）
> r2 <- nls(y ~ f(x, b), weights = w, start = list(b = b0))
> sm <- summary(r2)
> print(r2) # nls関数による計算結果の概要
Nonlinear regression model
model: y ~ f(x, b)
data: parent.frame()
b1 b2 b3
-10.9272 10.1219 0.2413
weighted residual sum-of-squares: 70.58
Number of iterations to convergence: 16
Achieved convergence tolerance: 5.171e-07
> print(sm$coef) # 推定した回帰係数,その標準誤差,t値,p値
Estimate Std. Error t value Pr(>|t|)
b1 -10.927232 3.53132490 -3.094372 0.0270250777
b2 10.121862 2.44947466 4.132258 0.0090650345
b3 0.241281 0.02632775 9.164514 0.0002593529
> print(sum(sm$resi ^ 2)) # 残差平方和
[1] 70.57934
> print(sum(w * (y - fitted(r2)) ^ 2))
[1] 70.57934
> print(sm$df) # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2]) # 残差分散
[1] 14.11587
> print(sm$sigma) # 回帰値の標準誤差
[1] 3.757109
> print(sm$sigma ^ 2 * sm$cov.unscaled) # 推定した回帰係数の分散共分散行列
b1 b2 b3
b1 12.47025554 -8.57361848 0.0908751718
b2 -8.57361848 5.99992612 -0.0642324913
b3 0.09087517 -0.06423249 0.0006931503
>
> # 図の描画
> plot(x, fitted(r1), type = "l", col = "blue", xlab = "x", ylab = "y")
> lines(x, fitted(r2), col = "red")
> points(x, y, pch = 20)
