# Prvni priklad - overovani
# nacteni dat
load("minuta.RData")

# histogram
h <- hist(x,probability=TRUE, breaks="Scott",col="linen",
        xlab="odhad munity",
        main="Histogram, Normal Curve",ylim=c(0,0.08))
xx <- seq(from=min(h$breaks),to=max(h$breaks), length.out=50)
mx <- mean(x)
sx <- sd(x)
lines(xx,dnorm(xx,mean=mx,sd=sx),col=2)

# qqplot
qqnorm(x)
qqline(x,col=2)

# vyberova distribucni funkce
n <- length(x)
z <- sort((x-mx)/sx)
emp_F <- (1:n)/n
teor_F <- pnorm(z)
plot(z,emp_F,type='l',ylab="distribucni fce")
lines(z,teor_F,type="l",col="red")
legend("topleft",legend=c("Empiricka","Teoreticka"),col=c("black","red"),lty=c(1,1))

# KS test
ks.test(z,"pnorm")

# Shapiro-Wilk
shapiro.test(x)

# Test dobre shody
p <- pnorm(c(h$breaks[c(2,3)],Inf),mean=mx,sd=sx) - pnorm(c(-Inf,h$breaks[c(2,3)]),mean=mx,sd=sx)
n*p
# !!! Pozor, data nesplnuji podminku !!!
chisq.test(h$counts,p=p)


# Druhy priklad
rm(list=ls())
load("voda.RData")
attach(voda)
mod <- lm(Y~x)
res <- mod$residuals
n <- length(res)
# graficky
plot(res[-n],res[-1],xlab="e_1,...,e_n-1",ylab="e_2,...,e_n",main="residual plot")

# odhad theta
theta1 <- as.numeric((t(res[-n])%*%res[-1])/(t(res[-n])%*%res[-n]))
# Durbin-Watson
#install.packages("lmtest")
library(lmtest)
ts <- dwtest(Y~x)
D <- ts$statistic
theta2 <- 1-D/2

# asymptoticky test
alpha <- 0.05
abs(sqrt(n)*theta1)
qnorm(1-alpha/2)  # Zamitame hypotezu H: theta = 0

# Durbin-Watson
ts$p.value # Taky zamita

# Odstraneni
Ynew <- Y[-1]-theta1*Y[-n]
xnew <- x[-1]-theta1*x[-n]
modnew <- lm(Ynew~xnew)
resnew <- modnew$residuals
# vykresleni novych residui
plot(resnew[-n+1],resnew[-1],xlab="e_1,...,e_n-1",ylab="e_2,...,e_n",main="residual plot")
# testovani pomoci DW
dwtest(Ynew~xnew)
# vykresleni
plot(xnew,Ynew)
curve(predict(modnew,data.frame(xnew=x)),add=T,col=2)

# podle vety
emod <- lm(res[-1]~res[-n])
s <- var(emod$residuals)
ind <- sapply (1:n, function (i) {sapply (1:n, function (j) {abs(i-j)})})
W <- theta1^ind
X <- cbind(rep(1,times=n),x)
Winv <- solve(W)
beta <- solve(t(X)%*%Winv%*%X)%*%t(X)%*%Winv%*%Y
res <- Y-X%*%beta
theta <- as.numeric((t(res[-n])%*%res[-1])/(t(res[-n])%*%res[-n]))

# asymptoticky test
alpha <- 0.05
abs(sqrt(n)*theta)
qnorm(1-alpha/2)  # Zamitame hypotezu H: theta = 0
plot(x,Y)
lines(x,X%*%beta)
modpuv <- lm(Y~x)
curve(predict(modpuv,data.frame(x=x)),add=T,col=2)

detach(voda)
# ---------------------------------------------------------------------------