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

Size: px

Start display at page:

Download "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"

けいしょうやまがた
5 years ago
Views:

1 Copright (c) 2004,2005 Hidetoshi Shimodaira 1. t F 4. 1 t- 1.1 = β0 + β βpp + ɛ = Xβ + ɛ :37:08 shimo S 2 ɛ ˆV ( ˆβi) = S 2 ɛ aii ŝe( ˆβi) = t-p- ti = ˆβi ŝe( ˆβi), ˆV ( ˆβi) = Sɛ aii pi = 2 PrT > ti T n p 1 t- βi = 0 ti t- βi = 0 βi 0 α α = pi < α βi 0 pi α βi = 0 βi = 0 pi [0, 1] pi < α βi 0 α β ˆβ = (X X) 1 X ɛ Nn(0, In) ˆβ Np+1(β, (X X) 1 ) (p + 1) (p + 1) A A = (aij) = (X X) 1, i, j = 0,..., p ˆβi N(βi, aii) V ( ˆβi) = aii ˆ ɛ Sɛ 2 1 n = e 2 i = e 2 n p 1 n p 1 β0, β1,..., βp p+1 n (p+1) 1 α α βi = 0 βi 0 βi 0 βi = 0 t- n p 1 t- # run0073.r # t- # fit <- lm(~.,data.frame(,)) # be <- coef(summar(fit)) # cat("# t-p- \n"); print(be) tval <- be[,"t value"] # t- m <- ma(abs(tval))+0.5; 0 <- seq(-m,m,len=300) # plot(0,dt(0,df.residual(fit)),tpe="l", lab="t-statistic",lab="densit") # t- rug(tval); tet(tval,0.01,names(tval),srt=90,adj=0) dev.cop2eps(file="run0073-d.eps") pv0 <- 2*pt(abs(0),df.residual(fit),lower.tail=F) # p- plot(0,pv0,tpe="l",log="",lab="t-statistic",lab="p-value") 2 abline(h=c(0.01,0.05),lt=3) rug(tval); tet(tval,pv0[1]*1.1,names(tval),srt=90,adj=0) dev.cop2eps(file="run0073-t.eps") > dat <- read.table("dat0002.tt") # (47 10 ) > <- dat[,-10]; <- dat[,10] > source("run0073.r") # t-p- Estimate Std. Error t value Pr(> t ) (Intercept) Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin n Z1,..., Zn N(0, 1) n n Zi 2 χ 2 n n ɛ 2 i σ 2 χ2 n ei N(0, (1 hii)) ei/σ N(0, 1) n n p 1 e 2 i χ2 n e = (In H)ɛ e 2 = ɛ (In H) 2 ɛ = ɛ (In H)ɛ = ɛ 2 ɛ Hɛ H Span(0,..., p) In H n n U = (u1, u2,..., un) H = (u1,..., up+1)(u1,..., up+1) densit Kfufu Ninzu Tandoku Tomo Kaku Konin X65Sai Ktan (Intercept) Zouka p value 5e 04 5e 03 5e 02 5e 01 Kfufu Ninzu Tandoku Tomo Kaku Konin X65Sai Ktan (Intercept) Zouka In H = (up+2,..., un)(up+2,..., un) zi = u iɛ/σ, i = 1,..., n N(0, 1) ɛ Hɛ/ = z z 2 p+1 χ 2 p+1 e 2 / = z 2 p z 2 n χ 2 n p 1 ˆβ z1,..., zp+1 zp+2,..., zn ˆβ S 2 ɛ = e 2 /(n p 1) t statistic run0073-d t statistic run0073-t 1. ˆβ 1.2 t- ˆβ Np+1(β, (X X) 1 ) (p + 1) (p + 1) A n p 1 n e 2 i σ 2 χ2 n p 1 3 A = (aij) = (X X) 1, i, j = 0,..., p ˆβi N(βi, aii) 4

