.3 ˆβ1 = S, S ˆβ0 = ȳ ˆβ1 S = (β0 + β1i i) β0 β1 S = (i β0 β1i) = 0 β0 S = (i β0 β1i)i = 0 β1 β0, β1 ȳ β0 β1 = 0, (i ȳ β1(i ))i = 0 {(i ȳ)(i ) β1(i ))

Size: px

Start display at page:

Download ".3 ˆβ1 = S, S ˆβ0 = ȳ ˆβ1 S = (β0 + β1i i) β0 β1 S = (i β0 β1i) = 0 β0 S = (i β0 β1i)i = 0 β1 β0, β1 ȳ β0 β1 = 0, (i ȳ β1(i ))i = 0 {(i ȳ)(i ) β1(i ))"

ちえこくまじ
6 years ago
Views:

1 Copright (c) 004,005 Hidetoshi Shimodaira :15:07 shimo cat(" 1: "); c(mea(), mea()) cat(" : "); mmea <- fuctio(v) sum(v)/legth(v) c(mmea(), mmea()) cat(" 3:\"); mea(dat) cat(" 1: "); c(var(), var()) cat(" : "); mvar <- fuctio(v) sum( (v-mmea(v))^ )/(legth(v)-1) c(mvar(), mvar()) cat(" 3:\"); var(dat) 1 (1, 1), (, ),..., (, ) ) ( ) =, ȳ = S = (1 ) + + ( ) = 1 (i ) cat(" 1: "); cov(,) # cat(" : ") mcov <- fuctio(u,v) sum( (u-mmea(u))*(v-mmea(v)) )/(legth(v)-1) mcov(,) cat(" 1: "); cov(,)/sqrt(var()*var()) cat(" : "); cor(,) cat(" 3:\"); cor(dat) ( ) ( ) S = (1 ȳ) + + ( ȳ) S = S, S = S = 1 (i ȳ) (1 )(1 ȳ) + + ( )( ȳ) S = = 1 (i )(i ȳ) ) # ru0031.r r = S SS = S S S dat <- read.table("dat0001.tt") (47 ) <- dat[,1]; <- dat[,] plot(,,pch=16) # plot(dat,pch=16) dev.copeps(file="ru0031-s.eps") ru0031-s > source("ru0031.r",prit=t) X11 1: [1] : [1] : Gakureki Shushou : [1] : [1] : Gakureki Shushou Gakureki Shushou : [1] : [1] : [1] : [1] : Gakureki Shushou Gakureki Shushou :,, :,, = β 0 + β 1 + ɛ optim # ru003.r ( ) dat <- read.table("dat0001.tt") (47 ) <- dat[,1]; <- dat[,] rss <- fuctio(be) sum( (be[1]+be[]* - )^ ) be <- optim(c(0,0),rss)$par # cat(": "); prit(be) # cat(" : "); prit(rss(be)) # plot(,,pch=16) # plot(dat,pch=16) ablie(be[1],be[],col=) dev.copeps(file="ru003-s.eps") th <- list(beta0=seq(1,.5,legth=100), beta1=seq(-0.1,0.1,legth=100)) z <- matri(appl(epad.grid(th),1,rss), legth(th[[1]])) image(th$beta0,th$beta1,-z) cotour(th$beta0,th$beta1,z,levels=c(0:9,(1:10)*10),add=t) poits(be[1],be[],col=4) dev.copeps(file="ru003-cp.eps") > source("ru003.r") : [1] : [1] ɛ: β0, β1:,. β0, β1 i = 1,..., i = β0 + β1i + ɛi β0, β1 th$beta ɛ i = (β0 + β1i i) ru003-s th$beta0 ru003-cp 3 4

2 .3 ˆβ1 = S, S ˆβ0 = ȳ ˆβ1 S = (β0 + β1i i) β0 β1 S = (i β0 β1i) = 0 β0 S = (i β0 β1i)i = 0 β1 β0, β1 ȳ β0 β1 = 0, (i ȳ β1(i ))i = 0 {(i ȳ)(i ) β1(i ))(i ) = 0 = ˆβ0 + ˆβ1 S β1s = 0 ȳ = ˆβ0 + ˆβ1 # ru0033.r ( ) dat <- read.table("dat0001.tt") (47 ) <- dat[,1]; <- dat[,] rss <- fuctio(be) sum( (be[1]+be[]* - )^ ) beta1 <- cov(,)/var() # beta1 beta0 <- mea() - beta1*mea() # beta0 be <- c(beta0,beta1); cat(": "); prit(be) # cat(" : "); prit(rss(be)) # plot(,,pch=16) # plot(dat,pch=16) ablie(be[1],be[],col=) ablie(v=mea(),col=3,lt=) # () ablie(h=mea(),col=3,lt=) # () dev.copeps(file="ru0033-s.eps") > source("ru0033.r") : [1] : [1] ru0033-s.4 R lsfit() lm() R # ru0034.r ( ) dat <- read.table("dat0001.tt") (47 ) <- dat[,1]; <- dat[,] cat("# lsfit \") fit1 <- lsfit(,) QR prit(fit1$coef) cat("# lm \") fit <- lm( ~, data.frame(,)) prit(fit$coef) fit3 <- lm(shushou ~ Gakureki,dat) prit(fit3$coef) cat("# lsfit \") ls.prit(fit1) 5 6 cat("# lm \") prit(summar(fit3)) > source("ru0034.r") # lsfit Itercept X # lm Gakureki # lsfit Residual Stadard Error=0.09 R-Square=0.534 F-statistic (df=1, 45)= p-value=0 Estimate Std.Err t-value Pr(> t ) Itercept X # ru0035.r fuc0035 <- fuctio(a) { plot(0,0,lab="",lab="",lim=c(0,1),lim=c(0,1),tpe="") a <- locator(tpe="p") <- a$; <- a$ beta1 <- cov(,)/var() # beta1 beta0 <- mea() - beta1*mea() # beta0 be <- c(beta0,beta1); prit(be) # ablie(be[1],be[],col=) ablie(v=mea(),col=3,lt=) # () ablie(h=mea(),col=3,lt=) # () tet(0.0,0.95,pos=4, paste("beta0=",sigif(be[1],5),", beta1=",sigif(be[],5), ", =",legth(),sep="")) # beta dev.copeps(file=paste("ru0035-s",a,".eps",sep="")) ivisible(a) prit # lm Call: lm(formula = Shushou ~ Gakureki, data = dat) > source("ru0035.r") > fuc0035(1) [1] > fuc0035() [1] Residuals: Mi 1Q Media 3Q Ma Coefficiets: < e-16 *** Gakureki e-09 *** --- Sigif. codes: 0 *** ** 0.01 * Residual stadard error: o 45 degrees of freedom Multiple R-Squared: 0.534,Adjusted R-squared: 0.5 F-statistic: 51.4 o 1 ad 45 DF, p-value: 5.943e-09 beta0= , beta1=0.867, =9 ru0035-s1 beta0= , beta1= , =4 ru0035-s 7 8

3 .6 ɛ 0 σ ɛ N(0, σ ) ɛ1,..., ɛ (1,..., ) N(0, σ ) (1, 1),..., (, ) f(1, 1,...,, β0, β1, σ) = f(1,..., )f(1,..., 1,...,, β0, β1, σ) 1,..., 1,..., ɛi = i β0 β1i f(1,..., 1,...,, β0, β1, σ) = { 1 πσ ep ( ) (i β0 β1i) σ θ = (β0, β1, σ) l(β0, β1, σ) = ) (i β0 β1i) log(πσ ) ( σ β0 β1 l l = 0, β0 l = 0 β1 ˆβ1 = S, S ˆβ0 = ȳ ˆβ1 σ l (σ ) = ) σ + (i β0 β1i) ( σ 4 0 ˆσ = 1 (i ˆβ0 ˆβ1i) 9.7 i β0 + β1i ŷi = ˆβ0 + ˆβ1i (residual) ei = i ŷi ɛ Sɛ = 1 e i ˆσ = 1 e i = # ru0036.r fuc0036 <- fuctio(,) { beta1 <- cov(,)/var() # beta1 beta0 <- mea() - beta1*mea() # beta0 beta <- c(beta0,beta1) pred <- beta0 + beta1 * # resid <- - pred # se <- sqrt(sum(resid^)/(legth()-)) list(beta=beta,=,=,pred=pred,resid=resid,se=se) dat <- read.table("dat0001.tt") (47 ) fit <- fuc0036(dat[,1],dat[,]) plot(fit$pred,fit$) ablie(0,1,col=,lt=) tet(1.,1.8,pos=4,paste("beta0=",sigif(fit$beta[1],5), ", beta1=",sigif(fit$beta[],5),", =",legth(fit$),sep="")) dev.copeps(file="ru0036-s1.eps") plot(fit$pred,fit$resid) ablie(h=0,col=,lt=) # e=0 ablie(h=fit$se*,col=3,lt=3) # e=*se ablie(h=-fit$se*,col=3,lt=3) # e=-*se 10 tet(1.,0.3,pos=4,paste("se=",sigif(fit$se,5),sep="")) dev.copeps(file="ru0036-s.eps") source("ru0044.r") # drawhist drawhist(fit$resid,0,"resid","ru0036-") fit$ beta0=1.745, beta1= , = fit$pred ru0036-s1 resid mea=.0787e 16, sd= ru0036-resid fit$resid se= fit$pred ru0036-s (i, i) (X, Y ) (i, i) 1. (1, 1),..., (, ) ( 1, 1),..., (, ). ( ˆβ 0, ˆβ 1, Sɛ, Sɛ ) 3. 1, (10000 ) # ru0050.r # fit <- fuc0036(,) fuc0050 <- fuctio(f) # fit ulist(list(beta0=f$beta[1],beta1=f$beta[],sigma=f$se,sigma=f$se^)) cat("\ \"); prit(fuc0050(fit)) a <- fuc0036([i],[i]); fuc0050(a) <- legth() b < cat("\ \"); prit(boot1(1:)) simi <- matri(0,,b) for(j i 1:b) simi[,j] <- sample(1:,replace=t) simt <- appl(simi,,boot1) cat("\ \") a <- appl(simt,1,fuctio() ulist(list(mea=mea(),sd=sd()))) for(k i rowames(simt)) drawhist(simt[k,],0,k,"ru0050-") prit(summar(lm(~,data.frame(,)))$coef) ( ˆβ0, ˆβ1, Sɛ) (1, 1),..., (, ) 11 > librar(mass) > source("ru0044.r") > dat <- read.table("dat0001.tt") (47 ) > <- dat[,1]; <- dat[,] > source("ru0050.r") 1

4 mea sd beta0 mea= , sd= ru0050-beta e e-09 beta1 mea= , sd= ru0050-beta1 0 5 sigma mea= , sd= ru0050-sigma sigma mea= , sd= ru0050-sigma lm() Std.Error lm() 1,..., ɛ lm() lm() 3. ( ˆβ0, ˆβ1, Sɛ, S ɛ ) (1, 1),..., (, ) (i, i) i i Y 1,..., 1,..., Y # ru0051.r 1. ˆβ0, ˆβ1. ŷi = ˆβ0 + ˆβ1i ei = i ŷi 3. e1,..., e e 1,..., e i i = ŷi + e i 4. (1, 1),..., (, ) ( ˆβ 0, ˆβ 1, Sɛ, Sɛ ) 5. 3,4 (10000 ) # fit <- fuc0036(,) cat("\ \"); prit(fuc0050(fit)) pred <- fit$pred # resid <- fit$resid # i <- pred + resid[i] a <- fuc0036(,i); fuc0050(a) <- legth() b < cat("\ \"); prit(boot1(1:)) simi <- matri(0,,b) for(j i 1:b) simi[,j] <- sample(1:,replace=t) simt <- appl(simi,,boot1) cat("\ \") a <- appl(simt,1,fuctio() ulist(list(mea=mea(),sd=sd()))) for(k i rowames(simt)) drawhist(simt[k,],0,k,"ru0051-") prit(summar(lm(~,data.frame(,)))$coef) > source("ru0051.r") mea sd beta0 mea= 1.74, sd= ru0051-beta e e-09 beta1 mea= 0.088, sd= ru0051-beta1 0 5 sigma mea= , sd= ru0051-sigma sigma mea= , sd= ru0051-sigma lm() E(S ) < S lm() Std.Error S S 3.3 ( ˆβ0, ˆβ1, Sɛ, S ɛ ) hat() hii, i = 1,..., ri, i = 1,..., ri = ei/ 1 hii 15 # ru005.r # fit <- fuc0036(,) cat("\ \"); prit(fuc0050(fit)) pred <- fit$pred # resid <- fit$resid/sqrt(1-hat()) plot(,1/sqrt(1-hat())) dev.copeps(file="ru005-hat.eps") i <- pred + resid[i] a <- fuc0036(,i); fuc0050(a) <- legth() b < cat("\ \"); prit(boot1(1:)) simi <- matri(0,,b) for(j i 1:b) simi[,j] <- sample(1:,replace=t) simt <- appl(simi,,boot1) cat("\ \") a <- appl(simt,1,fuctio() ulist(list(mea=mea(),sd=sd()))) for(k i rowames(simt)) drawhist(simt[k,],0,k,"ru005-") prit(summar(lm(~,data.frame(,)))$coef) > source("ru005.r") mea sd

5 e e > fuc0035(3) [1] > fuc0035(4) [1] /sqrt(1 hat()) ru005-hat beta0= , beta1=0.891, =10 beta0= , beta1= , =11 beta0 beta1 sigma sigma ru0035-s3 4. ru0035-s4 mea= 1.74, sd= ru005-beta0 mea= , sd= ru005-beta1 mea= , sd= ru005-sigma mea= , sd= ru005-sigma (β0 + β1i i) 1/ 1 hii, i = 1,..., 1 ei ei lm() E(S ) S S E(S ) < S S = S ( ) 17 wi (β0 + β1i i) w1 = = w = 1 (i, i) wi wi = 0 wi 1. ( ˆβ0, ˆβ1) ˆσ. ŷi = ˆβ0 + ˆβ1i, i = 1,..., ei = i ŷi ˆσ ei wi ei ˆσ wi 1 4. ˆβ0, ˆβ1, ˆσ 5.,3,4 4.3 R MASS rlm lqs Veables ad Riple Moder Applied Statistics with S fourth editio 156 lm rlm lqs a <- locator(tpe="p") fit <- rkaiki(a$,a$) cat("\\") prit(t(sappl(fit,fuctio(f) f$coef))) dev.copeps(file=paste("ru0053-s",a,".eps",sep="")) ivisible(a) prit > source("ru0053.r") > fuc0053(1) LM M MM LTS LMS > fuc0053() # ru0053.r librar(mass) # rlm,lqs MASS rkaiki <- fuctio(,,led=c(mi(),ma()),pl=t) { a <- data.frame(,) fit <- list() fit[[1]] <- lm( ~, a) LM M MM LTS LMS fit[[]] <- rlm( ~, a) # M (method="m") fit[[3]] <- rlm( ~, a, method="mm") # MM fit[[4]] <- lqs( ~, a) LTS (method="lts") fit[[5]] <- lqs( ~, a, method="lms") LMS ames(fit) <- c("lm","m","mm","lts","lms") if(pl) { plot(,) # for(i i 1:5) ablie(fit[[i]],lt=i,col=i) leged(led[1],led[],ames(fit),lt=1:5,col=1:5,bt="") # ivisible(fit) LM M MM LTS LMS LM M MM LTS LMS fuc0053 <- fuctio(a) { plot(0,0,lab="",lab="",lim=c(0,1),lim=c(0,1),tpe="") ru0053-s1 ru0053-s 19 0

6 rlm lm (robust) rlm help(rlm) 1. rlm(...,method="m") M rlm Huber (ˆσ) MAD. rlm(...,method="mm") MM M Tuke s biweight (breakdow poit 0.5 efficiet, > dat <- read.table("dat0001.tt") (47 ) > <- dat[,1]; <- dat[,] > source("ru0054.r") LM M MM LTS LMS lqs 1. lqs(...,method="lts") LTS (least trimmed squares) lqs e 1,..., e ( α 0.5 ). lqs(...,method="lms") LMS (least media of squares) LM M MM LTS LMS e 1,..., e ru0054-s MM LTS # ru0054.r # fit <- rkaiki(,) cat("\\") prit(t(sappl(fit,fuctio(f) f$coef))) dev.copeps(file="ru0054-s.eps") for(i i c(,3)) { resid <- fit[[i]]$resid; weight <- fit[[i]]$w plot(resid,weight) dev.copeps(file=paste("ru0054-w",i,".eps",sep="")) weight resid ru0054-w weight resid ru0054-w ˆβ0, ˆβ1 (1, 1),..., (, ) (i, i) (X, Y ) (i, i) dat <- read.table("dat0001.tt") (47 ) <- dat[,1]; <- dat[,] source("ru0055.r") > source("ru0055bat.r") # ru0055.r # fit <- rkaiki(,,pl=f) fuc0055 <- fuctio(fit) { # fit a <- sappl(fit,fuctio(f) f$coef) a <- as.vector(outer(c("beta0","beta1"), dimames(a)[[]],paste,sep="-")) a <- as.vector(a); ames(a) <- a a cat("\\"); prit(fuc0055(fit)) fit <- rkaiki([i],[i],pl=f); fuc0055(fit) <- legth() b < cat("\ \"); prit(boot1(1:)) simi <- matri(0,,b) prit(date()) for(j i 1:b) simi[,j] <- sample(1:,replace=t) prit(date()) simt <- appl(simi,,boot1) prit(date()) cat("\ \") a <- appl(simt,1,fuctio() ulist(list(mea=mea(),sd=sd()))) for(k i rowames(simt)) drawhist(simt[k,],0,k,"ru0055-") prit(summar(lm(~,data.frame(,)))$coef) # ru0055bat.r source("ru0044.r") source("ru0053.r") 3 beta0-lm beta1-lm beta0-m beta1-m beta0-mm beta1-mm beta0-lts beta1-lts beta0-lms beta1-lms beta0-lm beta1-lm beta0-m beta1-m beta0-mm beta1-mm beta0-lts beta1-lts beta0-lms beta1-lms [1] "Tue Sep 8 17:14:9 004" [1] "Tue Sep 8 17:14:30 004" [1] "Tue Sep 8 17:34:36 004" beta0-lm beta1-lm beta0-m beta1-m beta0-mm beta1-mm mea sd beta0-lts beta1-lts beta0-lms beta1-lms mea sd e e-09 There were 18 warigs (use warigs() to see them) beta0 LM beta0 M beta0 MM beta0 LTS beta0 LMS mea= , sd= mea= 1.746, sd= mea= , sd= mea= , sd= mea= 1.743, sd= ru0055-beta0-lm ru0055-beta0-m ru0055-beta0-mmru0055-beta0-ltsru0055-beta0-lms 4

7 beta1 LM beta1 M beta1 MM beta1 LTS beta1 LMS mea= , sd= mea= , sd= mea= , sd= mea= , sd= mea= , sd= ru0055-beta1-lm ru0055-beta1-m ru0055-beta1-mmru0055-beta1-ltsru0055-beta1-lms LMS ˆβ0, ˆβ1 lqs(..., method = lms ) media() optim() ( ) 0 LM M MM << LTS LMS 1. lm. rlm 3. lqs MM ( ) = 0 50 LM, M, MM, LTS, LMS ( ˆβ0, ˆβ1) LM LM, M, MM, LTS, LMS ) 5 6

2. S 2 ɛ 3. ˆβ S 2 ɛ (n p 1)S 2 ɛ χ 2 n p 1 Z N(0, 1) S 2 χ 2 n T = Z/ S 2 /n n t- Z T = S2 /n t- n ( ) (n+1)/2 Γ((n + 1)/2) f(t) = 1 + t2 nπγ(n/2) n

2. S 2 ɛ 3. ˆβ S 2 ɛ (n p 1)S 2 ɛ χ 2 n p 1 Z N(0, 1) S 2 χ 2 n T = Z/ S 2 /n n t- Z T = S2 /n t- n ( ) (n+1)/2 Γ((n + 1)/2) f(t) = 1 + t2 nπγ(n/2) n Copright (c) 2004,2005 Hidetoshi Shimodaira 1. t- 2. 3. F 4. 1 t- 1.1 = β0 + β11 + + βpp + ɛ = Xβ + ɛ 2005-04-13 17:37:08 shimo S 2 ɛ ˆV ( ˆβi) = S 2 ɛ aii ŝe( ˆβi) = t-p- ti = ˆβi ŝe( ˆβi), ˆV ( ˆβi)