MINAKA Nobuhiro — Jun 26, 2014, 3:32 PM
# ーーー【ランダム誤差と回帰】ーーー
#【0】直線「y = x」に正規分布誤差を付加する
# 区間 [-1, +1] からランダムに100個の変数値を選ぶ（定義域）
n <- 100
x <- runif(n, min = -1, max = 1)
# ランダム誤差：正規分布 N(0, 0.1)
y <- x + rnorm(n, 0, 0.1)
plot(x, y)
# 回帰直線とそのグラフ描画
一次関数 <- lm(y ~ x)
abline(一次関数, lty=1, lwd=3, col="red")
# 回帰直線の係数
一次関数
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
-0.00838 1.01195
# 回帰直線とデータとのずれ（残差）
names(一次関数)
[1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "xlevels" "call" "terms" "model"
一次関数$model
y x
1 -0.35852 -0.33492
2 -0.07079 -0.02602
3 0.91555 0.95281
4 0.79529 0.76055
5 0.55257 0.51180
6 0.33522 0.28548
7 -0.94529 -0.82717
8 -0.96960 -0.82899
9 0.22960 0.20065
10 -0.54989 -0.50108
11 -0.84961 -0.73087
12 -0.21420 -0.01726
13 -0.44414 -0.49982
14 -0.64574 -0.53294
15 -1.11760 -0.96354
16 0.72881 0.77916
17 0.55585 0.59002
18 -0.14807 0.03840
19 0.75096 0.85227
20 0.76323 0.77895
21 -0.87517 -0.86019
22 -0.74306 -0.72071
23 0.69338 0.81584
24 0.78240 0.86429
25 0.31043 0.40376
26 -0.28608 -0.39326
27 -0.78560 -0.68612
28 -0.48722 -0.45481
29 0.64619 0.56000
30 0.56673 0.73452
31 -0.38878 -0.64531
32 0.57643 0.54830
33 0.23053 0.23982
34 0.30289 0.27879
35 -0.85048 -0.92710
36 0.73998 0.89237
37 0.68655 0.53020
38 0.41003 0.43020
39 0.85073 0.84653
40 0.05793 0.06814
41 0.73839 0.78231
42 -0.32811 -0.27692
43 0.51435 0.46392
44 -0.52885 -0.66913
45 -0.85200 -0.64663
46 -0.86279 -0.83766
47 -0.90867 -0.71797
48 -0.77692 -0.66792
49 0.72822 0.75770
50 0.49789 0.35528
51 0.75949 0.65480
52 0.60764 0.50931
53 -0.70968 -0.51721
54 -0.60180 -0.63026
55 0.14815 0.01920
56 -0.90084 -0.89746
57 -0.31352 -0.27847
58 1.02635 0.88324
59 -0.97194 -0.92574
60 -0.45261 -0.33635
61 -0.13461 -0.17024
62 0.12180 0.23598
63 0.30040 0.28418
64 -0.86991 -0.98516
65 0.12076 0.11456
66 0.79624 0.66769
67 -0.40928 -0.37816
68 -0.16352 0.03153
69 -0.57525 -0.60534
70 -0.34387 -0.27599
71 0.20740 0.16551
72 -0.19963 -0.20896
73 -0.02343 -0.10779
74 0.47008 0.46589
75 -0.30181 -0.34277
76 0.59724 0.55365
77 0.76725 0.70910
78 -0.65103 -0.71501
79 0.24730 0.26887
80 -0.91200 -0.93674
81 0.16208 0.19778
82 0.05062 -0.11220
83 -0.39939 -0.37575
84 -0.42033 -0.30397
85 0.26113 0.34288
86 -0.36946 -0.43873
87 -0.61013 -0.50062
88 0.51534 0.38245
89 -0.23969 -0.48950
90 -0.82850 -0.73703
91 -0.15051 -0.26029
92 0.36008 0.53309
93 -0.14781 -0.24096
94 -0.80938 -0.89220
95 -0.53752 -0.52870
96 0.45444 0.30776
97 -0.62224 -0.58966
98 -0.09383 -0.01060
99 0.71852 0.63696
100 0.50195 0.55740
fitted <- predict(一次関数)
for (i in 1:100)
lines (
c(x[i],x[i]),
c(y[i],fitted[i])
)
#【1】曲線「y = 1 - x^2」に正規分布誤差を付加する
n <- 100
x <- runif(n, min = -1, max = 1)
y <- (1 - x^2) + rnorm(n, 0, 0.04)
plot(x, y)
# 回帰直線とそのグラフ描画
一次関数 <- lm(y ~ x)
abline(一次関数, lty=1, lwd=3, col="red")
# 回帰二次関数とそのグラフ描画
x2 <- x^2
二次関数 <- lm(y ~ x + x2)
xx <- seq(-1, 1, 0.01)
yy <- predict(二次関数, list(x=xx, x2=xx^2))
lines(xx, yy, lty=1, lwd=3, col="blue")
# 回帰二次関数の係数
二次関数
Call:
lm(formula = y ~ x + x2)
Coefficients:
(Intercept) x x2
0.99839 0.00238 -0.99833
#【3】曲線「y = (x - 0.9)*(x - 0.6)*(x - 0.1)*(x + 0.3)*(x + 1)」に正規分布誤差を付加する
n <- 100
x <- runif(n, min = -1, max = 1)
y <- (x - 0.9)*(x - 0.6)*(x - 0.1)*(x + 0.3)*(x + 1) + rnorm(n, 0, 0.04)
plot(x, y)
# 回帰直線とそのグラフ描画
一次関数 <- lm(y ~ x)
abline(一次関数, lty=1, lwd=3, col="red")
# 回帰二次関数とそのグラフ描画
x2 <- x^2
二次関数 <- lm(y ~ x + x2)
xx <- seq(-1, 1, 0.01)
yy <- predict(二次関数, list(x=xx, x2=xx^2))
lines(xx, yy, lty=1, lwd=3, col="blue")
# 回帰三次関数とそのグラフ描画
x3 <- x^3
三次関数 <- lm(y ~ x + x2 + x3)
xx <- seq(-1, 1, 0.01)
yy <- predict(三次関数, list(x=xx, x2=xx^2, x3=xx^3))
lines(xx, yy, lty=1, lwd=3, col="green")
# 回帰四次関数とそのグラフ描画
x4 <- x^4
四次関数 <- lm(y ~ x + x2 + x3 + x4)
xx <- seq(-1, 1, 0.01)
yy <- predict(四次関数, list(x=xx, x2=xx^2, x3=xx^3, x4=xx^4))
lines(xx, yy, lty=1, lwd=3, col="pink")
# 回帰五次関数とそのグラフ描画
x5 <- x^5
五次関数 <- lm(y ~ x + x2 + x3 + x4 + x5)
xx <- seq(-1, 1, 0.01)
yy <- predict(五次関数, list(x=xx, x2=xx^2, x3=xx^3, x4=xx^4, x5=xx^5))
lines(xx, yy, lty=1, lwd=3, col="black")
# 回帰六次関数とそのグラフ描画
x6 <- x^6
六次関数 <- lm(y ~ x + x2 + x3 + x4 + x5 + x6)
xx <- seq(-1, 1, 0.01)
yy <- predict(六次関数, list(x=xx, x2=xx^2, x3=xx^3, x4=xx^4, x5=xx^5, x6=xx^6))
lines(xx, yy, lty=1, lwd=3, col="gray")
logLik(一次関数)
'log Lik.' 132.6 (df=3)
logLik(二次関数)
'log Lik.' 150.9 (df=4)
logLik(三次関数)
'log Lik.' 152.1 (df=5)
logLik(四次関数)
'log Lik.' 152.1 (df=6)
logLik(五次関数)
'log Lik.' 177.8 (df=7)
logLik(六次関数)
'log Lik.' 177.9 (df=8)
AIC(一次関数)
[1] -259.2
AIC(二次関数)
[1] -293.8
AIC(三次関数)
[1] -294.2
AIC(四次関数)
[1] -292.3
AIC(五次関数)
[1] -341.7
AIC(六次関数)
[1] -339.8