Stat 6304: Computational Statistics

Code used in class:

mat <- structure(c(1,3,7,9,2,8), dim=c(3,2))
mat
attributes(mat)

abc <- ts(c(18,23,76,98,198,203,208,237,239,243,278), start = 2000, frequency = 1)
abc
attribute(abc)

y <- c(1.2, 3.4, 5.6, 12.34)
storage.mode(y)
storage.mode(y) <- "integer"
y
storage.mode(y) <- "double"
as.integer(y)
storage.mode(as.integer(y))

NA <- 12324
pi <- 3.4454354
rm(pi)

x <- c(4,8,15)
y <- c("Texas", "New York", "Ohio")

library(tseries)
library(TSA)

search()
find("acf")
acf <- 7843265432
find("acf")
acf
rm(acf)

mat1 <- 23432
mat2 <- 12121
objects(pattern = "m", pos = 1)
ls()

x <- scan()
scan("e:/stat6304_10/aaa.txt")
y <- scan("e:/stat6304_10/aaa.txt")
xx <- scan("http://faculty.smu.edu/ngh/sta6104/ch5eg.txt")

read.fwf("e:/stat6304_10/temp.txt", widths = c(2,6,3,2,2))

site <- "http://www1.appstate.edu/~arnholta/PASWR/CD/data/Baberuth.txt"
Baberuth <- read.table(file=url(site), header = T)
Baberuth

num <- 1:8
num <- 1:100
num <- seq(1, 8, length = 8)
num <- seq(1, 8, length = 10)
num <- seq(1, 8, by = .4)
num <- seq(1, 8.1, by = .4)
xx <- seq(-3, 3, length = 1000)
yy <- exp(-xx*xx/2)/sqrt(2*pi)
plot(xx, yy)

rep(1, 100)
rep(c(1,2,4), 100)
rep(1:10, 9)
rep(1:5, 5:1)

cl <- list(dept="STAT", number="6304", title="Computational Statistics",
enroll=8, classif=c(1,1,2,3,2,4) )
cl$title
cl$classif[4]
mode(cl)
storage.mode(cl)
unlist(cl)
aa <- unlist(cl)
mode(aa)

mat1 <- matrix(1:15, nrow = 3, ncol = 5)
mat1 <- matrix(1:15, nrow = 3, ncol = 5, byrow = T)
mat1 <- matrix(1:15, nrow = 3)
mat1 <- matrix(1:18, nrow = 3)
mat1 <- matrix(1:19, nrow = 3)
mat1 <- matrix(1:9, nrow = 3, ncol = 6)

history()
 

x <- seq(-3,3, length = 1000)
y <- exp(-x*x/2)/sqrt(2*pi)

plot(x,y)
dev.off()

motif()
x11()
window()

pdf.graph("temp.pdf" , horizontal = T)
dev.off()

pdf("temp.pdf")
plot(Cars93$MPG.city, Cars93$Weight)

pdf("temp.pdf", height = 11, width = 8.5)
plot(Cars93$MPG.city, Cars93$Weight, pch = 16)

pdf("c:/downloads/abcdef.pdf", paper = "letter")
plot(Cars93$MPG.city, Cars93$Weight, pch = 16)

pdf("c:/downloads/abcdef.pdf", height = 11, width = 8.5)
plot(Cars93$MPG.city, Cars93$Weight, pch = 16)


dev.set(dev.prev())
plot(x,y, type = "l")
dev.set(dev.prev())
plot(car.miles, car.gals)

plot(Cars93$MPG.city, Cars93$Weight)

dev.cur()
dev.list()
dev.next()
dev.prev()
dev.off()
dev.set()


plot(Cars93$MPG.city, type = "p")
plot(Cars93$MPG.city, type = "l")
plot(Cars93$MPG.city, type = "b")
plot(Cars93$MPG.city, type = "o")
plot(Cars93$MPG.city, type = "h")
plot(Cars93$MPG.city, type = "s")
plot(Cars93$MPG.city, type = "n")
plot(Cars93$MPG.city, type = "c")