2 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 ˆβi βi σ aii R t- βi = 0 ˆβi βi σ N(0, 1) aii (n p 1)S 2 ɛ χ 2 n p 1 / Sɛ 2 σ = ˆβi βi t- 2 n p 1 aiisɛ ˆβi βi ŝe( ˆβi) t- n p 1 ti = ˆβi ŝe( ˆβi) ti t- n p 1 βi > 0 t- ti βi < 0 t- ti 1.3 # run0074.r # source("run0044.r") # drawhist 5 #n <- 11 # #be <- c(0,1) # (beta0,beta1) # <- seq(0,1,len=n) # #filename <- "run0074-" p <- length(be)-1 # p+1 X <- cbind(1,) # colnames(x) <- names(be) <- c("beta0","beta1") A <- solve(t(x) %*% X) # A=(X X)^-1 B <- A %*% t(x) # B = (X X)^-1 X sqra <- sqrt(diag(a)) # A func0074 <- function() be <- B %*% # pred <- X %*% be # resid <- -pred # se2 <- sum(resid^2)/(n-p-1) # sigma^2 se <- sqrt(se2) # sigma tval <- be/(se*sqra) # t- pval <- 2*pt(abs(tval),n-p-1,lower.tail=F) # p- list(be=be,se2=se2,tval=tval,pval=pval) 0 <- X %*% be # ) sigma0 <- 0.3 # cat("# start simulation: ",date(),"\n") b < # sim <- matri(0,n,b) # simbe <- matri(0,length(be),b) # simse2 <- rep(0,b) # se2 simtval <- matri(0,length(be),b) # t simpval <- matri(0,length(be),b) # p- for(i in 1:b) sim[,i] <- 0 + rnorm(n,mean=0,sd=sigma0) fit <- func0074(sim[,i]) simbe[,i] <- fit$be; simse2[i] <- fit$se2 simtval[,i] <- fit$tval; simpval[,i] <- fit$pval cat("# end simulation: ",date(),"\n") cat("# \n"); print(func0074(sim)) cat("# pval0 < 0.05 = ",sum(simpval[1,]<0.05),"\n") cat("# pval1 < 0.05 = ",sum(simpval[2,]<0.05),"\n") if(!is.null(filename)) plot(,0); abline(be) dev.cop2eps(file=paste(filename,"s0.eps",sep="")) 6 plot(,sim); abline(simbe) dev.cop2eps(file=paste(filename,"s1.eps",sep="")) plot(simbe[1,],simbe[2,],pch=".",lab="beta0",lab="beta1") dev.cop2eps(file=paste(filename,"th1.eps",sep="")) plot(simse2,simbe[2,],pch=".",lab="se2",lab="beta1") dev.cop2eps(file=paste(filename,"th2.eps",sep="")) for(i in 1:2) drawhist(simtval[i,],30,paste("tval",i-1,sep="")) t0 <- seq(min(simtval[i,]),ma(simtval[i,]),len=300) lines(t0,dt(t0,n-p-1),col=4,lt=2) dev.cop2eps(file=paste(filename,"tval",i-1,".eps",sep="")) drawhist(simpval[i,],20,paste("pval",i-1,sep=""),filename) 0 run0074-s0 sim[, 1] run0074-s1 > n <- 11 # > be <- c(0,1) # (beta0,beta1) > <- seq(0,1,len=n) # > filename <- "run0074-" > source("run0074.r") # start simulation: Wed Oct 6 11:45: # end simulation: Wed Oct 6 11:45: # $be beta beta $se2 [1] $tval beta beta $pval beta beta # pval0 < 0.05 = 504 # pval1 < 0.05 = 8748 beta beta0 run0074-th1 tval beta se2 run0074-th2 tval mean= , sd= run0074-tval0 mean= 3.83, sd= run0074-tval1 7 8

