Rikon. # x y Rikon Zouka. Ninzu -.3 Kaku. Tomo -1.7 Tandoku.7 X5Sai -1. Kfufu -.5 Ktan.13 Konin.7 # x Zouka Ninzu

Transcription

1 Copyright (c),5 Hidetoshi Shimodaira -1-13:1: shimo R [ ] x11 x1p y1 = x1,..., xp, y =. xn1 xnp yn i = 1,..., n i (xi1, xi,..., xip, yi) (x1, y1),..., (xn, yn) p 1 (x11, x1,..., x1p, y1),..., (xn1, xn,..., xnp, yn) i = 1,..., n (xi1, xi,..., xip, yi) x11 x1p y1 x1 =,..., xp =, y = xn1 xnp j = 1,..., p i j xij x1j xj = xnj [ ] x11 x1p X = x1,..., xp =. 1 xn1 xnp yn # run5.r # dat <- read.table("dat.txt") # cat("\n# \n"); dim(dat) cat("\n#,,3 \n"); dat[1:3,] cat("\n# j=1,,3 \n"); dat[,1:3] x <- dat[,-1]; # 1,...,9 y <- dat[,1]; # 1 cat("\n# x,,3 \n"); x[1:3,] cat("\n# y,,3 \n"); y[1:3] A <- cbind(x,y) # x y cat("\n# A,,3 ( \n"); A[1:3,] y <- dat[,1,drop=f]; # 1 cat("\n# y,,3 (drop=f )\n"); y[1:3,,drop=f] A <- cbind(x,y) # x y cat("\n# A,,3 ( \n"); A[1:3,] > source("run5.r",print=t) # [1] 7 1 #,,3 Rikon Hokkaido Aomori Iwate # j=1,,3 Zouka Ninzu Kaku Hokkaido...5 Aomori Iwate Miyagi Kumamoto Ooita Miyazaki Kagoshima Okinawa # x,,3 Hokkaido Aomori Iwate # y,,3 [1] # A,,3 ( y Hokkaido Aomori Iwate # y,,3 (drop=f ) Rikon Hokkaido. Aomori 1.9 Iwate 1. # A,,3 ( Rikon Hokkaido Aomori Iwate ( ) : j = 1,..., p xj xj xj = 1 n xij n xj xk () j, k = 1,..., p 1 n Sxjxk = Sjk = n 1 (xij xj)(xik xk) 3 xj xk ) j, k = 1,..., p rxjxk = rjk = Sjk SjjSkk # run57.r # pairs(x) ; dev.copyeps(file="run57-s.eps") pr <- function(a) print(round(a,3)) cat("\n# x \n"); pr(mean(x)) cat("\n# y \n"); pr(mean(y)) cat("\n# x \n"); pr(var(x)) cat("\n# y \n"); pr(var(y)) cat("\n# x y \n"); pr(var(x,y)) cat("\n# x \n"); pr(cor(x)) cat("\n# y?\n"); pr(cor(y)) cat("\n# x y \n"); pr(cor(x,y)) > source("run57.r") # x # y Rikon 1. # x Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # y Rikon

2 Rikon. # x y Rikon Zouka. Ninzu -.3 Kaku. Tomo -1.7 Tandoku.7 X5Sai -1. Kfufu -.5 Ktan.13 Konin.7 # x Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # y? Rikon Rikon 1 # x y Rikon Zouka.31 Ninzu -. Kaku.5 Tomo -.39 Tandoku.57 X5Sai -.71 Kfufu -.53 Ktan.5 Konin Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin run57-s 1.3 xj = 1 1 [ ] n 1 nxj 1n =., 1 n = A A ( x1,..., xp) = 1 n 1 nx S = (Sij) = 1 n 1 Z Z Z = X 1 n 1n1 nx S = 1 n 1 X J nx J n = In 1 n 1n1 n In n n J nj n = (J n) = J n > ## > n <- nrow(x) # > n [1] 7 > rep(1,n) # 1 7 [1] [39] > pr((1/n)*rep(1,n) %*% x) # x Error in rep(1, n) %*% x : requires numeric matrix/vector arguments > x <- as.matrix(x) # > pr(x[1:3,]) # Hokkaido Aomori Iwate > pr((1/n)*rep(1,n) %*% x) # [1,] > pr(mean(x)) # [1] > pr(apply(x,,mean)) # > xm <- (1/n)*rep(1,n) %*% (rep(1,n) %*% x) # > pr(xm[1:3,]) [1,] [,] [3,] > xc <- x - xm # > pr(xc[1:3,]) # Hokkaido Aomori Iwate > pr(rep(1,n) %*% xc) # [1,] > v <- (1/(n-1)) * (t(xc) %*% xc); pr(v) # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin > pr(var(x)) # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin > round(v - var(x),1) # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin > J <- diag(n) - (1/n) * rep(1,n) %o% rep(1,n) # J_n > pr((1/(n-1))* t(x) %*% J %*% x) #

