［R］例2.3 直流発電量の予測と予測区間（「線形回帰分析」（朝倉書店）pp.72-75）

> dtf <- read.csv("table2_1.csv", header = TRUE)
> n <- nrow(dtf)
> xxi <- 1 / dtf$WV
> yyi <- dtf$DC
> k <- 2
> degf <- n - k
> xxm <- mean(xxi)
> yym <- mean(yyi)
> ssx2 <- sum((xxi - xxm) ^ 2)
> ssxy <- sum((xxi - xxm) * (yyi - yym))
> bh <- ssxy / ssx2
> ah <- mean(yyi) - bh * xxm
> sse2 <- sum((yyi - (ah + bh * xxi)) ^ 2)
> s <- sqrt(sse2 / (n - k))
> xx0i <- c(xxi, 1 / 20)
> yy0i <- c(yyi, 2.632)
> yyh0i <- ah + bh * xx0i
> #
> tl2 <- qt(0.025, degf, lower.tail = FALSE)
> yycon <- tl2 * s * sqrt(1 / n + (xx0i - xxm) ^ 2 / ssx2)
> yypre <- tl2 * s * sqrt(1 + 1 / n + (xx0i - xxm) ^ 2 / ssx2)
> dtf <- data.frame(yy0i, yyh0i, yyh0i - yycon, yyh0i + yycon, yyh0i - yypre, yyh0i + yypre)
> names(dtf) <- c("DC", "DCの予測値", "E(DC)予測下限", "E(DC)予測上限", "DC予測下限", "DC予測上限")
> print(dtf)
      DC DCの予測値 E(DC)予測下限 E(DC)予測上限  DC予測下限 DC予測上限
1  1.582  1.5919507    1.55297389     1.6309275  1.39328146  1.7906199
2  1.822  1.8231023    1.78198202     1.8642225  1.62400142  2.0222031
3  1.057  0.9392874    0.88252512     0.9960497  0.73637800  1.1421968
4  0.500  0.4105093    0.32701897     0.4939996  0.19856376  0.6224548
5  2.236  2.2854054    2.22839667     2.3424142  2.08242693  2.4883839
6  2.386  2.2639584    2.20790651     2.3200103  2.06124655  2.4666702
7  2.294  2.2527296    2.19717273     2.3082864  2.05015405  2.4553051
8  0.558  0.7052381    0.63727035     0.7732058  0.49891337  0.9115628
9  2.166  2.1279955    2.07762556     2.1783654  1.92678065  2.3292103
10 1.866  1.8603848    1.81847394     1.9022956  1.66111915  2.0596504
11 0.653  0.5876369    0.51361863     0.6616553  0.37924072  0.7960332
12 1.930  1.8868055    1.84426818     1.9293428  1.68740713  2.0862038
13 1.562  1.4713499    1.43146889     1.5112309  1.27250127  1.6701985
14 1.737  1.7832486    1.74284605     1.8236511  1.58429470  1.9822024
15 2.088  2.0417592    1.99457587     2.0889425  1.84131831  2.2422000
16 1.137  1.0525970    1.00068770     1.1045063  0.85099133  1.2542027
17 2.179  2.0954783    2.04635313     2.1446036  1.89457150  2.2963852
18 2.112  2.1908434    2.13793584     2.2437510  1.98897840  2.3927085
19 1.800  1.9882106    1.94280548     2.0336156  1.78818081  2.1882403
20 1.501  1.7064662    1.66705050     1.7458818  1.50771036  1.9052220
21 2.303  2.2168220    2.16281924     2.2708248  2.01466717  2.4189768
22 2.310  2.2990026    2.24137972     2.3566255  2.09585075  2.5021544
23 1.194  1.2875072    1.24378717     1.3312272  1.08785318  1.4871611
24 1.144  1.2232786    1.17762784     1.2689293  1.02319292  1.4233642
25 0.123  0.1484327    0.05037867     0.2464867 -0.06966102  0.3665264
26 2.632  2.6321328    2.55808466     2.7061810  2.42372598  2.8405396
> # 図2.13
> idx <- order(xx0i)
> png("fig2_13.png", width = 512, height = 512)
> plot(xx0i, yy0i, type = "n")
> lines(xx0i[idx], yyh0i[idx] - yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] + yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] - yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx] + yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx], lty = "solid", col = "gray", lwd = 4)
> points(xx0i, yy0i, pch = 20)
> dev.off()
null device 
          1 