3 pval0 mean= , sd= run0074-pval pval1 mean= , sd= run0074-pval1 b = Pentium 1G vmware linu ) v = w0β0 + w1β1 + + wpβp ˆv = w0 ˆβ0 + w1 ˆβ1 + + wp ˆβp w = (w0, w1,..., wp) ˆv wj = 1, wk = 0, k j ˆv = ˆβj w0 = i0,... wp = ip ˆv = ŷi v = w β, ˆv = w ˆβ ˆβ Np+1(β, (X X) 1 ) ˆv N(w β, w (X X) 1 w) 2.2 S 2 ɛ v = w β, ˆv N(v, v) v = w (X X) 1 w ˆσ v 2 = Sɛ 2 w (X X) 1 w = Sɛ 2 σv 2 σ = e 2 / 2 n p 1 σ2 v (n p 1) ˆσ2 v χ 2 σv 2 n p 1 ˆv v σv / ˆσ v 2 t- σv 2 n p 1 ˆv v t- n p 1 R n p 1 t- T t pt(t) PrT t = pt(t,n-p-1) PrT > t = 1-pt(t,n-p-1) = pt(t,n-p-1,lower.tail=f) t- βi = 0 v = v0 v v0 p- ( ( )) p- (v0) = Pr T > ˆv v0 ˆv v0 = 2 1 pt v = v0 p- (v0) p- (v0) < α v = v0 v v0 v ˆv v σv / ˆσ v 2 t- σv 2 n p 1 p- (v) [0, 1] Prp- (v) < α = α v = v0 p- (v0) < α v = v0 v v0 v0 (1 α) = v0 p- (v0) α 1 α α = v 1 α Prv (1 α) = 1 α v (1 α)p- (v) α p- (v) [0,1] 1 α a = pt(b) b = qt(a) n p 1 t- T t a t qt(a) R PrT qt(a) = a qt(a) = qt(a,n-p-1) 0 < a < 1 ( ( )) ˆv v0 2 1 pt α ( ) ˆv v0 pt 1 α/2 ˆv v0 qt(1 α/2) ˆv qt(1 α/2) v0 ˆv + qt(1 α/2) (1 α) = [ˆv qt(1 α/2), 11 ˆv + qt(1 α/2)] 2.4 βj w wj = 1, wk = 0, k j ˆv = ˆβj V (ˆv) ˆσ v 2 = Sɛ 2 w (X X) 1 w = Sɛ 2 ajj βj (1 α) = [ ˆβj Sɛ ajjqt(1 α/2), ˆβj + Sɛ ajjqt(1 α/2)] 2.5 = (1, 1,..., p) ŷ = ˆβ = β v = w = V (ˆv) ˆσ v 2 = Sɛ 2 w (X X) 1 w = Sɛ 2 (X X) 1 (1 α) = [ ˆβ Sɛ (X X) 1 qt(1 α/2), ˆβ+Sɛ (X X) 1 qt(1 α/2)] # run0075.r # # = = p= # a <- func0075a(,p); b <- func0075b(,0.05,a) pow <- function(a,p) a^(0:p) # c(a^0,a^1,...,a^p) calcx <- function(,p) # X X <- matri(0,length(),p+1) for(i in 1:length()) X[i,] <- pow([i],p) X calcq0 <- function(n,p,alpha) qt(1-alpha/2,n-p-1) # calcq1 <- function(n,p,alpha) sqrt((p+1)*qf(1-alpha,p+1,n-p-1)) # func0075a <- function(,p) # X <- calcx(,p) # colnames(x) <- paste("beta",0:p,sep="") A <- solve(t(x) %*% X) # A=(X X)^-1 B <- A %*% t(x) # B = (X X)^-1 X sqra <- sqrt(diag(a)) # A 0 <- seq(min(),ma(),len=300) #