3 Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin > pr(diag(v)) # > pr(sqrt(diag(v) %o% diag(v))) # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin > pr(v / sqrt(diag(v) %o% diag(v)) ) # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin > pr(cor(x)) # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin > round(v / sqrt(diag(v) %o% diag(v)) - cor(x),1) # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin.1 (x1,..., xp) y y = β + β 1 x β p x p + ɛ x1,..., xp: p y: ɛ: β, β1,..., βp:. ) β, β1,..., βp i = 1,..., n yi = β + β1xi1 + + βpxip + ɛi β, β1,..., βp n ɛ i = (β + β1xi1 + + βpxip yi) 1 β (p + 1) 1 X 1n n (p + 1) β β1 [ ] 1 x11 x1p β =, X = 1n, x1,..., xp =.. 1 xn1 xnp βp y ɛ n 1 y1 ɛ1 y =, ɛ = yn ɛn y = Xβ + ɛ ɛ = Xβ y ˆβ ŷ ŷ = X ˆβ optim # run5.r # ( ) # x y n <- nrow(x); p <- ncol(x) xx <- as.matrix(cbind(1,x)) # 1 rss <- function(be) sum((xx %*% be - y)^) # be <- rep(,p+1) # be1 <- optim(be,rss)$par # names(be1) <- dimnames(xx)[[]] cat("# \n"); print(be1) # cat("# \n"); print(rss(be1)) # pred1 <- xx %*% be1 # \hat y plot(pred1,y,pch=1) # abline(,1,col=) # y 11 dev.copyeps(file="run5-s1.eps") be <- optim(be,rss,method="bfgs")$par # names(be) <- dimnames(xx)[[]] cat("# \n"); print(be) # cat("# \n"); print(rss(be)) # pred <- xx %*% be # \hat y be3 <- optim(be,rss,method="bfgs", control=list(reltol=1e-1))$par # names(be3) <- dimnames(xx)[[]] cat("# \n"); print(be3) # cat("# \n"); print(rss(be3)) # pred3 <- xx %*% be3 # \hat y plot(pred3,y,pch=1) # abline(,1,col=) # y dev.copyeps(file="run5-s3.eps") plot(pred1,pred3,pch=1); abline(,1,col=) dev.copyeps(file="run5-s13.eps") plot(pred,pred3,pch=1); abline(,1,col=) dev.copyeps(file="run5-s3.eps") > dat <- read.table("dat.txt") # (7 x 1 ) > x <- dat[,-1]; y <- dat[,1] > source("run5.r") # 1 Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # [1].913 # 1 Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # [1].73 # 1 Zouka Ninzu Kaku Tomo Tandoku 1

4 y y y X5Sai Kfufu Ktan Konin # [1] ˆβ = (X X) 1 X y A A ɛ = (y Xβ) (y Xβ) = y y β X y + β (X X)β β ɛ β = X y + X Xβ = (normal equation) X Xβ = X y pred1 run5-s pred3 run5-s3 ˆβ = (X X) 1 X y X X pred pred pred1 pred run5-s13 run5-s3 optim ˆβ = X + y A A + X + = (X X) 1 X ) solve() # run59.r # ( ) # x y n <- nrow(x); p <- ncol(x) xx <- as.matrix(cbind(1,x)) # 1 rss <- function(be) sum((xx %*% be - y)^) # A <- t(xx) %*% xx # A = X X be <- solve(a) %*% (t(xx) %*% y) # beta = A^-1 (X y) cat("# \n"); print(be) # cat("# \n"); print(rss(be)) # pred <- xx %*% be # \hat y plot(pred,y,pch=1) # abline(,1,col=) # y dev.copyeps(file="run59-s.eps") 13 1 > source("run59.r") # [,1] Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan.1555 Konin # [1].9979 > be3 - be [,1] e-3 Zouka -5.3e- Ninzu 1.35e-3 Kaku 5.197e-5 Tomo.7e- Tandoku.1577e-5 X5Sai 1.e-5 Kfufu e-5 Ktan -.733e-5 Konin 7.131e-5 > rss(be3) - rss(be) [1].9557e pred run59-s optim optim optim. R lsfit() lm() # run.r # ( ) # x y cat("# lsfit \n") fit1 <- lsfit(x,y) # QR print(fit1$coef) # cat("# lm \n") fit <- lm(y ~., data.frame(x,y)) # print(fit$coef) # cat("# lsfit \n") ls.print(fit1) # cat("# lm \n") print(summary(fit)) # 15 1

5 > source("run.r") # lsfit Intercept Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # lm (Intercept) Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # lsfit Residual Standard Error=.13 R-Square=.15 F-statistic (df=9, 37)=1.539 p-value= Estimate Std.Err t-value Pr(> t ) Intercept Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin Zouka * Ninzu ** Kaku Tomo Tandoku * X5Sai Kfufu ** Ktan Konin Signif. codes: ***.1 **.1 * Residual standard error:.13 on 37 degrees of freedom Multiple R-Squared:.15,Adjusted R-squared:.773 F-statistic: 1.53 on 9 and 37 DF, p-value: 3.51e-11 > names(fit1) # lsfit() [1] "coefficients" "residuals" "intercept" "qr" > names(fit) # lsfit() [1] "coefficients" "residuals" "effects" "rank" [5] "fitted.values" "assign" "qr" "df.residual" [9] "xlevels" "call" "terms" "model" > names(ls.print(fit1,print=f)) # lsfit() [1] "summary" "coef.table" > names(summary(fit)) # lm() [1] "call" "terms" "residuals" "coefficients" [5] "aliased" "sigma" "df" "r.squared" [9] "adj.r.squared" "fstatistic" "cov.unscaled" # lm lm(formula = y ~., data = data.frame(x, y)) (Intercept) * 17.5 (x1,..., xp) ŷ ŷ = ˆβ + ˆβ1x1 + + ˆβpxp i = 1,..., n ŷi ŷi = ˆβ + ˆβ1xi1 + + ˆβpxip ŷ = X ˆβ ˆβ = (X X) 1 X y 1 ŷ = X ˆβ ŷ = X(X X) 1 X y H = X(X X) 1 X ŷ = Hy H n n H = H H y ŷ i = 1,..., n ei ei = yi ŷi e1 e = = y ŷ en e = (In H)y ē = ŷ ȳ ŷ = 1 n ŷi = 1 n (yi ei) = ȳ ē = ȳ n n ē = HX = X(X X) 1 X X = X e = (In H)y X e X e = X (In H)y = (X X H)y = (X HX) y = X H = X H = (HX) X 1n X e = ē = 1 n 1 ne = n ei = 1 ne n = 1 n (ei )(ŷi n 1 ŷ) = 1 n eiŷi = e ŷ n 1 n 1 = X e = ŷ = X ˆβ e ŷ = e X ˆβ = (X e) ˆβ = ˆβ = 19 n e i = e (RSS, residuals sum of squares) ɛ σ Sɛ 1 n = e i n p 1 β, β1,..., βp p+1 n (p+1) # run1.r # ( ) # x y n <- nrow(x); p <- ncol(x) xx <- as.matrix(cbind(1,x)) # 1 rss <- function(be) sum((xx %*% be - y)^) # A <- t(xx) %*% xx # A = X X H <- xx %*% solve(a) %*% t(xx) # pred <- H %*% y # resid <- y - pred # plot(pred,resid,pch=1) # se <- sum(resid^)/(n - p - 1) # se <- sqrt(se) # epsilon abline(h=) # abline(h=c(se,-se),lty=,col=) # +-se abline(h=c(*se,-*se),lty=3,col=3) # +-se title(sub=paste("sd=",signif(se,5)),cex.sub=) dev.copyeps(file="run1-e.eps") source("run.r") # drawhist drawhist(resid,,"resid","run1-") # > source("run1.r") > mean(resid) # [1].1111e-1 > sum(resid) # [1] e-11 > cov(resid,pred) # [,1] [1,].195e-13 > sum(resid*pred) #