plot(1:6,1:6,type ="b", cex = 1:6)

plot(1:18,1:18, pch = 1:18)
plot(1:18,1:18, pch = 1:18, cex = 3)
plot(32:150,32:150, pch = 32:150)

pie(rep(1,15), col=1:15)}
plot(1:50, 1:50, col = 1:50, pch = 16, cex = 3)
points(1:50, 1:50, col = 1:50, pch = 15, cex = 3)
colors()
plot(1:50, 1:50, col = colors()[1:50], pch = 16, cex = 3)
points(1:50, 50:1, col = colors()[51:100], pch = 15, cex = 3)

plot(Cars93$MPG.city, Cars93$Weight, type = "p", log="x")

plot(Cars93$MPG.city, type = "n")
points(Cars93$MPG.city, pch = 5, col = 6, type = "b")
points(Cars93$MPG.city+10, pch = 0, col = 9, type = "p")
points(Cars93$MPG.city-10, pch = 16, col = 18, type = "p")
points(Cars93$MPG.city, type ="l", lwd = 5)


x <- seq(-3, 3, length = 100)
plot(x, exp(-x*x/2)/(sqrt(2*pi)), type="b", pch=8, xlab = "x axis label", ylab = "y axis label", main = "Main Title", sub = "Subtitle")

x <- 1:10
y <- x*x
plot(x,y)
identify(x,y, n =2)
z <- 100 - y
points(x,z, col = 5, pch = "A")

plot(car.miles, car.gals, xlab=" ", ylab=" ")
title(xlab="Miles per trip", ylab="Gallons per trip", main="Mileage Data", sub="Gallons vs. Miles")

plot(car.miles, car.gals, axes=F)
axis(1)
axis(2)
axis(3)
axis(4)

plot(Cars93$MPG.city, Cars93$Weight, xlab=" ", ylab=" ")
title(xlab="City MPG", ylab="Weight (kg)", main="Car Data", sub="Weight vs. MPG")

plot(Cars93$MPG.city, Cars93$Weight, axes=F)
axis(1)
axis(2)
axis(3)
axis(4)
box()

help(box)
box(lty = 5, col = 'red')

plot(Cars93$MPG.city, Cars93$Weight, xlim = c(0, 400), ylim = c(0, 30))
plot(Cars93$MPG.city, Cars93$Weight, xlim = c(0, 400), ylim = c(0, 30), axes = F)
plot(Cars93$MPG.city, Cars93$Weight, axes=F)
axis(1, at = c(20, 30, 40), labels = c("A", "B", "C"))

mtext("text", side = 3, line = 0, outer = F, at = 40)
mtext("lower", side = 3, line = -2, outer = F, at = 40)
mtext("upper", side = 3, line = 2, outer = F, at = 40)

abline(lm(Cars93$Weight~Cars93$MPG.city))
text(34, 2100, "Regression Line" , srt = -46, crt = -46)
arrows(20, 2000, 25, 2400, col = 3, lwd = 4)
segments(40, 3000, 35, 2500, col = 4)

abline(lm(car.gals~car.miles))
text(300, 17, "Regression Line" , srt = 35, crt = 35)
arrows(250, 20, 300, 15, col = 9)
segments(150, 10, 250, 20, col = 4)

locator(n = 2, type = "l")
locator(n = 4, type = "l", col = 2)
locator(n = 6, type = "b", col = 5)

text(locator(1), "outlier")

x <- c(200, 250, 200, 250)
y <- c(10, 10, 21, 21)
polygon(x,y)

x <- c(15, 25, 20, 25)
y <- c(2000, 2000, 3100, 3100)
polygon(x,y)
polygon(x,y, col = 6)

pretty(Cars93$MPG.city)
pretty(Cars93$MPG.city, 10)