4 X0 <- calcx(0,p) # 0 sqrxax <- appl(x0,1,function() sqrt(t() %*% A %*% )) # sqrt( A) list(x=x,a=a,b=b,sqra=sqra,0=0,x0=x0,sqrxax=sqrxax) func0075b <- function(,alpha,a,calcq=calcq0) # n <- nrow(a$x); p <- ncol(a$x)-1 q0 <- calcq(n,p,alpha) be <- a$b %*% # pred <- a$x %*% be # resid <- -pred # rss <- sum(resid^2) # se2 <- rss/(n-p-1) # sigma^2 se <- sqrt(se2) # sigma bese <- se*a$sqra # rsq <- 1-rss/sum((-mean())^2) # tval <- be/(se*a$sqra) # t- pval <- 2*pt(abs(tval),n-p-1,lower.tail=F) # p- beconf <- cbind(be-q0*se*a$sqra,be+q0*se*a$sqra) # pred0 <- a$x0 %*% be # (0) pred0conf <- cbind(pred0-q0*se*a$sqrxax,pred0+q0*se*a$sqrxax) list(be=be,bese=bese,se=se,rss=rss,rsq=rsq,tval=tval,pval=pval, beconf=beconf,pred0=pred0,pred0conf=pred0conf) func0075c <- function(,,a,b,col=2,lt=2,add=f) if(!add) plot(,) lines(a$0,b$pred0) lines(a$0,b$pred0conf,col=col,lt=lt) lines(a$0,b$pred0conf[,2],col=col,lt=lt) coef <- cbind(b$be,b$bese,b$tval,b$pval,b$beconf) colnames(coef) <- c("estimate","stderr", "t-value","p-value","lower","upper") invisible(list(coef=coef,se=b$se,rss=b$rss,rsq=b$rsq)) # run0076.r # p # = = p= source("run0075.r") for(i in 0:p) cat("# =",i,"\n") a <- func0075a(,i) c1 <- func0075c(,,a,func0075b(,0.05,a)) c2 <- func0075c(,,a,func0075b(,0.01,a),col=3,add=t) colnames(c1$coef)[5:6] <- c("lo05","up05") colnames(c2$coef)[5:6] <- c("lo01","up01") coef <- cbind(c1$coef,c2$coef[,5:6,drop=f]) cat("rss=",c1$rss,", RSQ=",c1$rsq,"\n") print(round(coef,3)) dev.cop2eps(file=paste("run0076-s",i,".eps",sep="")) > dat <- read.table("dat0001.tt") # (47 2 ) > <- dat/100; <- dat[,2] > p <- 3 > source("run0076.r") # = 0 RSS= , RSQ= 0 Estimate StdErr t-value p-value Lo05 Up05 Lo01 Up01 beta # = 1 RSS= , RSQ= Estimate StdErr t-value p-value Lo05 Up05 Lo01 Up01 beta beta # = 2 RSS= , RSQ= Estimate StdErr t-value p-value Lo05 Up05 Lo01 Up01 beta beta beta # = 3 RSS= , RSQ= Estimate StdErr t-value p-value Lo05 Up05 Lo01 Up01 beta beta beta beta p = 2 β0 = 0 β0 0 β1 = β1 = 0 β2 = 0 p = 0 4. p = 3 β0 0 β1 = β2 = β3 = 0 p = 0 t- p = 1 p = 2 β2 = 0 or β2 0 3 F run0076-s0 run0076-s1 3.1 p = β0 + β βpp + ɛ = Xβ + ɛ k p k = β0 + β βkk + ɛ = X (k) β (k) + ɛ run0076-s2 run0076-s3 X (k) = (0, 1,..., k), β (k) = (β0, β1,..., βk) p = 0, 1, 2, 3 α = 0.05 α = 0.01 RSS = e 2 R 2 p = 1 t 1. p = 0 β0 = 0 β0 0 p 0 X ( k) = (k+1, k+2,..., p), β ( k) = (βk+1, βk+2,..., βp) [ ] X = (X (k), X ( k) ), β = β (k) β ( k) β (k) ˆβ (k) = (X (k) X (k) ) 1 X (k) i = 0,..., k ˆβ (k) i ˆβ i 2. p = 1 β0 = 0 β0 0 β1 0 p