6 [1].5155e-1 > resid*pred # [,1] Hokkaido.7575 Aomori Iwate Miyagi Ooita Miyazaki.339 Kagoshima Okinawa.535 > dim(h) # n * n [1] 7 7 > H[1:3,1:3] # 3 * 3 Hokkaido Aomori Iwate Hokkaido Aomori Iwate > diag(h)[1:3] # Hokkaido Aomori Iwate > hat(x) # diag(h) [1] [7] [13] [19] [5] [31] [37] [3] resid pred sd=.135 run1-e. 1 3 y = ŷ + y ŷ resid mean=.11e 1, sd=.1 run1-resid y ŷ ŷ y ŷ ŷ e = y ŷ ŷ e = y = ŷ + e y = ŷ + e = ŷ + ŷ e + e = ŷ + e ȳ ȳ = 1nȳ = 1 n 1n1 ny y ŷ + ŷ ȳ = y ȳ ȳ y ŷ ȳ ŷ y ŷ ŷ e = y ŷ y ȳ = ŷ ȳ + e y ȳ = ŷ ȳ + e = ŷ ȳ + (ŷ ȳ) e + e = ŷ ȳ + e (ŷ ȳ) e = 1 shaei1x13 shaeiy311 e e.997 R R y ŷ n (yi ȳ)(ŷi ȳ) R = ry ŷ = n (yi ȳ) n (ŷi ȳ) 1 R 1 R 1 R x = 1 x = 1n X = [x,..., xp] x, x1,..., xp Span = ) Span(x, x1,..., xp) = { p j= βjxj β,..., βp R } R n y Span(X) = Span(x, x1,..., xp) Proj Span(X) (y) min y p p Span(X) p Projection= ) p Span(X) β p = Xβ y p = y Xβ Proj Span(X) (y) = ŷ H ŷ = Hy Proj Span(X) (y) = Hy H Span(X).7 ŷ = ˆβ + ˆβ1x ŷ ȳ ˆβ1(x x) n = (ŷi n ȳ) ˆβ = 1(xi x) ry ŷ = ry x x y R = (y ȳ) (ŷ ȳ) y ȳ ŷ ȳ = cos φ x x n (xi x) y ȳ ŷ ȳ φ a a/ a u v φ u v = cos φ R (y ȳ) (ŷ ȳ) = (y ȳ) H(y ȳ) = (y ȳ) H (y ȳ) = H(y ȳ) = ŷ ȳ H = H R ŷ ȳ R = y ȳ ŷ ȳ = ŷ ȳ y ȳ R = cos φ (ŷ ȳ) R = (y ȳ) y = Xβ + ɛ 3

7 R ŷ ȳ = y ȳ R = y ȳ e y ȳ = 1 e y ȳ R = cos φ = 1 sin φ R e e y ȳ # run.r # # x y fit <- lm(y ~., data.frame(x,y)) # resid <- fit$resid # pred <- y - resid # bary <- mean(y) # y cat("\n# \n"); print(fit$coef) cat("# R = ",cor(pred,y),"\n") cat("# R^ = ",cor(pred,y)^,"\n") cat("# R^ = ",sum((pred-bary)^)/sum((y-bary)^),"\n") cat("# R^ = ",1-sum(resid^)/sum((y-bary)^),"\n") cat("# R^ = ",summary(fit)$r.squared,"\n") > source("run.r") # (Intercept) Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # R =.97 # R^ =.117 # R^ =.117 # R^ =.117 # R^ =.117. ɛ1,..., ɛn σ ɛ1,..., ɛn N(, σ ) l(β, β1,..., βp, σ) = n n ) (yi β β1xi1 βpxip) log(πσ ) ( σ β,..., βp, σ l β =,..., l =, βp l = β1 ˆβ = (X X) 1 y, ˆσ = e n l( ˆβ, ˆβ1,..., ˆβp, ˆσ) = n { 1 + log(πˆσ ) } 1. β. σ ˆσ = n p 1 S n ɛ ˆβ,..., ˆβp Sɛ R (xi1, xi,..., xip, yi), i = 1..., n (xi1, xi,..., xip, yi) p + 1 (X1,..., Xp, Y ) 5 # run3.r # # x y source("run.r") # drawhist fit <- lsfit(x,y) # QR func3 <- function(f,y) { # fit coef <- f$coef # resid <- f$resid # se <- sum(resid^)/(length(resid) - length(coef)) # se <- sqrt(se) # (Se) rsq <- cor(y-resid,y)^ # ( y!) names(se) <- "sigma"; names(rsq) <- "rsquared" c(coef,se,rsq) } cat("\n# \n"); print(func3(fit,y)) boot1 <- function(i) # i = n func3(lsfit(x[i,],y[i]),y[i]) # n <- length(y) b <- 1 # cat("\n# \n"); print(boot1(1:n)) simi <- matrix(,n,b) # print(date()) for(j in 1:b) simi[,j] <- sample(1:n,replace=t) # print(date()) simt <- apply(simi,,boot1) # print(date()) cat("\n# \n") a <- apply(simt,1,function(x) unlist(list(mean=mean(x),sd=sd(x)))) print(a) for(k in rownames(simt)) drawhist(simt[k,],,k,"run3-") # cat("\n# lm() \n") print(summary(lm(y~.,data.frame(x,y)))$coef) > dat <- read.table("dat.txt") # (7 x 1 ) > x <- dat[,-1]; y <- dat[,1] > source("run3.r") # Intercept Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin sigma rsquared # Intercept Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin sigma rsquared [1] "Fri Oct 1 1:3:1 " [1] "Fri Oct 1 1:3:19 " [1] "Fri Oct 1 1:3: " # Intercept Zouka Ninzu Kaku Tomo Tandoku mean sd X5Sai Kfufu Ktan Konin sigma rsquared mean sd # lm() (Intercept) Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin Intercept mean= 15.13, sd=.33 run3-intercept Zouka mean= 1.57, sd=.75 run3-zouka Ninzu mean=.95, sd= run3-ninzu 1 Kaku mean=.3119, sd=.397 run3-kaku

