# Problem 1.
library(MASS)

hitters <- read.csv("hitters.csv")
hitters2 <- hitters[,-1]
hitters2 <- hitters2[c(19, 1:18, 20)]
hitters.names <- names(hitters2[-1])
n <- nrow(hitters2)
# 1)
y <- as.matrix(hitters2$Salary, ncol=1)
X <- as.matrix(hitters2[-1])
dim(X)
# [1] 263 19
x.star <- numeric(0) # keeps track of variables in M_k
# bic <- n * log(rss/n) + log(n) * (k+1)

ones <- matrix(rep(1, n), ncol=1)
k <- ncol(X)
# beta <- solve(t(X)%*%X)%*%t(X)%*%y
# e <- y - X%*%beta
# rss <- sum(e*e)
# n * log(rss/n) + log(n) * (k+1)

selected <- numeric(0)
bics <- numeric(19)
if (exists("X_start")) rm(X_star)

for (k in 1:19) {
bic <- 9999999
for (j in 1:19) {
if (length(selected)==0) {
   X.star <- cbind(ones)
} else {
   if (j %in% selected) {
   next
   } else {
   X.star <- cbind(ones, X[, selected])
   }
}
X.star <- cbind(X.star, X[,j])
beta <- solve(t(X.star)%*%X.star)%*%t(X.star)%*%y
e <- y - X.star%*%beta
rss <- sum(e*e)
if (n * log(rss/n) + log(n) * (k+1) < bic) {
   bic <- n * log(rss/n) + log(n) * (k+1)
   j.star <- j
}
}
bics[k] <- bic
print(bic)
selected <- c(selected, j.star)
}

names(bics) <- hitters.names[selected]
bics
# CRBI    Hits PutOuts Division AtBat Walks CWalks CRuns CAtBat Assists
# 3122.922 3084.842 3078.141 3071.959 3069.697 3065.851 3068.996 3066.386 3068.325 3070.552
# League    Runs Errors HmRun CHits    RBI NewLeague Years CHmRun
# 3074.908 3079.895 3084.991 3090.124 3095.550 3100.951 3106.415 3111.904 3117.464

plot(bics, type="l",
main="BIC vs Number of Variables",
xlab="Number of Variables",

1. Let us try independently a random event A until the first succes...