5 p βk+1 = βk+2 = = βp = 0 k k k Span(0,..., k) Span(0,..., p) 3.2 F (k) (k) (k) H (k) = X (k) (X (k) X (k) ) 1 X (k) ŷ (k) = H (k) e (k) = ŷ (k) = (In H (k) ) RSS (k) = e (k) 2 RSS (k) χ 2 n k 1 RSS (k) RSS (p) = e(k) 2 e(p) 2 χ 2 p k RSS (p) χ 2 n p 1 (k) RSS(k) σ 2 (k) F = (RSS(k) RSS (p) )/(p k) RSS (p) /(n p 1) p k n p 1 p k, n p 1 F F Fp k,n p 1 17 F R pf X m, n F PrX = pf(, m, n) PrX > = 1 pf(, m, n) = pf(, m, n, lower.tail = F) pf qf PrX qf(a, m, n) = a # run0077.r # fkentei <- function(fitk,fitp) rssk <- sum(resid(fitk)^2) # RSS^(k) rssp <- sum(resid(fitp)^2) # RSS^(p) dfk <- df.residual(fitk) # n-k-1 dfp <- df.residual(fitp) # n-p-1 bunshi <- (rssk - rssp)/(dfk-dfp) # (RSS^(k)-RSS^(p))/(p-k) bunbo <- rssp/dfp # RSS^(p) / (n-p-1) fvalue <- bunshi/bunbo # F pvalue <- pf(fvalue,dfk-dfp,dfp,lower.tail=f) list(fvalue=fvalue,df1=dfk-dfp,df2=dfp,pvalue=pvalue) > source("run0077.r") > # > dat <- read.table("dat0001.tt") # (47 2 ) > <- dat/100; <- dat[,2] > fit0 <- lm( ~ 1, data.frame(,)) # k=0 > fit1 <- lm( ~, data.frame(,)) # k=1 > fit2 <- lm( ~ + I(^2), data.frame(,)) # k=2 > summar(fit2) # k=2 Call: lm(formula = ~ + I(^2), data = data.frame(, )) Residuals: Min 1Q Median 3Q Ma Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) <2e-16 *** I(^2) Signif. codes: 0 *** ** 0.01 * Residual standard error: on 44 degrees of freedom Multiple R-Squared: ,Adjusted R-squared: F-statistic: on 2 and 44 DF, p-value: 4.7e-08 > unlist(fkentei(fit0,fit2)) # F! fvalue df1 df2 pvalue e e e e-08 > unlist(fkentei(fit1,fit2)) # t (^2 )! fvalue df1 df2 pvalue > sqrt( ) # t (^2 ) [1] F (k) (p) 1. k = 0 = (p) F p-value < α p 2. k = p 1 p βp = 0 or βp 0 βp t- ( ) t 2 = F F t- 3.3 F F (k) (p) l1(θ1) l2(θ2) θ1 θ2 p1 = dim θ1,p2 = dim θ2 θ1 = (β0, β1,..., βk, σ) θ2 = (β0, β1,..., βp, σ) p1 = k + 2 p2 = p ˆθ1, ˆθ2 l1(ˆθ1) l2(ˆθ2) l1(ˆθ1) = l( ˆβ0, ˆβ1,..., ˆβk, ˆσ) = n 2 l2(ˆθ1) = l( ˆβ0, ˆβ1,..., ˆβp, ˆσ) = n log(2π RSS(k) ) n 1 + log(2π RSS(p) ) n 1 p2 p1 n # run0081.r # 2 (l2(ˆθ1) l1(ˆθ2)) χ 2 p2 p1 2 (l2(ˆθ1) l1(ˆθ2)) = n log RSS (k) n log RSS (p) χ 2 p k udohikentei <- function(fitk,fitp) rssk <- sum(resid(fitk)^2) # RSS^(k) rssp <- sum(resid(fitp)^2) # RSS^(p) dfk <- df.residual(fitk) # n-k-1 dfp <- df.residual(fitp) # n-p-1 n <- length(resid(fitk)) # n chisqvalue <- n*(log(rssk)-log(rssp)) # pvalue <- pchisq(chisqvalue,dfk-dfp,lower.tail=f) list(chisqvalue=chisqvalue,df=dfk-dfp,n=n,pvalue=pvalue) > dat <- read.table("dat0001.tt") # (47 2 ) > <- dat/100; <- dat[,2] > fit0 <- lm( ~ 1, data.frame(,)) # k=0 > fit1 <- lm( ~, data.frame(,)) # k=1 > fit2 <- lm( ~ + I(^2), data.frame(,)) # k=2 > source("run0081.r") > unlist(fkentei(fit0,fit2)) # F fvalue df1 df2 pvalue e e e e-08 > unlist(udohikentei(fit0,fit2)) # chisqvalue df n pvalue e e e e-08 > unlist(fkentei(fit1,fit2)) # F 20