8 Tandoku X5Sai # あらかじめ x に説明変数の行列 y に目的変数のベクトルを設定しておく Kfufu 5 1 Tomo 3 pred <- y - fit$resid # 予測値 resid <- fit$resid/sqrt(1-hat(x)) # 修正残差 mean=.15991, sd= mean=.179, sd=.519 run3-tomo run3-kfufu sigma Konin.1 mean=.171, sd=.73 run3-x5sai rsquared dev.copyeps(file="run-hat.eps") yi <- pred + resid[i] # ブートストラップ標本 func3(lsfit(x,yi),yi) # 重回帰分析 } # i = サイズ n の添え字ベクトル boot1 <- function(i) { 1 Ktan.1 mean=.93, sd=.11 run3-tandoku plot(pred,1/sqrt(1-hat(x))) 1 cat("\n# 統計量\n"); print(func3(fit,y)) 3 fit <- lsfit(x,y) # 最小二乗法の計算 QR 分解 mean=.1193, sd= mean=.115, sd=.15 mean=.533, sd=.131 run3-ktan mean=.51, sd=.19 run3-sigma run3-konin b <- 1 # シミュレーションの繰り返し回数 n <- length(y) run3-rsquared cat("\n# 統計量をオリジナルデータへ適用\n"); print(boot1(1:n)) simi <- matrix(,n,b) # ブートストラップ標本の添え字アレイを準備 print(date()) ブートストラップの繰り返し回数=1 だが ３０秒程度で終了した for(j in 1:b) simi[,j] <- sample(1:n,replace=t) # ブートストラップ法 回帰係数 β,..., β 1 のバラツキ (sd) は lm() の返す Std.Error と大体同じ値 print(date()) ブートストラップ法では 他の統計量 (S² と R ) のバラツキも簡単に計算できる S² =.135 ±.1 R =.1 ±. だった simt <- apply(simi,,boot1) # 統計量を繰り返し適用 print(date()) cat("\n# 統計量の平均 標準偏差\n") a <- apply(simt,1,function(x) unlist(list(mean=mean(x),sd=sd(x)))) print(a) ブートストラップ法 その２ 3. for(k in rownames(simt)) drawhist(simt[k,],,k,"run-") # ヒストグラム 推定した回帰係数 β,..., β p および誤差標準偏差 S² のバラツキを調べる 決定係数 R のバラツキも調べる 残差のリサンプリング ハット行列で補正したバージョンを適用 データ (xi1, xi,..., xip, yi ), cat("\n# lm() を利用してチェック\n") print(summary(lm(y~.,data.frame(x,y)))$coef) > source("run.r") i = 1...,n # 統計量 をリサンプリング 各時点の (xi1, xi,..., xip, yi ) のうち 説明変数 (xi1, xi,..., xip ) は定 Intercept Zouka 数と考える yi だけが 確率変数 修正残差 ri = ei / 1 hii, i = 1,..., n をリサンプリングして ブートストラップ標本 r1,..., rn を作成 これから y1,..., yn のブートストラップ標本 y1,..., yn を次式でつくる yi = y i + ri, Ninzu Kaku X5Sai Kfufu Ktan Intercept Konin sigma rsquared Tomo Tandoku Zouka 1.97 Ninzu X5Sai # 重回帰分析 ブートストラップ法 その２ Kaku Kfufu Ktan # run.r Konin sigma rsquared [1] "Fri Oct 1 1:5:5 " [1] "Fri Oct 1 1:5:5 " ::19 " Intercept Ninzu Kaku Zouka Tandoku.9739 Konin sigma mean= 15.9, sd=.3 run-intercept mean= , sd=.5593 mean= 3.71, sd= 1.15 run-zouka run-ninzu Tomo Tandoku mean=.7, sd=.373 run-kaku X5Sai Kfufu 3 # lm() を利用してチェック sd rsquared mean Ktan Kfufu X5Sai sd. mean Tomo. Kaku.3 Ninzu.3 Zouka # 統計量の平均 標準偏差 Intercept Tandoku # 統計量をオリジナルデータへ適用 i = 1,..., n y i = β + β 1 xi1 + + β p xip [1] "Fri Oct Tomo Pr(> t ) mean=.35, sd=.339 run-x5sai Ktan..3 Konin..1. mean=.17, sd=.195 run-kfufu sigma rsquared Tandoku mean=.1579, sd=.59 run-ktan mean=.1375, sd=.1775 mean=.7, sd=.1173 run-sigma run-konin mean=.1, sd=.5 run-rsquared 1.7 ブートストラップの繰り返し回数=1 だが ３０秒程度で終了した 1. 回帰係数 β,..., β 1 のバラツキ (sd) は lm() の返す Std.Error とほとんど同じ値 1.5 ブートストラップ法では 他の統計量 (S² と R ) のバラツキも簡単に計算できる S² = ±.1 R =.1 ±.5 だった 正規回帰モデル 回帰係数のバラツキを計算する理論を学ぶ (lm や lsfit の結果からＲが出力する 1.1 1/sqrt(1 hat(x)) Konin. Ktan Kfufu X5Sai Tomo. mean=.1195, sd=.539 run-tandoku Kaku run-tomo 1. mean=.151, sd=.97 Ninzu Zouka t value Std. Error 1 Estimate (Intercept) Std.Error のこと) ブートストラップ法を知っていれば 理論は知らなくても結果は出せる でもやはり理論. pred の勉強はしておくべき run-hat まず 説明変数は定数であるとみなす 残差のリサンプリングと同じ仮定 正規回帰モデルでは 誤差 ²1,..., ²n が互いに独立に そして説明変数に依存せずに 平 均０ 分散 σ の正規分布に従うと仮定する ²1,..., ²n N (, σ ) 31 3

