# Prvni priklad - LRM
# nacteni dat
load("cviceni5.RData")
X <- data5$pr1$X
Y <- data5$pr1$Y
d <- dim(X)

# vypocet
XX<- solve(t(X)%*%X)
XY<- t(X)%*%Y
beta <- XX%*%XY
Yhat <- X%*%beta
Se <- t(Y-Yhat)%*%(Y-Yhat)
s <- Se/(d[1]-d[2])

# Druhy priklad
# nacteni dat
load("cviceni5.RData")
x <- data5$pr2$X
Y <- data5$pr2$y
d <- length(x)

# vypocet
# model: beta_0 + beta_1 x
model <- lm(Y~1+x)
(prehled <- summary(model))
beta <- model$coefficients
Yhat <- model$fitted.values
s <- prehled$sigma^2
ID <- prehled$r.squared
F <- prehled$fstatistic["value"]
M <- prehled$coefficients
M[,4] # p-hodnoty testu pro jednotlive koeficienty
# vykresleni
x11(w=12,h=9)
plot(x,Y,type="p")
curve(predict(model,data.frame(x=x)),add=T)
#abline(model)

# model: beta_1 x
model <- lm(Y~0+x)
(prehled <- summary(model))
beta <- model$coefficients
Yhat <- model$fitted.values
s <- prehled$sigma^2
ID <- prehled$r.squared
F <- prehled$fstatistic["value"]
M <- prehled$coefficients
M[,4] # p-hodnoty testu pro jednotlive koeficienty
# vykresleni
x11(w=12,h=9)
plot(x,Y,type="p")
curve(predict(model,data.frame(x=x)),add=T,col=2)

# model: beta_0 + beta_1 x + beta_2 x^2
model <- lm(Y~1+x+I(x^2))
(prehled <- summary(model))
beta <- model$coefficients
Yhat <- model$fitted.values
s <- prehled$sigma^2
ID <- prehled$r.squared
F <- prehled$fstatistic["value"]
M <- prehled$coefficients
M[,4] # p-hodnoty testu pro jednotlive koeficienty
# vykresleni
x11(w=12,h=9)
plot(x,Y,type="p")
curve(predict(model,data.frame(x=x)),add=T,col=3)

# model: beta_1 x + beta_2 x^2
model <- lm(Y~0+x+I(x^2))
(prehled <- summary(model))
beta <- model$coefficients
Yhat <- model$fitted.values
s <- prehled$sigma^2
ID <- prehled$r.squared
F <- prehled$fstatistic["value"]
M <- prehled$coefficients
M[,4] # p-hodnoty testu pro jednotlive koeficienty
# vykresleni
x11(w=12,h=9)
plot(x,Y,type="p")
curve(predict(model,data.frame(x=x)),add=T,col=4)

# model: beta_0 + beta_1 x + beta_2 exp(x)
model <- lm(Y~1+x+I(exp(x)))
(prehled <- summary(model))
beta <- model$coefficients
Yhat <- model$fitted.values
s <- prehled$sigma^2
ID <- prehled$r.squared
F <- prehled$fstatistic["value"]
M <- prehled$coefficients
M[,4] # p-hodnoty testu pro jednotlive koeficienty
# vykresleni
x11(w=12,h=9)
plot(x,Y,type="p")
curve(predict(model,data.frame(x=x)),add=T,col=5)

# Treti priklad
# nacteni dat
load("cviceni5.RData")
x <- data5$pr3$X
Y <- data5$pr3$Y

# vypocet
# model: beta_0 + beta_1 x
model <- lm(Y~1+x)
(prehled <- summary(model))
# vykresleni
x11(w=12,h=9)
plot(x,Y,type="p")
curve(predict(model,data.frame(x=x)),add=T)

# model: beta_1 x
model <- lm(Y~0+x)
(prehled <- summary(model))
# vykresleni
curve(predict(model,data.frame(x=x)),col="red",add=T)

# teploty
rm(list=ls())
load("teploty.RData")
attach(teploty)
n1<-length(Y1)
n2<-length(Y2)
alpha <- 0.05
#x<-x-1977
mod1 <- lm(Y1~1+x)
mod2 <- lm(Y2~1+x)
sum1 <- summary(mod1)
sum2 <- summary(mod2)
b1 <- sum1$coefficients[2]
b2 <- sum2$coefficients[2]
M1<-model.matrix(mod1)
library(Matrix)
X<-as.matrix(bdiag(M1,M1))
XX <- t(X)%*%X
XXinv <- solve(XX)
# vypocet sigma - 1.zpusob
s1<-sum1$sigma
s2<-sum2$sigma
s<-((n1-2)*s1^2+(n2-2)*s2^2)/(n1+n2-4)
# vypocet sigma - 2.zpusob
mod <- lm(c(Y1,Y2)~X)
(sum <- summary(mod))
s<-sum$sigma
# Testovani rovnobeznosti
# t-statistika a kvantil
t <- abs((b1-b2)/(s*sqrt(2*XXinv[2,2])))
t
qt(1-alpha/2,n1+n2-4)
# Testovani shodnosti primek
beta1 <- sum1$coefficients[,1]
beta2 <- sum2$coefficients[,1]
W <- 2*solve(t(M1)%*%M1)
K1 <- t(beta1-beta2)%*%solve(W)%*%(beta1-beta2)
K2 <- s^2
(f <- K1/(2*K2))
(obor_prijeti <- c(qf(alpha/2,2,n1+n2-4),qf(1-alpha/2,2,n1+n2-4)))
# Testovani shodnosti rozptylu
(f <- s1^2/s2^2)
(obor_prijeti <- c(qf(alpha/2,n1-2,n2-2),qf(1-alpha/2,n1-2,n2-2)))

# vykresleni
x11(w=12,h=9)
plot(c(x[1],tail(x,1)),c(min(c(Y1,Y2)),max(c(Y1,Y2))),type="n",xlab="roky",ylab="teplota")
points(x,Y1,pch=1,col="red")
points(x,Y2,pch=20,col="green")
curve(predict(mod1,data.frame(x=x)),add=T,col="red")
# abline(mod1, col = "red")   # alternativa
curve(predict(mod2,data.frame(x=x)),add=T,col="green")
# abline(mod2, col = "green)  # alternativa
legend("bottomright",legend=c("Praha","V. Pavlovice"),col=c("red","green"),lty=c(1,1))
detach(teploty)
# -----------------------------------------------------------------