［R］重み付き非線形最小二乗法を行う

nls関数を使う。weightsオプションに重みを指定すると、その重みを考慮した計算を行う。

以下の例では、与えられた8組の測定値を使用して、非線形最小二乗法により回帰モデルy = b1 + b2 * x * 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] * x * exp(b[3] * x)
> b0 <- c(1, 1, 1)
> 
> # 重みを考慮しない計算（図の青実線）
> r1 <- nls(y ~ f(x, b), start = list(b = b0))
> sm <- summary(r1)
> print(sm)       # nls関数による計算結果の概要
Formula: y ~ f(x, b)
Parameters:
   Estimate Std. Error t value Pr(>|t|)   
b1 -3.89078    2.46437  -1.579  0.17521   
b2  3.39494    0.74297   4.569  0.00601 **
b3  0.10847    0.02318   4.679  0.00544 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.502 on 5 degrees of freedom
Number of iterations to convergence: 21 
Achieved convergence tolerance: 3.21e-06
> print(sm$coef)  # 推定した回帰係数,その標準誤差,t値,p値
     Estimate Std. Error   t value    Pr(>|t|)
b1 -3.8907795 2.46436971 -1.578813 0.175211732
b2  3.3949412 0.74297094  4.569413 0.006005113
b3  0.1084657 0.02318143  4.678990 0.005438508
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 31.30051
> print(sum((y - fitted(r1)) ^ 2))
[1] 31.30051
> print(sm$df)                           # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 6.260101
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 2.502019
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            b1          b2            b3
b1  6.07311807 -1.56568510  0.0449342593
b2 -1.56568510  0.55200581 -0.0170136719
b3  0.04493426 -0.01701367  0.0005373789
> 
> # 重みを考慮した計算（図の赤実線）
> r2 <- nls(y ~ f(x, b), weights = w, start = list(b = b0))
> sm <- summary(r2)
> print(sm)       # nls関数による計算結果の概要
Formula: y ~ f(x, b)
Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b1 -0.03045    1.18131  -0.026 0.980435    
b2  2.06079    0.33571   6.139 0.001666 ** 
b3  0.15929    0.02080   7.659 0.000604 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.741 on 5 degrees of freedom
Number of iterations to convergence: 14 
Achieved convergence tolerance: 3.676e-06
> print(sm$coef)  # 推定した回帰係数,その標準誤差,t値,p値
      Estimate Std. Error     t value     Pr(>|t|)
b1 -0.03044686 1.18131366 -0.02577373 0.9804348365
b2  2.06079018 0.33571296  6.13854825 0.0016664488
b3  0.15929178 0.02079896  7.65864353 0.0006043864
> print(sum(sm$resi ^ 2))                # 残差平方和
[1] 69.95679
> print(sum(w * (y - fitted(r2)) ^ 2))
[1] 69.95679
> print(sm$df)                           # 推定した回帰係数の数と自由度
[1] 3 5
> print(sum(sm$resi ^ 2) / sm$df[2])     # 残差分散
[1] 13.99136
> print(sm$sigma)                        # 回帰値の標準誤差
[1] 3.740502
> print(sm$sigma ^ 2 * sm$cov.unscaled)  # 推定した回帰係数の分散共分散行列
            b1           b2            b3
b1  1.39550197 -0.371789063  0.0219612500
b2 -0.37178906  0.112703189 -0.0069215507
b3  0.02196125 -0.006921551  0.0004325965
> 
> # 図の描画
> plot(x, fitted(r1), type = "l", col = "blue", xlab = "x", ylab = "y")
> lines(x, fitted(r2), col = "red")
> points(x, y, pch = 20)