9 ɛ σ In n y ɛ ɛ Nn(, σ In) E(ɛ) =, V (ɛ) = σ In y = Xβ + ɛ ɛ y E(y) = E(Xβ) + E(ɛ) = Xβ, V (y) = σ In y Nn(Xβ, σ In) ˆβ = (X X) 1 X y y ˆβ E( ˆβ) = (X X) 1 X E(y) = (X X) 1 X Xβ = β V ( ˆβ) = (X X) 1 X V (y)x(x X) 1 = σ (X X) 1 X X(X X) 1 = σ (X X) 1 (p + 1) (p + 1) A ˆβ Np+1(β, σ (X X) 1 ) A = (aij) = (X X) 1, ˆβi N(βi, σ aii) i, j =,..., p i j ˆβi ˆβj cov( ˆβi, ˆβj) = σ aij, ρ( ˆβi, ˆβj) = aij aiiajj 33 V ( ˆβi) = σ aii ˆσ S ɛ ˆV ( ˆβi) = S ɛ aii se( ˆβi) = t-p-value ti = ˆβi, Sɛ aii T n p 1 t- # run5.r # ˆV ( ˆβi) = Sɛ aii pi = Pr{T > ti } # x y n <- nrow(x); p <- ncol(x) X <- as.matrix(cbind(1,x)) # 1 A <- solve(t(x) %*% X) # A = (X X)^-1 coef <- A %*% (t(x) %*% y) # pred <- X %*% coef # resid <- y - pred # se <- sum(resid^)/(n-p-1) # se <- sqrt(se) # coefsd <- se*sqrt(diag(a)) # tval <- coef/coefsd # t- pval <- *pt(abs(tval),df=n-p-1,lower.tail=f) cat("# t-p- \n"); a <- cbind(round(coef,),round(coefsd,),round(tval,3),round(pval,5)) colnames(a) <- c("coef","sd","t-val","p-val"); print(a) coefcor <- A / sqrt(diag(a) %o% diag(a)) # cat("# \n"); print(round(coefcor,)) fit <- lm(y~.,data.frame(x,y)) # cat("# lm() \n") print(summary(fit,cor=t)) > source("run5.r") # t-p- coef sd t-val p-val Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin # lm() lm(formula = y ~., data = data.frame(x, y)) (Intercept) * Zouka * Ninzu ** Kaku Tomo Tandoku * X5Sai Kfufu ** Ktan Konin Signif. codes: ***.1 **.1 * Residual standard error:.13 on 37 degrees of freedom Multiple R-Squared:.15,Adjusted R-squared:.773 F-statistic: 1.53 on 9 and 37 DF, p-value: 3.51e-11 Correlation of (Intercept) Zouka Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Zouka. Ninzu Kaku Tomo Tandoku X5Sai Kfufu Ktan Konin e = y ŷ = (In H)y HX = X e = (In H)(Xβ + ɛ) = (In H)ɛ E(e) = (In H)E(ɛ) = V (e) = (In H)V (ɛ)(in H) = σ (In H) = σ (In H) e Nn(, σ (In H)) ei N(, σ (1 hii)), i = 1..., n H = (hij) i j ei ej ɛi ɛj cov(ei, ej) = σ hij hij ρ(ei, ej) = (1 hii)(1 hjj) 3