6 fvalue df1 df2 pvalue > unlist(udohikentei(fit1,fit2)) # chisqvalue df n pvalue QR QR X = QR X n (p + 1) Q n (p + 1) R (p + 1) (p + 1) Q Q Q = Ip+1 R X X = R Q QR = R R (X X) 1 = (R R) 1 = R 1 R 1 QR ˆβ ˆβ = (X X) 1 X = R 1 R 1 R Q = R 1 Q Q ˆβ = R 1 (Q ) R 1 R Q R ˆβ R 1 > X <- matri(c(1^(1:5),2^(1:5),3^(1:5)),5,3) > X [,2] [,3] [1,] [2,] [3,] [4,] [5,] > QR <- qr(x) > R <- qr.r(qr) > Q <- qr.q(qr) > R [,2] [,3] [1,] [2,] [3,] > Q [,2] [,3] 21 [1,] [2,] [3,] [4,] [5,] > Q %*% R [,2] [,3] [1,] [2,] [3,] [4,] [5,] > t(q) %*% Q [,2] [,3] [1,] e e e-17 [2,] e e e-17 [3,] e e e F U = (u1, u2,..., un) H (k) = (u1,..., uk+1)(u1,..., uk+1), k = 0,..., p 0, 1,..., p QR X = QR ui = q i, i = 1,..., p + 1 up+2,..., un (p) = Xβ + ɛ ɛ Nn(0, In) = X (k) β (k) + X ( k) β ( k) + ɛ (k) e (k) = (In H (k) ) = (In H (k) )(X ( k) β ( k) + ɛ) β ( k) = 0 e (k) = (In H (k) ) = (In H (k) )ɛ = (uk+2,..., un)(uk+2,..., un) ɛ 22 zi = u ɛ/σ, i = 1,..., n (k) e (k) 2 = z 2 k z 2 n RSS (k) = e(k) 2 = z 2 k zn 2 χ 2 n k 1 RSS (p) = e(p) 2 = z 2 p zn 2 χ 2 n p 1 RSS (k) RSS (p) = e(k) 2 e(p) 2 = z 2 k zp+1 2 χ 2 p k RSS (p) RSS (k) RSS (p) zi ˆβ ˆβ Np+1(β, (X X) 1 ) F = z 2 /(p + 1) Sɛ 2 / Fp+1,n p 1 p + 1, n p 1 F z 2 = X( ˆβ β) 2 / F = X( ˆβ β) 2 (p + 1)S 2 ɛ ˆβ or 1 α (1 α) = qf(a) = qf(a, p + 1, n p 1) PrF qf(1 α) = 1 α β : (β ˆβ) (X X)(β ˆβ) (p + 1)Sɛ 2 qf(1 α) ˆβ Prβ (1 α) = 1 α 2. S 2 ɛ (n p 1)S 2 ɛ ˆβ (p + 1) (p + 1) R X X = R R χ 2 n p 1 QR X = QR R (X X) 1 = (R R) 1 = R 1 R 1 γ = Rβ, ˆγ = R ˆβ γ β = R 1 γ γ γ : γ ˆγ 2 (p + 1)Sɛ 2 qf(1 α) ˆγ Sɛ (p + 1)qf(1 α) p + 1 z = (z1, z2,..., zp+1) z = 1 σ R( ˆβ β) z E(z) = 0 V (z) = 1 σ Rσ2 (X X) 1 R 1 σ = RR 1 R 1 R = Ip+1 z Np+1(0, Ip+1) z1, z2,..., zp+1 N(0, 1) 23 # run0078.r # # = = source("run0075.r") p <- 1 # =1 n <- length() # alpha < # =1-alpha a <- func0075a(,p) # a$qr <- qr(a$x); a$r <- qr.r(a$qr) # QR a$ir <- solve(a$r) # R^-1 b0 <- func0075b(,alpha,a) # 24

