# # # Written by: # -- # John L. Weatherwax 2009-04-21 # # email: wax@alum.mit.edu # # Please send comments and especially bug reports to the # above email address. # #----- # Figure 1.1: # a=c(5,5.6,1,4,-5) a[1] b=a[2:4] b %% 3 d %/% 2.4 e = 3/d log(d * e) t(d) e[d==5] a[-3] is.vector(d) # Ex 1.4: # x = runif(5) x order(x) x[order(x)] # is a sorted list rank(x) # gives the spot in the sorted list of each element of x x = c(1,2,3) rep(x,times=5) rep(x,length.out=5) rep(x,each=5) # Ex 1.5: # n = 10 1:n-1 seq(1,n-1,by=1) n=1 1:n-1 # Figure 1.2: # x1 = matrix(1:20,nrow=5) x2 = matrix(1:20,nrow=5,byrow=T) x3 = t(x2) b = x3 %*% x2 b[-2,] rbind(x1,x2) cbind(x1,x2) apply(x1,1,sum) apply(x1,2,sum) as.matrix(1:10) t(as.matrix(1:10)) # Figure 1.3: # state = c("tas","tas","sa","sa","wa") statef = factor(state) levels(statef) incomes = c(60,59,40,42,23) tapply(incomes,statef,mean) statef=factor(state,levels=c("tas","sa","wa","yo")) table(statef) # Figure 1.4: # li = list(num=1:5,y="color",a=T) a = matrix(c(6,2,0,2,6,0,0,0,36),nrow=3) res = eigen(a,symmetric=T) names(res) res$vectors diag(res$values) res$vec %*% diag(res$values) %*% t(res$vec) x = list( a = 1:10, beta = exp(-3:3), logic = c(T,F,F,T) ) lapply( x, mean ) sapply( x, mean ) # Figure 1.5: # v1 = sample( 1:12, 30, rep=T ) v2 = sample( LETTERS[1:10], 30, rep=T ) v3 = runif(30) v4 = rnorm(30) xx = data.frame(v1,v2,v3,v4) # Ex. 1.6: # x = rnorm(20) y = 3*x + 5 + rnorm(20,sd=0.3) reslm = lm(y~x) summary(reslm) b = matrix(1:9,ncol=3) var(b) sd(b)^2 x = rnorm(25) t.test(x) attach(faithful) cor.test(faithful[,1],faithful[,2]) ks.test(jitter(faithful[,1]),pnorm) shapiro.test(faithful[,2]) wilcox.test(faithful[,1]) Nit = c(0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,6,6,6) AOB = c(4.26,4.15,4.68,6.08,5.87,6.92,6.87,6.25,6.84,6.34,6.56,6.52,7.39,7.38,7.74,7.76,8.14,7.22) AOBm = tapply(AOB,Nit,mean) Nitm = tapply(Nit,Nit,mean) plot(Nit,AOB,xlim=c(0,6),ylim=c(min(AOB),max(AOB)),pch=19) library(splines) fitAOB = lm(AOBm~ns(Nitm,df=2)) xmin=min(Nit); xmax=max(Nit) lines(seq(xmin,xmax,0.5), predict(fitAOB,data.frame(Nitm=seq(xmin,xmax,0.5)))) fitAOB2 = loess(AOBm~Nitm,span=1.25) lines(seq(xmin,xmax,0.5), predict(fitAOB2,data.frame(Nitm=seq(xmin,xmax,0.5)))) arima(diff(EuStockMarkets[,1]),order=c(0,0,5)) acf(ldeaths,plot=F) acf(ldeaths,plot=T) acf(ldeaths,plot=T,type="partial") y = c(4.313, 4.513, 5.489, 4.265, 3.641, 5.106, 8.006, 5.087) x = seq(-3,3,le=5) y = 2 + 4 * x + rnorm(5) lm(y~x) fit = lm(y~x) Rdata = fit$residuals nBoot = 2000 B=array(0,dim=c(nBoot,2)) for(i in 1:nBoot){ ystar = y + sample(Rdata,replace=T) Bfit = lm(ystar~x) B[i,] = Bfit$coefficients } x=rnorm(1) for( t in 2:10^3) x = c(x,0.09*x[t-1]+rnorm(1)) plot(x,type="l",xlab="time",ylab="x",lwd=2,lty=2,col="steelblue",ylim=range(cumsum(x))) lines(cumsum(x),lwd=2,col="orange3") par(mar=c(2,2,2,2)) x=matrix(0,ncol=100,nrow=10^4) for( t in 2:10^4 ) x[t,]=x[t-1,]+rnorm(100)*10^(-2) plot(seq(0,1,le=10^4),x[,1],ty="n",ylim=range(x),xlab="",ylab="") polygon(c(1:10^4,10^4:1)/10^4,c(apply(x,1,max),rev(apply(x,1,min))),col="gold",bor=F) polygon(c(1:10^4,10^4:1)/10^4,c(apply(x,1,quantile,0.95),rev(apply(x,1,quantile,0.05))),col="brown",bor=F) # Wait! You can get all of this stuff for free and don't have to type in in yourself! # install.packages("mcsm") #The downloaded packages are in #/tmp/RtmpOC1IJe/downloaded_packages library(mcsm) demo(Chapter.1)