10 ei ri =, 1 hii ri N(, σ ) i = 1,..., n ɛi N(, σ ) ei ri Sɛ (standerdized residuals) MASS stdres internally Studentized residuals) si = ri (Studentized residuals) externally Studentized residuals, t-residuals) MASS studres ti Sɛ ei ti = s(i) 1 hii s (i) = (n p 1)S ɛ e i /(1 hii) n p s (i) i Sɛ ti n p t-.1 # run.r # if(!exists("x.data")) { } X.data <- read.table("xdata.txt") # X X.item <- read.table("xitem.txt") # mygetdat <- function(a,namae=null,echo=t) { # a= dat <- X.data[,a,drop=F] # imitan <- X.item[a,c("Imi","Tani"),drop=F] # if(!is.null(namae)) { # names(dat) <- namae # namae code <- rownames(imitan) # imitan <- cbind(code,imitan) 37 rownames(imitan) <- namae # =namae } if(echo) { cat("# =",nrow(dat),"*",ncol(dat),"\n") cat("# \n") print(dat[1:3,]) cat("# \n") print(imitan) } invisible(dat) } # run7.r # source("run.r") dat <- mygetdat(c("e95","a53"),c("gakureki","shushou")) fit <- lm(shushou~gakureki,dat) print(summary(fit)) plot(shushou~gakureki,dat,pch=1) # abline(fit,col=,lty=) # dev.copyeps(file="run7-s1.eps") plot(shushou~predict(fit),dat,pch=1) # = abline(,1,col=,lty=) # y= dev.copyeps(file="run7-s.eps") plot(predict(fit),resid(fit)/sqrt(1-hatvalues(fit))) # se <- sqrt(sum(resid(fit)^)/(nrow(dat)-)) # Se abline(h=,lty=3,col=); abline(h=*c(se,-se),lty=3,col=3) dev.copyeps(file="run7-z1.eps") plot(predict(fit),sqrt(abs(resid(fit)/(se*sqrt(1-hatvalues(fit)))))) abline(h=,lty=3,col=); abline(h=sqrt(),lty=3,col=3) dev.copyeps(file="run7-z.eps") plot(fit,which=1) # dev.copyeps(file="run7-z3.eps") plot(fit,which=3) # S-L dev.copyeps(file="run7-z.eps") predict(fit) fitted(fit) ŷ1,..., ŷn resid(fit) e1,..., en hatvalues(fit) > source("run7.r") 3 # = 7 * # Gakureki Shushou Hokkaido Aomori Iwate # code Imi Tani Gakureki E95 ( ) Shushou A53 lm(formula = Shushou ~ Gakureki, data = dat) resid(fit)/sqrt(1 hatvalues(fit)) predict(fit) sqrt(abs(resid(fit)/(se * sqrt(1 hatvalues(fit))))) predict(fit) (Intercept) < e-1 *** Gakureki e-9 *** Signif. codes: ***.1 **.1 * Residual standard error:.95 on 5 degrees of freedom Multiple R-Squared:.53,Adjusted R-squared:.5 F-statistic: 51. on 1 and 5 DF, p-value: 5.93e-9 Residuals.... run7-z1 Residuals vs Fitted Okinawa Saga Hokkaido Fitted values lm(formula = Shushou ~ Gakureki, data = dat) run7-z3 Standardized residuals run7-z Scale Location plot Okinawa Hokkaido Saga Fitted values lm(formula = Shushou ~ Gakureki, data = dat) run7-z Shushou Gakureki run7-s1 Shushou predict(fit) run7-s. # run.r # source("run.r") dat <- mygetdat(c("e95","a13","c131","b11","a53"), c("gakureki","mikonw","shotoku","kion","shushou")) fit <- lm(shushou~.,dat) print(summary(fit)) plot(shushou~predict(fit),dat,pch=1) # = abline(,1,col=,lty=) # y= dev.copyeps(file="run-s1.eps") 39

11 plot(fit,which=3) # S-L dev.copyeps(file="run-z1.eps") > source("run.r") # = 7 * 5 # Gakureki MikonW Shotoku Kion Shushou Hokkaido Aomori Iwate # code Imi Tani Gakureki E95 ( ) MikonW A13 ( ) Shotoku C131 ( :thousandyen) Kion B11 ( C) Shushou A53 lm(formula = Shushou ~., data = dat) (Intercept) 3.3e+.1e e-7 *** Gakureki -1.5e-.57e * MikonW -.7e- 7.39e *** Shotoku 1.9e-5.3e Kion.13e- 5.9e *** Signif. codes: ***.1 **.1 * Residual standard error:.7 on degrees of freedom Multiple R-Squared:.731,Adjusted R-squared:.717 F-statistic: 3.37 on and DF, p-value:.7e-1 Shushou Standardized residuals predict(fit) Fitted values lm(formula = Shushou ~., data = dat) run-s1.3 Scale Location plot Hokkaido Nagano Fukui run-z1 y = β + β1x + βx + + βpx p + ɛ # run9.r # source("run.r") dat <- mygetdat(c("b3","b1"),c("yuki","hkion")) cat("# 1 \n") x <- dat[,"yuki"] <- dat[,"yuki"]/1 fit1 <- lm(hkion~yuki,dat) # fit <- lm(hkion~yuki+i(yuki^),dat) # fit3 <- lm(hkion~yuki+i(yuki^)+i(yuki^3),dat) # cat("# \n"); print(summary(fit1)) cat("# \n"); print(summary(fit)) cat("# \n"); print(summary(fit3)) plot(hkion~yuki,dat,pch=1) # abline(fit1,col=) # x <- seq(min(x),max(x),length=3) # x 3 newdat <- data.frame(yuki=x) # lines(x,predict(fit,newdat),col=3) # lines(x,predict(fit3,newdat),col=) # dev.copyeps(file="run9-s.eps") > source("run9.r") 1 # = 7 * # Yuki Hkion Hokkaido 15.3 Aomori Iwate # code Imi Tani Yuki B3 ( ) Hkion B1 ( C) # 1 # lm(formula = Hkion ~ Yuki, data = dat) (Intercept) < e-1 *** Yuki ** Signif. codes: ***.1 **.1 * Residual standard error: 1.3 on 5 degrees of freedom Multiple R-Squared:.199,Adjusted R-squared:.111 F-statistic: on 1 and 5 DF, p-value:.179 # lm(formula = Hkion ~ Yuki + I(Yuki^), data = dat) (Intercept) < e-1 *** Yuki ** I(Yuki^) e-5 *** Signif. codes: ***.1 **.1 * Residual standard error: 1.71 on degrees of freedom Multiple R-Squared:.719,Adjusted R-squared:.79 F-statistic: 19. on and DF, p-value: 7.931e-7 # lm(formula = Hkion ~ Yuki + I(Yuki^) + I(Yuki^3), data = dat) (Intercept) < e-1 *** Yuki *** I(Yuki^) *** I(Yuki^3) * Signif. codes: ***.1 **.1 * Residual standard error: 1.7 on 3 degrees of freedom Multiple R-Squared:.53,Adjusted R-squared:.51 F-statistic: 17.9 on 3 and 3 DF, p-value: 1.7e-7 > model.matrix(fit3) (Intercept) Yuki I(Yuki^) I(Yuki^3) Hokkaido Aomori Iwate Miyagi Ooita Miyazaki Kagoshima