7 th <- seq(0,2*pi,len=300) # 0..2*pi 300 ga <- b0$se*calcq1(n,p,alpha)*rbind(cos(th),sin(th)) # gamma be <- as.vector(b0$be) + a$ir %*% ga # beta plot(t(be),tpe="l") # points(t(b0$be)) # abline(v=b0$beconf[1,],lt=2,col=2) # beta0 abline(h=b0$beconf[2,],lt=2,col=2) # beta1 b1 <- func0075b(,alpha,a,calcq1) # calcq1 abline(v=b1$beconf[1,],lt=3,col=4) # beta0 abline(h=b1$beconf[2,],lt=3,col=4) # beta1 dev.cop2eps(file="run0078-b.eps") coef <- cbind(b0$be,b0$beconf,b1$beconf) colnames(coef) <- c("estimate","lo","up","jointlo","jointup") print(coef) > dat <- read.table("dat0001.tt") # (47 2 ) > <- dat/100; <- dat[,2] > source("run0078.r") Estimate Lo Up JointLo JointUp beta beta Prβ1 1 = 1 α Prβ0 0, β1 1 < 1 α Prβ0 0, β1 1 1 α β 4.2 v = w β (1 α) = [ˆv qt(1 α/2), ˆv + qt(1 α/2)] = Sɛ w (X X) 1 w beta qt(a) = qt(a,n-p-1) Prv (1 α) = 1 α w W Prw β w (1 α), w W 1 α beta0 run0078-b β = (β0, β1) α = 0.05 Prβ = 1 α βi, i = 0, 1 t- ( ) β0 β1 W Prβ0 0, β1 1 1 α Prβ0 0 = 1 α w (1 α) = [ˆv (p + 1)qf(1 α), ˆv + (p + 1)qf(1 α)] qf(a) = qf(a, p + 1, n p 1) # run0079.r # # = = p=, alpha source("run0075.r") func0079 <- function(,,a,alpha,output=t) b0 <- func0075b(,alpha,a) b1 <- func0075b(,alpha,a,calcq1) coef <- cbind(b0$be,b0$bese,b0$tval,b0$pval,b0$beconf,b1$beconf) colnames(coef) <- c("estimate","stderr","t-value","p-value", "Lo","Up","JointLo","JointUp") if(output) func0075c(,,a,b0); func0075c(,,a,b1,col=4,lt=3,add=t) cat("rss=",b0$rss,", RSQ=",b0$rsq,"\n") print(round(coef,3)) list(b0=b0,b1=b1,coef=coef) func0079b <- function(,,p,alpha,filename="run0079-") for(i in 0:p) cat("# =",i,"\n") a <- func0075a(,i) func0079(,,a,alpha) if(!is.null(filename)) dev.cop2eps(file=paste(filename,"s",i,".eps",sep="")) beta # = 1 RSS= , RSQ= beta beta # = 2 RSS= , RSQ= beta beta beta # = 3 RSS= , RSQ= beta beta beta beta run0079-s0 run0079-s1 > source("run0079.r") > dat <- read.table("dat0001.tt") # (47 2 ) > <- dat/100; <- dat[,2] > p <- 3; alpha < > func0079b(,,p,alpha) # = 0 RSS= , RSQ=

