第1章　線形回帰（問題1～20）

3

N <- 100  # サンプル数をNとする

# 直線の係数をランダムに生成
a <- rnorm(1)
b <- rnorm(1)

# 直線の周りの点をランダムに生成
x <- rnorm(N)
y <- a*x + b + rnorm(N)

plot(y ~ x)  # 点のプロット

abline(lm(y ~ x), col = "red")  # あてはめ直線
abline(h = 0)  # x軸
abline(v = 0)  # y軸

# N組の点の重心を原点にうつす
x <- x - ## 空欄(1) ##
y <- y - ## 空欄(2) ##

abline(lm(y ~ x), col = "blue")  # あてはめ直線
abline(h = 0)  # x軸
abline(v = 0)  # y軸

legend("topleft", c("中心化前", "中心化済み"), lty = 1, col = c("red", "blue"))

12

N <- 100
x <- rnorm(N)
y <- rnorm(N)
x.bar <- mean(x)
y.bar <- mean(y)
beta.0 <- sum(y.bar * sum(x^2) - x.bar * sum(x*y)) / sum((x - x.bar) ^ 2)
beta.1 <- sum((x - x.bar) * (y - y.bar)) / sum((x - x.bar) ^ 2)
RSS <- sum((y - beta.0 - beta.1*x) ^ 2)
RSE <- sqrt(RSS / (N - 1 - 1))
B.0 <- sum(x^2) / N / sum((x - x.bar)^2)
B.1 <- 1 / sum((x - x.bar) ^ 2)
se.0 <- RSE * sqrt(B.0)
se.1 <- RSE * sqrt(B.1)
t.0 <- beta.0 / se.0
t.1 <- beta.1 / se.1
p.0 <- 2 * (1 - pt(abs(t.0), N - 2))  # p値（その値より外側にある確率）
p.1 <- 2 * (1 - pt(abs(t.1), N - 2))  # p値（その値より外側にある確率）

beta.0
se.0
t.0
p.0

beta.1
se.1
t.1
p.1

lm(y ~ x)
summary(lm(y ~ x))

13

N <- 100
r <- 1000
T <- NULL
for (i in 1:r) {
  x <- rnorm(N)
  y <- rnorm(N)
  x.bar <- mean(x)
  y.bar <- mean(y)
  fit <- lm(y ~ x)
  beta <- fit$coefficients
  RSS <- sum((y - fit$fitted.values) ^ 2)
  RSE <- sqrt(RSS / (N - 1 - 1))
  B.1 <- 1 / sum((x - x.bar) ^ 2)
  se.1 <- RSE * sqrt(B.1)
  T <- c(T, beta[2] / se.1)
}

hist(T, breaks = sqrt(r), probability = TRUE, xlab = "tの値", ylab = "確率密度",
     main = "tの値のヒストグラムと理論値（赤）")
curve(dt(x, N - 2), -3, 3, type = "l", col = "red", add = TRUE)

15(d)

R2 <- function(x, y) {
  y.hat <- lm(y ~ x)$fitted.values
  y.bar <- mean(y)
  RSS <- sum((y - y.hat) ^ 2)
  TSS <- sum((y - y.bar) ^ 2)
  return(1 - RSS / TSS)
}

N <- 100
m <- 2
x <- matrix(rnorm(m * N), ncol = m)
y <- rnorm(N)
R2(x, y)

16

library(MASS)
X <- as.matrix(Boston)
p <- ncol(X)
T <- NULL
for (j in 1:p)
  T <- c(T, ## 空欄##)
T

18

# データを生成
N <- 100
p <- 1
X <- matrix(rnorm(N * p), ncol = p)
X <- cbind(rep(1, N), X)
beta <- rnorm(p + 1)
epsilon <- rnorm(N)
y <- X %*% beta + epsilon

# 関数f(x)を定義。Uは t(X) %*% X の逆行列
U <- solve(t(X) %*% X)
beta.hat <- U %*% t(X) %*% y
RSS <- sum((y - X %*% beta.hat) ^ 2)
RSE <- sqrt(RSS / (N - p - 1))
alpha <- 0.05

f <- function(x) {
  x <- cbind(1, x)
  range <- qt(df = N - p - 1, 1 - alpha / 2) * RSE * sqrt(x %*% U %*% t(x))
  return(list(lower = x %*% beta.hat - range, upper = x %*% beta.hat + range))
}

# グラフで信頼区間を表示
x.seq <- seq(-10, 10, 0.1)
lower.seq <- NULL
for (x in x.seq)
  lower.seq <- c(lower.seq, f(x)$lower)
upper.seq <- NULL
for (x in x.seq)
  upper.seq <- c(upper.seq, f(x)$upper)
x.lim <- c(min(x.seq), max(x.seq))
y.lim <- c(min(lower.seq), max(upper.seq))
plot(x.seq, lower.seq, col = "blue", xlim = x.lim, ylim = y.lim,
     xlab = "x", ylab = "y", type = "l")
par(new = TRUE)
plot(x.seq, upper.seq, col = "red", xlim = x.lim, ylim = y.lim,
     xlab = "", ylab = "", type = "l", axes = FALSE)
abline(beta.hat[1], beta.hat[2])