12 Okinawa 1... attr(,"assign") [1] 1 3 # run7.r # dat$yuki <- x^; dat$yuki3 <- x^3 fit3b <- lm(hkion~yuki+yuki+yuki3,dat) cat("# \n"); print(summary(fit3b)) plot(hkion~yuki,dat,pch=1) # newdat$yuki <- x^; newdat$yuki3 <- x^3 pred3b <- as.matrix(cbind(1,newdat)) %*% coef(fit3b) Hokkaido Aomori Iwate > dim(newdat) [1] 3 3 > newdat[1:3,] Yuki Yuki Yuki3 1..e+.e e e e e-7 lines(x,pred3b,col=) # dev.copyeps(file="run7-s.eps") > source("run7.r") # lm(formula = Hkion ~ Yuki + Yuki + Yuki3, data = dat) Hkion Hkion Yuki run9-s Yuki run7-s (Intercept) < e-1 *** Yuki *** Yuki *** Yuki * Signif. codes: ***.1 **.1 * y = β + β1x1 + βx + ɛ 1. y. x1 3. x = 1 x = Residual standard error: 1.7 on 3 degrees of freedom Multiple R-Squared:.53,Adjusted R-squared:.51 F-statistic: 17.9 on 3 and 3 DF, p-value: 1.7e-7 > dim(dat) [1] 7 > dat[1:3,] Yuki Hkion Yuki Yuki3 5 # run71.r # y = (β + β) + β1x1 + ɛ y = β + β1x1 + ɛ source("run.r") source("japan.pref.r") cat("# japan.pref \n") print(japan.pref[1:5,]) dat <- mygetdat(c("c131","c"),c("shotoku","chokin")) dat <- cbind(dat,japan.pref) cat("# Shotoku 1 \n") dat[,"shotoku"] <- dat[,"shotoku"]/1 cat("# \n") fit1 <- lm(chokin~shotoku,dat) print(summary(fit1)) plot(chokin~shotoku,dat,type="n") text(dat$shotoku,dat$chokin,rownames(dat)) abline(fit1,col=) dev.copyeps(file="run71-s1.eps") cat("# \n") fit <- lm(chokin ~ Shotoku + I(Chihou=="Shikoku"), dat) print(summary(fit)) plot(predict(fit),dat$chokin,type="n") text(predict(fit),dat$chokin,rownames(dat)) abline(,1,col=) dev.copyeps(file="run71-s.eps") > source("run71.r") # japan.pref Jpref Jcity Longitude Latitude Umi Chihou Hokkaido Sonota Tohoku Aomori 1..9 Sonota Tohoku Iwate Taihei Tohoku Miyagi Taihei Tohoku Akita Nihon Tohoku # = 7 * # Shotoku Chokin Hokkaido Aomori Iwate 19.7 # code Imi Tani Shotoku C131 ( :thousandyen) Chokin C ( :1thousandyen) # Shotoku 1 7 # lm(formula = Chokin ~ Shotoku, data = dat) (Intercept) Shotoku e- *** Signif. codes: ***.1 **.1 * Residual standard error: 9. on 5 degrees of freedom Multiple R-Squared:.37,Adjusted R-squared:.35 F-statistic: 5.39 on 1 and 5 DF, p-value:.93e- # lm(formula = Chokin ~ Shotoku + I(Chihou == "Shikoku"), data = dat) (Intercept) Shotoku e- *** I(Chihou == "Shikoku")TRUE *** Signif. codes: ***.1 **.1 * Residual standard error: 79.7 on degrees of freedom Multiple R-Squared:.539,Adjusted R-squared:.5 F-statistic:. on and DF, p-value:.13e- > model.matrix(fit)

