A=matrix(c(1,2,3,4),ncol=2)
svd(A)

svd(A)$u%*%diag(svd(A)$d)%*%t(svd(A)$v)

#zgadza się - to jest równe A



library(faraway)


summary(meatspec)

help(meatspec)

model1=lm(fat~.,meatspec)
summary(model1)

step(model1)

model2=lm(formula = fat ~ V1 + V2 + V4 + V5 + V7 + V8 + V11 + V12 + 
     V13 + V14 + V15 + V17 + V19 + V20 + V22 + V24 + V25 + V26 + 
     V28 + V29 + V30 + V32 + V34 + V36 + V37 + V39 + V40 + V41 + 
     V42 + V45 + V46 + V47 + V48 + V50 + V51 + V52 + V54 + V55 + 
     V56 + V59 + V60 + V61 + V63 + V64 + V65 + V67 + V68 + V69 + 
     V71 + V73 + V74 + V78 + V79 + V80 + V81 + V84 + V85 + V87 + 
     V88 + V92 + V94 + V98 + V99, data = meatspec)

summary(model2)
plot(model2)

X=model.matrix(model1)
head(X)
X=X[,2:101]
data=X
for(i in 1:100) {
  X[,i]=X[,i]-mean(X[,i])    #uśredniamy
  X[,i]=X[,i]/var(X[,i])}




d=svd(X)$d
plot(svd(X)$d)
plot(svd(X)$d[1:10],type="p")

V=svd(X)$v

pca<-prcomp(data,center=TRUE,scale.=TRUE)
plot(pca,type="l")
summary(pca)


help(prcomp)


# X%*%V[,j] - nowa zmienna objaśniająca (j-ty kierunek główny)

g1=X%*%V[,1]
g2=X%*%V[,2]
g3=X%*%V[,3]
g4=X%*%V[,4]
g5=X%*%V[,5]

kg=data.frame(g1=g1,g2=g2,g3=g3,g4=g4,g5=g5,fat=meatspec$fat)


model3=lm(fat~g1+g2+g3+g4+g5-1,kg)
summary(model3)

#sprawdzamy krzyżowo



data=meatspec


r1=0
r2=0
r3=0

for (j in (1:length(data)))
{
  training=data[-j,]
  test=data[j,]
  
  modelt1=lm(fat~.,training)
  modelt2=lm(formula = fat ~ V1 + V2 + V4 + V5 + V7 + V8 + V11 + V12 + 
               V13 + V14 + V15 + V17 + V19 + V20 + V22 + V24 + V25 + V26 + 
               V28 + V29 + V30 + V32 + V34 + V36 + V37 + V39 + V40 + V41 + 
               V42 + V45 + V46 + V47 + V48 + V50 + V51 + V52 + V54 + V55 + 
               V56 + V59 + V60 + V61 + V63 + V64 + V65 + V67 + V68 + V69 + 
               V71 + V73 + V74 + V78 + V79 + V80 + V81 + V84 + V85 + V87 + 
               V88 + V92 + V94 + V98 + V99, data = training)
  
  Xt=model.matrix(modelt1)
Xt=Xt[,2:101]
for(i in 1:100) 
{
  Xt[,i]=Xt[,i]-mean(Xt[,i])    #uśredniamy
  Xt[,i]=Xt[,i]/var(Xt[,i])
}
  
  
Vt=svd(Xt)$v

  g1t=Xt%*%Vt[,1]
  g2t=Xt%*%Vt[,2]
  g3t=Xt%*%Vt[,3]
  g4t=Xt%*%Vt[,4]
  g5t=Xt%*%Vt[,5]
  g6t=Xt%*%Vt[,6]
  g7t=Xt%*%Vt[,7]
  kgt=data.frame(g1t=g1t,g2t=g2t,g3t=g3t,g4t=g4t,g5t=g5t,g6t=g6t,g7t=g7t,fatt=training$fat)
  
  modelt3=lm(fatt~g1t+g2t+g3t+g4t+g5t+g6t+g7t-1,kgt)
  x=as.numeric(test[1:100])
  
  vtest=data.frame(g1t=sum(x*Vt[,1]),g2t=sum(x*Vt[,2]),g3t=sum(x*Vt[,3]),g4t=sum(x*Vt[,4]),g5t=sum(x*Vt[,5]),g6t=sum(x*Vt[,6]),g7t=sum(x*Vt[,7]))
  # obliczenie nowych zmiennych, kierunków głównych, potrzebnych do prognozy w modelu 3


  r1=r1+abs(predict(modelt1,test) - test$fat    )
  r2=r2+abs(predict(modelt2,test) - test$fat    )
  r3=r3+abs(predict(modelt3,vtest) -test$fat    )
  
  
}


rmse1=r1/length(data)
rmse2=r2/length(data)

rmse3=r3/length(data)


rmse1
rmse2
rmse3