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')

######################
# jednoduchy lrm - spatny
######################

m<-lm(s~d+lambda+i,data=B)
summary(m)
# plot(m)

p<-predict(m,newdata=C)
RMS_simple<-sqrt(sum((p-C$s)^2)/109)
sim<-C$s
plot(sim,p)

######################
# transformace - taky spatny
######################

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

# predikce
pred<-predict(m,newdata=C)
p<-pred^2
sim<-C$s
RMS_trans<-sqrt(sum((p-sim)^2)/109)
plot(sim,p)
lines(c(0,30000),c(0,30000))

######################
# polynomialni model pro s - pouzitelny
######################
mps<-lm(s^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(mps)

# predikce
pred<-predict(mps,newdata=C)
p<-pred^2
sim<-C$s
RMS_mps<-sqrt(sum((p-sim)^2)/109)
RSS<-sum((p-sim)^2)
TSS<-sum((sim-mean(sim))^2)
kod_mps<-1-RSS/TSS
plot(sim,p,xlab='s_sim',ylab='s_pred')
lines(c(0,30000),c(0,30000),type='l')

######################
# glm - se zaokrouhlenymi hodnotami
######################
gm<-glm(round(s)~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,family=poisson)
summary(gm)
gms<-step(gm)
# plot(gms)
summary(gms)

# predikce
pred<-predict(gms,newdata=C,type="response")
RMS_gms<-sqrt(sum((pred-C$s)^2)/109)
p<-pred
sim<-C$s
plot(sim,p,xlab='s_sim',ylab='s_pred')
lines(c(0,max(p)),c(0,max(p)))

RSS<-sum((p-sim)^2)
TSS<-sum((sim-mean(sim))^2)
kod_gms<-1-RSS/TSS

######################
# UNS
######################
A<-read.table("UNSs.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='s_sim',ylab='s_pred')
lines(c(0,max(A)),c(0,max(A)),type='l')