> #
> tl2 <- qt(0.005, degf, lower.tail = FALSE)
> yycon <- tl2 * s * sqrt(1 / n + (xx0i - xxm) ^ 2 / ssx2)
> yypre <- tl2 * s * sqrt(1 + 1 / n + (xx0i - xxm) ^ 2 / ssx2)
> dtf <- data.frame(yy0i, yyh0i, yyh0i - yycon, yyh0i + yycon, yyh0i - yypre, yyh0i + yypre)
> names(dtf) <- c("DC", "DCの予測値", "E(DC)予測下限", "E(DC)予測上限", "DC予測下限", "DC予測上限")
> print(dtf)
      DC DCの予測値 E(DC)予測下限 E(DC)予測上限 DC予測下限 DC予測上限
1  1.582  1.5919507    1.53905601     1.6448454  1.3223405  1.8615609
2  1.822  1.8231023    1.76729876     1.8789058  1.5529063  2.0932982
3  1.057  0.9392874    0.86225639     1.0163185  0.6639229  1.2146519
4  0.500  0.4105093    0.29720616     0.5238124  0.1228821  0.6981365
5  2.236  2.2854054    2.20803993     2.3627709  2.0099472  2.5608637
6  2.386  2.2639584    2.18789146     2.3400253  1.9888620  2.5390547
7  2.294  2.2527296    2.17733445     2.3281247  1.9778182  2.5276409
8  0.558  0.7052381    0.61300037     0.7974758  0.4252388  0.9852374
9  2.166  2.1279955    2.05963943     2.1963515  1.8549307  2.4010603
10 1.866  1.8603848    1.80350838     1.9172612  1.5899652  2.1308044
11 0.653  0.5876369    0.48718811     0.6880858  0.3048264  0.8704475
12 1.930  1.8868055    1.82907892     1.9445320  1.6162058  2.1574051
13 1.562  1.4713499    1.41722815     1.5254716  1.2014962  1.7412035
14 1.737  1.7832486    1.72841908     1.8380780  1.5132521  2.0532450
15 2.088  2.0417592    1.97772761     2.1057908  1.7697447  2.3137736
16 1.137  1.0525970    0.98215187     1.1230422  0.7790018  1.3261922
17 2.179  2.0954783    2.02881146     2.1621452  1.8228315  2.3681252
18 2.112  2.1908434    2.11904355     2.2626433  1.9168963  2.4647906
19 1.800  1.9882106    1.92659220     2.0498289  1.7167540  2.2596671
20 1.501  1.7064662    1.65297592     1.7599564  1.4367385  1.9761939
21 2.303  2.2168220    2.14353588     2.2901081  1.9424816  2.4911625
22 2.310  2.2990026    2.22080369     2.3772015  2.0233091  2.5746961
23 1.194  1.2875072    1.22817560     1.3468387  1.0165606  1.5584538
24 1.144  1.2232786    1.16132684     1.2852303  0.9517462  1.4948110
25 0.123  0.1484327    0.01536546     0.2814999 -0.1475381  0.4444035
26 2.632  2.6321328    2.53164348     2.7326221  2.3493079  2.9149577
> # 図2.14
> png("fig2_14.png", width = 512, height = 512)
> plot(xx0i, yy0i, type = "n")
> lines(xx0i[idx], yyh0i[idx] - yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] + yycon[idx], lty = "dotted")
> lines(xx0i[idx], yyh0i[idx] - yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx] + yypre[idx], lty = "dashed")
> lines(xx0i[idx], yyh0i[idx], lty = "solid", col = "gray", lwd = 4)
> points(xx0i, yy0i, pch = 20)
> dev.off()
null device 
          1