plot(Cars93$MPG.city)
rug(Cars93$MPG.city)

split.screen(c(3,1))
split.screen(c(1,3), screen = 1)
split.screen(c(1,2), screen = 3)
screen(4)
plot(x)
screen(5)
plot(x*x)
screen(6)
plot(1:10)
screen(2)
hist(Cars93$Length)
screen(7)
hist(Cars93$Weight)
screen(8)
hist(Cars93$Price)
close.screen(all = TRUE)

par(mfrow = c(3,3))
plot(x)
plot(x*x)
hist(x)
par(mfrow = c(3,1))
par(mfg = c(2,1))
plot(x)
par(mfrow = c(3,2))
par(mfg = c(3,1))
plot(x)
plot(x*x)

pdf(width = 3, height = 3)
plot(xx,yy)
mtext("Outer = F")
mtext("Outer = T", outer = T, line = 0)
mtext("Outer = T -1", outer = T, line = -1)
mtext("Outer = T -2", outer = T, line = -2)
dev.off()

a <- seq(0.1, 3, length = 1000)
plot(a, gamma(a), type = "l", xlab = "", ylab ="")
title(xlab = expression(alpha == sum(frac(x[i], n), i==1, n)))
title(ylab = expression(Gamma(alpha) == integral(x^{alpha-1} * e^{-x}*dx, 0, infinity)) , line = 1.5)

demo(plotmath)

x <- seq(0,1, length = 20)
y <- x*x
xx <- x
yy <- x
plot(x,y, pch=16)
points(xx,yy,pch=3)
legend(locator(1), c(expression(y == x^{2}), "y = x"), pch = c(16, 3))

x <- rnorm(100, 0,1) 
y <- rexp(100, 1) 
boxplot(x, y, names=c("Normal(0,1)", "Exp(1)")) 
 

hist(Cars93$MPG.city, probability = T) 
points(density(Cars93$MPG.city), type = "l") 
dimnames(Cars93)[[1]] <- Cars93$Make  
stars(Cars93[, 1:7]) 
stars(Cars93[1:20, 1:7])

library(TeachingDemos) 
faces2(matrix( runif(18*10), nrow=10)) 

pie(rep(1, 24), col = rainbow(24), radius = 0.9) 
pie.sales <- c(0.12, 0.3, 0.26, 0.16, 0.04, 0.12) 
names(pie.sales) <- c("Blueberry", "Cherry",  "Apple", "Boston Cream", "Other", "Vanilla Cream") 
pie(pie.sales)  
pie(pie.sales,   col = c("purple", "violetred1", "green3", "cornsilk", "cyan", "white"))

cars <- c(1, 3, 6, 4, 9)  
trucks <- c(2, 5, 4, 5, 12)  
suvs <- c(4,4,6,6,16)  
autos_data <- data.frame(cars, trucks, suvs)

barplot(as.matrix(autos_data), main="Autos", ylab= "Total", beside=TRUE, col=rainbow(5)) 
abline(h = 0)

legend("topleft", c("Mon","Tue","Wed","Thu","Fri"), cex=1.6, bty="n", fill=rainbow(5));

barplot(t(autos_data), main="Autos", ylab="Total", col=heat.colors(3), space=0.1, cex.axis=0.8, las=1, 
names.arg=c("Mon","Tue","Wed","Thu","Fri"), cex=0.8)

barplot(t(autos_data), main="Autos", ylab="Total",  
col=heat.colors(3), space=0.1, cex.axis=0.8, las=1, 
names.arg=c("Mon","Tue","Wed","Thu","Fri"), cex=0.8)  
legend(1, 30, names(autos_data), cex=1.2, fill=heat.colors(3)) 
 

Trellis in R:

library(lattice)  
bwplot(Cars93$Passengers~Cars93$MPG.city)  
 

3D scatterplot in R:

