B <- read.table("M.txt") # nacteni matice M
colnames(B)<-cbind('d','lambda','i','k','s')
B<-cbind(B,B$d^2)
B<-cbind(B,B$d^3)
B<-cbind(B,B$lambda^2)
B<-cbind(B,B$lambda^3)
B<-cbind(B,B$i^2)
B<-cbind(B,B$i^3)
colnames(B)<-cbind('d','lambda','i','k','s','d2','d3','lambda2','lambda3','i2','i3')

C<-read.table("T.txt") # nacteni matice T
colnames(C)<-cbind('d','lambda','i','k','s')
C<-cbind(C,C$d^2)
C<-cbind(C,C$d^3)
C<-cbind(C,C$lambda^2)
C<-cbind(C,C$lambda^3)
C<-cbind(C,C$i^2)
C<-cbind(C,C$i^3)
colnames(C)<-cbind('d','lambda','i','k','s','d2','d3','lambda2','lambda3','i2','i3')

# histogramy
k<-B$k
hist(k)
s<-B$s
hist(s)

######################
# jednoduchy lin. model
######################

m<-lm(k~d+lambda+i,data=B)
summary(m)
plot(m)
RSS<-sum(resid(m)^2)
TSS<-sum((B$k-mean(B$k))^2)
# kod = koeficient determinace
kod<-1-RSS/TSS 

######################
# po transformaci (Box-Cox)
######################

m<-lm(k^-0.5~d+lambda+i,data=B)
summary(m)
plot(m)

RSS<-sum(resid(m)^2)
SST<-sum((B$k^-0.5-mean(B$k^-0.5))^2)
kod_bc<-1-RSS/SST

# predikce
pred<-predict(m,newdata=C)
p<-pred^-2
RMS_bc<-sqrt(sum((p-C$k)^2)/109)

# velky plot
sim<-C$k
plot(sim,p)
lines(c(0,max(C$k)),c(0,max(C$k)),type='l')

# omezeny plot 
D<-cbind(p,sim)
D<-D[D[,2]<3,]
sim<-D[,2]
p<-D[,1]
plot(sim,p)
lines(c(0,3),c(0,3),type='l')

######################
# polynomialni model pro k
######################

mp<-lm(k^-0.5~I(d)+I(d2)+I(d3)+
         I(lambda)+I(lambda2)+I(lambda3)+
         I(i)+I(i2)+I(i3)+
         I(d*lambda)+I(d*lambda2)+I(d2*lambda)+  
         I(d*i)+I(d*i2)+I(d2*i)+
         I(i*lambda)+I(i*lambda2)+I(i2*lambda)+
         I(d*i*lambda),data=B)

summary(mp)
# vyrazeni promennych navic
mps<-step(mp)
summary(mps)
# plot(mps)

# predikce
pred<-predict(mps,newdata=C)
p<-pred^-2
RMS_mps<-sqrt(sum((p-C$k)^2)/109)

RSS<-sum((p-C$k)^2)
TSS<-sum((C$k-mean(C$k))^2)
kod_mps<-1-RSS/TSS

# velky plot
sim<-C$k
plot(sim,p,xlab="k_sim",ylab="k_pred")
lines(c(0,max(C$k)),c(0,max(C$k)),type='l')

# zmenseni scalu plotu
D<-cbind(p,C$k)
D<-D[D[,2]<3,]
sim<-D[,2]
p<-D[,1]
plot(sim,p,xlab="k_sim",ylab="k_pred")
lines(c(0,3),c(0,3),type='l')

######################
# glm pro k
######################

mg<-glm(k~I(d)+I(d2)+I(d3)+
          I(lambda)+I(lambda2)+I(lambda3)+
          I(i)+I(i2)+I(i3)+
          I(d*lambda)+I(d*lambda2)+I(d2*lambda)+  
          I(d*i)+I(d*i2)+I(d2*i)+
          I(i*lambda)+I(i*lambda2)+I(i2*lambda)+
          I(d*i*lambda),family=Gamma(link=log),data=B)
summary(mg)

mgs<-step(mg)
# protoze nebyla vyrazena i2 a i3, musim do modelu vratit i
mgs<-update(mgs, . ~ . + I(i), data=B)
summary(mgs)

# plot(mgs)

# predikce
pred<-predict(mgs,newdata=C,type="response")
RMS_mgs<-sqrt(sum((pred-C$k)^2)/109)

RSS<-sum((pred-C$k)^2)
TSS<-sum((C$k-mean(C$k))^2)
kod_mgs<-1-RSS/TSS

# velky plot
sim<-C$k
p<-pred
plot(sim,p,xlab="k_sim",ylab="k_pred")
lines(c(0,max(C$k)),c(0,max(C$k)),type='l')

# zmenseni scalu plotu
D<-cbind(p,C$k)
D<-D[D[,2]<3,]
sim<-D[,2]
p<-D[,1]
plot(sim,p,xlab="k_sim",ylab="k_pred")
lines(c(0,3),c(0,3),type='l')

######################
# UNS
######################

A<-read.table("UNSk.txt")
colnames(A)<-cbind('sim','p')

RSS<-sum((A[,1]-A[,2])^2)
TSS<-sum((A$sim-mean(A$sim))^2)
kod_UNS<-1-RSS/TSS

RMS_UNS<-sqrt(sum((A[,1]-A[,2])^2)/109)

# velky plot
plot(A,xlab="k_sim",ylab="k_pred")
lines(c(0,max(A)),c(0,max(A)),type='l')

# omezeny plot
plot(A[A[,1]<3,],xlab="k_sim",ylab="k_pred")
lines(c(0,3),c(0,3),type='l')