［R］例2.2 風速と直流発電量（「線形回帰分析」（朝倉書店）pp.62-67）その2

pp.64-67に掲載されている、p.23の式(2.30)による結果を流れに沿って計算している。表2.1（p.62）の値をCSV形式で入力したtable2_1.csvをカレントディレクトリに置いておくこと。

> dtf <- read.csv("table2_1.csv", header = TRUE)
> xxi <- as.double(1 /dtf$WV)
> yyi <- as.double(dtf$DC)
> n <- nrow(dtf)
> k <- 2
> ssxx <- sum(xxi)
> ssyy <- sum(yyi)
> ssxx2 <- sum(xxi ^ 2)
> ssyy2 <- sum(yyi ^ 2)
> ssxxyy <- sum(xxi * yyi)
> xxm <- mean(xxi)
> yym <- mean(yyi)
> ssx2 <- ssxx2 - ssxx ^ 2 / n
> ssy2 <- ssyy2 - ssyy ^ 2 / n
> ssxy <- ssxxyy - ssxx * ssyy / n
> bh <- ssxy / ssx2
> ah <- yym - bh * xxm
> ssyh2 <- bh * ssxy
> r2 <- ssyh2 / ssy2
> sse2 <- ssy2 - ssyh2
> s2 <- sse2 / (n - k)
> s <- sqrt(s2)
> sah2 <- s2 * (1 / n + xxm ^ 2 / ssx2)
> sah <- sqrt(sah2)
> ta <- ah / sah
> sbh2 <- s2 / ssx2
> sbh <- sqrt(sbh2)
> tb <- bh / sbh
> yyhi <- ah + bh * xxi
> ei <- yyi - yyhi
> ri <- 100 * ei / yyi
> ai2 <- 100 * (ei ^ 2 / sse2)
> sh2 <- sum(ei ^ 2) / n
> print(ah)  # 推定されたパラメーター α^
[1] 2.97886
> print(bh)  # 推定されたパラメーター β^
[1] -6.934547
> print(sh2)  # 誤差項の分散の最尤推定量 σ^^2
[1] 0.008158786
> dtf <- data.frame(yyi, yyhi, ei, ri, ai2)
> colnames(dtf) <- c("DC", "DC推定値", "残差", "誤差率(%)", "平方残差率(%)")
> print(dtf)
      DC  DC推定値         残差    誤差率(%) 平方残差率(%)
1  1.582 1.5919507 -0.009950703  -0.62899514   0.048544719
2  1.822 1.8231023 -0.001102282  -0.06049844   0.000595689
3  1.057 0.9392874  0.117712578  11.13647849   6.793290635
4  0.500 0.4105093  0.089490699  17.89813971   3.926361222
5  2.236 2.2854054 -0.049405439  -2.20954556   1.196696381
6  2.386 2.2639584  0.122041615   5.11490423   7.302143226
7  2.294 2.2527296  0.041270439   1.79906010   0.835050285
8  0.558 0.7052381 -0.147238090 -26.38675451  10.628569260
9  2.166 2.1279955  0.038004532   1.75459519   0.708117342
10 1.866 1.8603848  0.005615206   0.30092206   0.015458445
11 0.653 0.5876369  0.065363052  10.00965576   2.094590369
12 1.930 1.8868055  0.043194527   2.23805841   0.914727865
13 1.562 1.4713499  0.090650121   5.80346482   4.028758432
14 1.737 1.7832486 -0.046248561  -2.66255390   1.048650839
15 2.088 2.0417592  0.046240820   2.21459864   1.048299802
16 1.137 1.0525970  0.084402980   7.42330521   3.492609407
17 2.179 2.0954783  0.083521654   3.83302682   3.420051394
18 2.112 2.1908434 -0.078843429  -3.73311692   3.047652653
19 1.800 1.9882106 -0.188210552 -10.45614178  17.366903572
20 1.501 1.7064662 -0.205466164 -13.68861852  20.697366094
21 2.303 2.2168220  0.086177997   3.74198856   3.641055032
22 2.310 2.2990026  0.010997410   0.47607833   0.059294615
23 1.194 1.2875072 -0.093507161  -7.83142048   4.286710893
24 1.144 1.2232786 -0.079278565  -6.92994450   3.081385381
25 0.123 0.1484327 -0.025432682 -20.67697705   0.317116449
> # 図2.8と2.9の作成
> png("fig2_8.png", width = 512, height = 512)
> plot(1 / xxi, yyi, xlim = c(2, 11), ylim = c(0, 2.5), xlab = "WV", ylab = "DC", pch = 20)
> idx <- order(1 / xxi)
> lines(1 / xxi[idx], yyhi[idx])
> dev.off()
null device 
          1 
> png("fig2_9.png", width = 512, height = 512)
> plot(xxi, yyi, xlim = c(0.09, 0.42), ylim = c(0, 2.5), xlab = "1 / WV", ylab = "DC", pch = 20)
> idx <- order(xxi)
> lines(xxi[idx], yyhi[idx])
> dev.off()
null device 
          1