8 4.4 run0079-s2 α = % (Up,Lo) (JointLo,JointUp) 4.3 v = w β, w W run0079-s3 γ = Rβ, a = R 1 w v = w β = a γ (Schwartz) a ˆγ γ ˆv v = a (ˆγ γ) a ˆγ γ γ γ : γ ˆγ Sɛ (p + 1)qf(1 α) γ ˆv v a Sɛ (p + 1)qf(1 α) γ 1 α Pr ˆv v a Sɛ (p + 1)qf(1 α) 1 α a = w R 1 R 1 w = w (X X) 1 w w = v : ˆv v Sɛ w (X X) 1 w (p + 1)qf(1 α) # run0080.r # # run0074.r # run0075.r,run0079.r source("run0079.r") func0080a <- function(v) # check if v1 in [v2,v3] (v[1] >= v[2]) && (v[1] <= v[3]) func0080b <- function(be,00,d) coef <- cbind(be,d$coef) be0 <- appl(coef[,c("be","lo","up")],1,func0080a) be1 <- appl(coef[,c("be","jointlo","jointup")],1,func0080a) pred0 <- all(appl(cbind(00,d$b0$pred0conf),1,func0080a)) pred1 <- all(appl(cbind(00,d$b1$pred0conf),1,func0080a)) list(be0=be0,be1=be1,pred0=pred0,pred1=pred1) a <- func0075a(,p) # 00 <- a$x0 %*% be # 300 cat("# \n") d <- func0079(,sim,a,alpha) # abline(be,col=3,lt=4) # cat("# \n") esno1 <- unlist(func0080b(be,00,d)); print(esno1) dev.cop2eps(file="run0080-s1.eps") cat("# \n") esno <- matri(0,length(esno1),b) # rownames(esno) <- names(esno1) for(i in 1:b) d <- func0079(,sim[,i],a,alpha,output=f) esno[,i] <- unlist(func0080b(be,00,d)) print(esno[,1:5]) cat("# \n") print(appl(esno,1,sum)) > n <- 11 # > be <- c(0,1) # (beta0,beta1) > <- seq(0,1,len=n) # > filename <- NULL > source("run0074.r") # start simulation: Wed Oct 6 11:46: # end simulation: Wed Oct 6 11:46: # $be beta beta $se2 [1] $tval beta beta $pval beta e-01 beta e-05 # pval0 < 0.05 = 508 # pval1 < 0.05 = 8751 > alpha < > source("run0080.r") # RSS= , RSQ= beta beta # be0.beta0 be0.beta1 be1.beta0 be1.beta1 pred0 pred1 TRUE TRUE TRUE TRUE FALSE TRUE # [,2] [,3] [,4] [,5] be0.beta be0.beta be1.beta be1.beta pred pred # be0.beta0 be0.beta1 be1.beta0 be1.beta1 pred0 pred > sum(esno[1,] & esno[2,]) # joint0 [1] 9235 > sum(esno[3,] & esno[4,]) # joint1 [1] run0080-s1 β0 = 0, β1 = =95%9500 #β0 0 = 9498 #β1 1 = 9453 β0 β1 95% #β0 0, β1 1 = % #β0 0, β1 1 = % #β0 + β1 (), = 9500, #β0 + β1 (), = %95% 32

9 #β0 + β1 (), = % X2000 p = 1, 2, 3 t p- p, 95% 95% run0075.r calcq0 <- function(n,p,alpha) qt(1-alpha/2,n-p-1) # calcq1 <- function(n,p,alpha) sqrt((p+1)*qf(1-alpha,p+1,n-p-1)) # n = 30, α = 0.05 calcq0(n,p,alpha) calcq1(n,p,alpha) p = 0, 1,..., 10 =p= ) ( ) 33

.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 ))

.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 )) Copright (c) 004,005 Hidetoshi Shimodaira 1.. 3. 4. 004-10-01 16:15:07 shimo cat(" 1: "); c(mea(), mea()) cat(" : "); mmea