# Exercise 1

euler <- function() {
	return(exp(1))
}
euler()

numpow <- function(x) x^2
numpow(10)

numpow2 <- function(x, p = 2) x^p
numpow2(3, 3)

adjr2 <- function(obj) summary(obj)$adj.r.squared
adjr2(out)

chgr2 <- function(obj1, obj2) summary(obj2)$r.squared - summary(obj1)$r.squared

out1 <- lm(ACT ~ SATQ, data = sat.act)
out2 <- lm(ACT ~ SATQ + SATV, data = sat.act)
chgr2(out1, out2)

describe(sat.act$SATQ)[,c("skew", "kurtosis")]

skewkur <- function(x) describe(x)[,c("skew", "kurtosis")]

skewkur(attitude[,1])

# Exercise 9.2

x <- 2
y <- "i"

if(is.numeric(y)) {
	print("This is a number.")
} else {
	print("This is NOT a number.")
}



x <- c(1:5, NA)

if(!is.vector(x)) {
	print("This is NOT a vector.")
} else {
	if(any(is.na(x))) {
		print(which(is.na(x)))
	} else {
		print("This is a vector that does not have any NA.")
	}
}

# any: check whether any elements is TRUE. If so, return TRUE. If all are FALSE, return FALSE.
# all: check whether all elements are TRUE. If so, return TRUE. If one of them is FALSE, return FALSE.
# which: return the position(s) of the elements that are TRUE.

raise <- function(x, pow) {
	if(pow < 0) {
		return(NA)
	} else {
		return(x^pow)
	}
}



attitude2 <- data.frame(attitude, f = factor(rep(1:2, each = 15)))


alpha <- function(dat) {
	checknumeric <- apply(dat, 2, is.numeric)
	if(all(checknumeric)) {
		p <- ncol(dat)
		term1 <- p / (p - 1)
		CM <- cov(dat)
		numerator <- sum(diag(CM))
		sumscore <- apply(dat, 1, sum)
		denominator <- var(sumscore)
		term2 <- 1 - numerator/denominator
		alpha <- term1 * term2
		return(alpha)
	} else {
		return(NA)
	}
}
alpha(attitude)
alpha(attitude2)

lengthtransform <- function(inches, unit = "cm") {
	result <- NA
	if(unit == "cm") {
		result <- inches * 2.54
	} else if (unit == "m") {
		result <- inches * 0.0254
	} else if (unit == "yards") {
		result <- inches * 0.0277778
	} else {
		stop("The specified units are not available.")
	}
	return(result)
}

lengthtransform(1, "cm")
lengthtransform(1, "m")
lengthtransform(1, "yards")
lengthtransform(1, "km")



# Exercise 9.3

library(boot)
FUN1 <- function(dat, i) cor(dat[i, 1], dat[i, 2])

bootout1 <- boot(attitude, FUN1, stype = "i", R = 1000)
bootci1 <- boot.ci(bootout1, conf.level = 0.95, type = c("perc", "bca"))

n.male.progun <- 60
n.female.progun <- 45
n.male.antigun <- 40
n.female.antigun <- 55

sex <- rep(c("male", "female", "male", "female"), c(n.male.progun, n.female.progun, n.male.antigun, n.female.antigun))
gun <- rep(c("pro", "pro", "anti", "anti"), c(n.male.progun, n.female.progun, n.male.antigun, n.female.antigun))
datgun <- data.frame(sex, gun)

FUN2 <- function(dat, i) {
	tab <- table(dat[i,])
	pfemale <- tab[1, 2] / (tab[1, 1] + tab[1, 2])
	pmale <- tab[2, 2] / (tab[2, 1] + tab[2, 2])
	pfemale - pmale
}

bootout2 <- boot(datgun, FUN2, stype = "i", R = 1000)
bootci2 <- boot.ci(bootout2, conf.level = 0.95, type = c("perc", "bca"))


# Assignment 9

fun1 <- function(d, h, piece) (pi * ((d/2)^2) * h) / piece

fun1(30, 10, 8)

fun2 <- function(v) v[!duplicated(v)]

fun2(c("a", "x", "e", "a", "s", "s"))

fun3 <- function(v) abs(v)

fun3(c(-1, 1, 3, -3))

library(psych)
library(car)

out <- aov(SATQ ~ education*gender, data = sat.act)
out2 <- Anova(out, type = 3)

fun4 <- function(obj) {
	ss <- obj[,1]
	ss <- ss[-1]
	sst <- sum(ss)
	names(ss) <- rownames(obj)[-1]
	result <- ss/sst
	result[-length(result)]
}
fun4(out2)



library(boot)
fun5 <- function(dat, i) {
	out <- aov(SATQ ~ education*gender, data = dat[i,])
	obj <- Anova(out, type = 3)
	ss <- obj[,1]
	ss <- ss[-1]
	sst <- sum(ss)
	names(ss) <- rownames(obj)[-1]
	result <- ss/sst
	result[-length(result)]
}

bootout5 <- boot(sat.act, fun5, stype = "i", R = 1000)
boot.ci(bootout5, conf.level = 0.95, type = c("perc", "bca"), index = 1)
boot.ci(bootout5, conf.level = 0.95, type = c("perc", "bca"), index = 2)
boot.ci(bootout5, conf.level = 0.95, type = c("perc", "bca"), index = 3)

fun6 <- function(dat, i) {
	out <- aggregate(SATQ ~ gender, data = dat[i,], FUN = median)
	out[1, 2] - out[2, 2]
}

bootout6 <- boot(sat.act, fun6, stype = "i", R = 1000)
boot.ci(bootout6, conf.level = 0.95, type = c("perc", "bca"))