library(scatterplot3d)  
scatterplot3d(Cars93$MPG.city, Cars93$Weight, Cars93$Horsepower, pch = 15)

library(rgl)  
plot3d(Cars93$MPG.city, Cars93$Weight, Cars93$Horsepower, col="blue", size=6)

library(Rcmdr)  
scatter3d(Cars93$MPG.city, Cars93$Weight, Cars93$Horsepower)

Programming in R:

ource("E:\\STAT6304_09fall\\source.R")

history()
savehistory("")
loadhistory("")

sink("E:\\STAT6304_09fall\\abcde.txt")
x <- rnorm(100)
x
sink()

1 + 1
3 - 5
6*7
16/3
2^3
16%/%3
19%%3

m1 <- matrix(1:6, nrow=2, ncol=3)
m2 <- matrix(7:12, nrow=3, ncol=2)
m1%*%m2
m1%o%m2


v1 <- c(1,2,3)

v1 > 10 & v1 > 15
v1 > 10 && v1 > 15
v1 < 2 && sgdfsg

if (all(v1 > 2)) {v3 <- 2*v1}
if (any(v1 > 2)) {v3 <- 2*v1}

- 5 * 8 / 4 + 6

round(4.5)
round(5.5)
round(4.555, digits = 2)
round(4.5545, digits = 3)

0.5 + -2:4
round(0.5 + -2:4)

zapsmall(4.5545, digits = 2)
round(0.15, 1)


x1 <- c(52938, 1928, 12, 2, 0.0053)
x2 <- c(52938, 1928, 12, 2, 0.005)
y <- c(52938, 1928, 12, 2, 0)

log(max(x1) / min(x1), 10)

all(zapsmall(x1, 6) - y == 0)
all(zapsmall(x1, 7) - y == 0)


x2 <- pi * 100^(-1:3)
print(x2 / 1000, digits=4)
zapsmall(x2 / 1000, digits=4)



x <- 4.5467468437687745654
round(x, digits = 4)
format(x, digits = 15)

 


signif(x, digits = 2)
signif(x, digits = 4)
signif(x, digits = 10)

xx <- c(5,4,2,6,9,10,3,8, NA)
sort(xx)
rev(sort(xx))
rev(list(a=1:2, b=5:6))

matx <- matrix(c(4,5,6,2,4,1,7,3,8,7,3,4,9,6,1,3,12,7,19,8,4,5,2,10), ncol=3, nrow=8)
sort(matx)
ordx <- order(matx[,2])
matx[ordx,]

order(matx[,2], matx[,1])

mat <- matrix(1:12, ncol = 3)
scale(mat)

xmat <- matrix(c(1,2,3,2,5,8,3,8,13), ncol = 3)
y <- eigen(mat)
y$vector %*% diag(y$values) %*% t(y$vector)

xm <- matrix(c(5,1,1,3), ncol = 2)
chol(xm)
t(chol(xm)) %*% chol(xm)
crossprod(chol(xm))

state.name
abbreviate(state.name)
abbreviate(state.name, 2)
abbreviate(state.name, 5, dot = T)

nchar(state.name)
paste(state.name, "USA")
substring(state.name, 2)

pmatch("Ala", state.name)
pmatch("Alas", state.name)
pmatch(c("Alas", "Alas"), state.name, duplicates.ok = T)
pmatch(c("Alas", "Alas"), state.name, duplicates.ok = F)


str(Cars93)
head(Cars93)
tail(Cars93)
 

x <- 1:10
replace(x, 3, 1.5)
append(x, c(3.4, 5.6))
append(x, c(3.4, 5.6), 6)
replace(x, c(3, 8), 0)

xx <- c(1, 1, 2.3, 3.5, 6, 3.5)
unique(xx)
duplicated(xx)
yy <- c(1, 2, 3.5, 6, 4, 4, 7)
yy[-duplicated(yy)]
yy[!duplicated(yy)]