# 1

# a

# Daten einlesen
data <- read.table("E:\\wasser.txt", header = TRUE)

# Designmatrix
X <- cbind(1, data$Luftdruck)

# Zielvariablen
Y <- data$Siedepunkt

# Parameter n und m
n <- length(Y)
m <- length(X[1,])

# MKQ-Schätzer für beta
# Lösung der Normalengleichung Satz 3.2.1
beta <- solve(t(X) %*% X, t(X) %*% Y)

# Schätzer für sigma^2 (Kapitel 3.2.2)
sigma2 <- 1 / (n - m) * t(Y - X %*% beta) %*% (Y - X %*% beta)

# Plot
plot(data$Luftdruck, data$Siedepunkt)
plot(function(x) beta[1] + beta[2] * x, min(data$Luftdruck), max(data$Luftdruck), add = TRUE)

# b

beta0 <- c(150, 2)

# Testgröße (Satz 3.2.8)
T <- t(beta - beta0) %*% t(X) %*% X %*% (beta - beta0) / (m * sigma2)
print(T)

# Quantil
print(qf(0.95, m, n - m))

# Testentscheidung
if (T > qf(0.95, m, n - m)) {
	print("H0 wird abgelehnt")
} else {
	print("H0 wird nicht abgelehnt")
}

# c

# siehe Satz 3.2.11

# obere Schranke
upper <- function(x) {
	yupper <- rep(0, length(x))
	for (i in 1:length(x)) {
		yupper[i] = beta[1] + beta[2] * x[i] + sqrt(m * qf(.95, m, n-m) * sigma2 * c(1, x[i]) %*% solve(t(X) %*% X, c(1, x[i])))
	}
	return(yupper)
}

# untere Schranke
lower <- function(x) {
	ylower <- rep(0, length(x))
	for (i in 1:length(x)) {
		ylower[i] = beta[1] + beta[2] * x[i] - sqrt(m * qf(.95, m, n-m) * sigma2 * c(1, x[i]) %*% solve(t(X) %*% X, c(1, x[i])))
	}
	return(ylower)
}

# Plot
plot(upper, min(data$Luftdruck), max(data$Luftdruck), add = TRUE)
plot(lower, min(data$Luftdruck), max(data$Luftdruck), add = TRUE)