13 (Intercept) Shotoku I(Chihou == "Shikoku")TRUE Hokkaido Aomori Yamaguchi Tokushima Kagawa Ehime Kochi Fukuoka Saga Okinawa attr(,"assign") [1] 1 attr(,"contrasts") attr(,"contrasts")$"i(chihou == "Shikoku")" [1] "contr.treatment" Chokin Kagawa Wakayama Tokushima Fukui Gifu Toyama Nara Kyoto Osaka Nagano Hyogo Yamanashi Mie OkayamaIshikawa Hiroshima Shizuoka Kochi Ehime Gumma Yamaguchi Shiga Tochigi Tottori Niigata Ibaraki Kanagawa Shimane Saga Ooita Yamagata Fukuoka Hokkaido Nagasaki Iwate Fukushima Kumamoto Kagoshima Akita Miyagi Aomori Miyazaki Okinawa Chiba Saitama Shotoku run71-s1 Aichi Tokyo.5 1. y. x1 dat$chokin Kagawa Wakayama Tokushima Fukui Gifu Toyama Nara Kyoto Osaka Nagano Hyogo Aichi Yamanashi Mie OkayamaIshikawa Hiroshima Shizuoka Gumma Kochi Ehime Yamaguchi Shiga Tochigi Tottori Niigata Ibaraki Kanagawa Shimane Saga Ooita Yamagata Fukuoka Hokkaido Nagasaki Iwate Fukushima Kumamoto Kagoshima Akita Miyagi Aomori Miyazaki Okinawa y = β1x1 + βx + + β9x9 + ɛ 9 Chiba Saitama 5 7 predict(fit) run71-s Tokyo 3. x = 1 x =. x3 = 1 x3 = 5. x = 1 x =. x5 = 1 x5 = 7. x = 1 x =. x7 = 1 x7 = 9. x = 1 x = 1. x9 = 1 x9 = y = β + β1x1 + ɛ y = β3 + β1x1 + ɛ β # run7.r # cat("# beta Chihou \n") fit3 <- lm(chokin ~ -1 + Shotoku + Chihou, dat) print(summary(fit3)) plot(predict(fit3),dat$chokin,type="n") text(predict(fit3),dat$chokin,rownames(dat)) abline(,1,col=) dev.copyeps(file="run7-s1.eps") # beta Chihou lm(formula = Chokin ~ -1 + Shotoku + Chihou, data = dat) Shotoku *** ChihouChugoku ChihouKanto ChihouKinki * 5 ChihouKyushu ChihouShikoku ** ChihouShinetsu * ChihouTohoku ChihouTokai Signif. codes: ***.1 **.1 * Residual standard error: 57.7 on 3 degrees of freedom Multiple R-Squared:.991,Adjusted R-squared:.99 F-statistic: 7.1 on 9 and 3 DF, p-value: <.e-1 > 1-sum(resid(fit3)^)/sum(dat$Chokin^) # summary(lm()) R^ [1] > fit <- lm(chokin ~ 1,dat) # > 1-sum(resid(fit3)^)/sum(resid(fit)^) # R^ [1] > fit <- lm(chokin ~ Shotoku + Chihou, dat) # > summary(fit) # R^ lm(formula = Chokin ~ Shotoku + Chihou, data = dat) (Intercept) Shotoku *** ChihouKanto ChihouKinki ChihouKyushu ** ChihouShikoku ** ChihouShinetsu ChihouTohoku *** ChihouTokai Signif. codes: ***.1 **.1 * Residual standard error: 57.7 on 3 degrees of freedom 51 Multiple R-Squared:.7,Adjusted R-squared:.7371 F-statistic: 17.1 on and 3 DF, dat$chokin Kagawa Wakayama Tokushima Fukui Gifu Toyama Nara Kyoto Osaka Nagano Hyogo Aichi Yamanashi Mie OkayamaIshikawa Hiroshima Shizuoka Gumma Kochi Ehime Yamaguchi Shiga Tochigi Tottori Niigata Ibaraki Kanagawa Shimane Saga Ooita Yamagata Fukuoka Hokkaido Nagasaki Iwate Fukushima Kumamoto Kagoshima Akita Miyagi Aomori Miyazaki Okinawa Chiba Saitama 5 7 predict(fit) run71-s Tokyo p-value: 1.937e-1 dat$chokin Saga Ooita Yamagata Fukuoka Hokkaido Nagasaki Iwate Fukushima Kumamoto Kagoshima AkitaMiyagi Aomori Miyazaki Okinawa Kagawa Wakayama Tokushima Gifu Fukui Toyama Nagano NaraKyoto Osaka Hyogo Aichi Yamanashi Mie Okayama Hiroshima Ishikawa Shizuoka Gumma Kochi Ehime Yamaguchi Shiga Tochigi Tottori Ibaraki KanagawaNiigata Shimane predict(fit3) Saitama Chiba run7-s1 β R ( y = β + ɛ y = ɛ ) > model.matrix(fit3) Shotoku ChihouChugoku ChihouKanto ChihouKinki ChihouKyushu Hokkaido 73.1 Aomori.9 Iwate 1.9 Miyagi 77. Akita 57. Yamagata.9 Fukushima 73.7 Ibaraki Tochigi Gumma 3. 1 Saitama 3. 1 Chiba Tokyo 3. 1 Kanagawa Niigata 9.1 Toyama 9. Ishikawa 99. Fukui 9. Yamanashi.5 1 Nagano Tokyo

14 Gifu 93.1 Shizuoka 37.3 Aichi 359. Mie 7. Shiga Kyoto Osaka Hyogo 3. 1 Nara.7 1 Wakayama 3. 1 Tottori. 1 Shimane.5 1 Okayama 7. 1 Hiroshima Yamaguchi Tokushima 71. Kagawa.1 Ehime 5. Kochi 35.7 Fukuoka Saga Nagasaki. 1 Kumamoto Ooita. 1 Miyazaki Kagoshima Okinawa ChihouShikoku ChihouShinetsu ChihouTohoku ChihouTokai Hokkaido 1 Aomori 1 Iwate 1 Miyagi 1 Akita 1 Yamagata 1 Fukushima 1 Ibaraki Tochigi Gumma Saitama Chiba Tokyo Kanagawa 53 Niigata 1 Toyama 1 Ishikawa 1 Fukui 1 Yamanashi Nagano 1 Gifu 1 Shizuoka 1 Aichi 1 Mie 1 Shiga Kyoto Osaka Hyogo Nara Wakayama Tottori Shimane Okayama Hiroshima Yamaguchi Tokushima 1 Kagawa 1 Ehime 1 Kochi 1 Fukuoka Saga Nagasaki Kumamoto Ooita Miyazaki Kagoshima Okinawa attr(,"assign") [1] 1 attr(,"contrasts") attr(,"contrasts")$chihou [1] "contr.treatment" > model.matrix(fit) (Intercept) Shotoku ChihouKanto ChihouKinki ChihouKyushu Hokkaido Aomori 1.9 Iwate Miyagi Akita Yamagata 1.9 Fukushima Ibaraki Tochigi Gumma Saitama Chiba Tokyo Kanagawa Niigata Toyama 1 9. Ishikawa Fukui 1 9. Yamanashi Nagano Gifu Shizuoka Aichi Mie 1 7. Shiga Kyoto Osaka Hyogo Nara Wakayama Tottori 1. Shimane 1.5 Okayama 1 7. Hiroshima Yamaguchi Tokushima Kagawa 1.1 Ehime 1 5. Kochi Fukuoka Saga Nagasaki 1. 1 Kumamoto Ooita 1. 1 Miyazaki Kagoshima Okinawa ChihouShikoku ChihouShinetsu ChihouTohoku ChihouTokai Hokkaido 1 Aomori 1 Iwate 1 Miyagi 1 Akita 1 Yamagata 1 Fukushima 1 Ibaraki Tochigi Gumma Saitama Chiba Tokyo Kanagawa Niigata 1 Toyama 1 Ishikawa 1 Fukui 1 Yamanashi Nagano 1 Gifu 1 Shizuoka 1 Aichi 1 Mie 1 Shiga Kyoto Osaka Hyogo Nara Wakayama Tottori Shimane Okayama Hiroshima Yamaguchi Tokushima 1 Kagawa 1 5

15 Ehime 1 Kochi 1 Fukuoka Saga Nagasaki Kumamoto Ooita Miyazaki Kagoshima Okinawa attr(,"assign") [1] 1 attr(,"contrasts") attr(,"contrasts")$chihou [1] "contr.treatment" X japan.pref R ŷ y ŷ = y ŷ ri±sɛ 57