# przykład 6.9

### Random forest
load("zip.dat")

zip.ucz<-zip.train
zip.test<-zip.test
rm(zip.train)
rm(zip.test)

set.seed(71)

lzm<-round((ncol(zip.ucz)-1)/2)

ks<-seq(1,128,by=2)

iter<-1

rf.blad.zip.test<-matrix(nrow=length(ks),ncol=iter)
#rf.blad.zip.ucz<-rf.blad.zip.test

for (j in 1:iter)
{
for (i in seq(along=ks))
{
lz<-ks[i] # liczba zmiennych
print(lz)

zip.rf<-randomForest(klasa~.,zip.ucz,ntree=50,mtry=lz)

#rf.blad.zip.ucz[i,j]<-1-sum(zip.rf$predicted==zip.ucz$klasa)/nrow(zip.ucz)

#zip.rf$err.rate[,1] # OOB dla kolejnych drzew od 1 do 50

Z.test<-predict(zip.rf,zip.test,type ="class")

rf.blad.zip.test[i,j]<-1-sum(Z.test==zip.test$klasa)/nrow(zip.test)

}
}

blad.zip<-rbind(zip.blad.test,rf.blad.zip.test)

# rys. 6.8a

plot(ks,blad.zip,type="n",xlab="Liczba zmiennych",ylab="Błąd predykcji",main="Zbiór Zip")
lines(ks,blad.zip,type="l",lty=1)
points(ks,blad.zip,pch=19)

#################
# Regresja
data(Boston)
bos.ucz<-Boston

y<-Boston$medv

lzm<-ncol(Boston)-1

ks<-seq(1,lzm,by=1)

iter<-1

rf.blad.bos<-matrix(nrow=length(ks),ncol=iter)
mse<-matrix(nrow=length(ks),ncol=iter)

for (j in 1:iter)
{
for (i in seq(along=ks))
{
lz<-ks[i] # liczba zmiennych
print(lz)

bos.rf<-randomForest(medv~ .,data=Boston,ntree=100,mtry=lz)

ypred<-predict(bos.rf,Boston)

rf.blad.bos[i,j]<-1-sum((y-ypred)^2)/(sum((y-mean(y))^2))
mse[i,j]<-sum((y-ypred)^2)/(nrow(Boston)-(lz+1))
}
}

############################
data(Boston)

m<-nrow(Boston)

lzm<-ncol(Boston)-1

ks<-seq(1,lzm,by=1)

iter<-100

rf.blad.bos<-matrix(nrow=length(ks),ncol=iter)
mse<-matrix(nrow=length(ks),ncol=iter)

for (j in 1:iter)
{
test<-sample(1:m,round(m/3),replace=FALSE)

bos.ucz<-Boston[-test,]
bos.test<-Boston[test,]

y<-bos.test$medv

for (i in seq(along=ks))
{
lz<-ks[i] # liczba zmiennych
print(lz)

bos.rf<-randomForest(medv~ .,data=bos.ucz,ntree=100,mtry=lz)

ypred<-predict(bos.rf,bos.test)

rf.blad.bos[i,j]<-1-sum((y-ypred)^2)/(sum((y-mean(y))^2))
mse[i,j]<-sum((y-ypred)^2)/(nrow(Boston)-(lz+1))
}
}

mse.100.bos<-numeric(length=lzm)
for (i in 1:lzm) mse.100.bos[i]<-mean(mse[i,])

# rys. 6.8b
plot(ks,mse.100.bos,type="n",xlab="Liczba zmiennych",ylab="MSE",main="Zbiór Boston")
lines(ks,mse.100.bos,type="l",lty=1)
points(ks,mse.100.bos,pch=19)