ZADANIE
Copyright by B. Kamiñski i M. Zawisza
 
Zbudowaæ 2 modele: na pe³nym zbiorze danych i na takim,
gdzie zmienne s¹ mocno skorelowane z zmienn¹ zale¿n¹.
Porównaæ jakoœæ modeli na zbiorze walidacyjnym. Na podstawie danych {http://kdd.ics.uci.edu/databases/tic/tic.html}.
 
library(car); 
set.seed(1)
 
names.file <- "http://kdd.ics.uci.edu/databases/tic/dictionary.txt"
var.names <- readLines(names.file)
 
var.names <- var.names[4:89]
var.names <- strsplit(var.names,split=" ")
var.names <- sapply(var.names, function(x) { return (x[2]) } )
 
data.flie <- "http://kdd.ics.uci.edu/databases/tic/ticdata2000.txt"
data.set <- read.delim(data.flie, header = F)
names(data.set) <- var.names
data.set <- data.set[,c(-2,-3)]
 
values <- rep(rep(c("training", "validation"), c(3,2)),
              len = nrow(data.set))
assign.split <- factor((sample(values,nrow(data.set))))
data.set.split <- split(data.set, assign.split)
 
min.level <- 0.000001
restricted.vars <- ncol(data.set)
for (i in 1:(ncol(data.set) - 1)) {
  if (cor.test(data.set.split$training$CARAVAN,
               data.set.split$training[,i])$p.value < min.level) {
    restricted.vars <- c(restricted.vars,i)
  }
}
 
 
 
name <- c("pełny", "selekcja")
model1 <- glm(CARAVAN ~ ., data=data.set.split$training,
              family="binomial")
model2 <- glm(CARAVAN ~ .,
              data=data.set.split$training[,restricted.vars],
              family="binomial")
 
 
valid1 <- predict(model1,
                  newdata=data.set.split$validation)
 
valid2 <- predict(model2,
                  newdata=data.set.split$validation)
 
res_valid1<-(data.set.split$validation$CARAVAN-valid1)
res_valid2<-(data.set.split$validation$CARAVAN-valid2)
 
 
 
 
 
 
 
 
ZADANIE
Copyright by B. Kamiñski i M. Zawisza
 
Jedn¹ z metod bezpoœredniej oceny b³êdu konkuruj¹cych klasyfikatorów
jest podzia³ zbioru pierwotnego na trzy zbiory:
  ucz¹cy, walidacyjny i testowy. Zbiór ucz¹cy s³u¿y do oszacowania konkuruj¹cych
modeli. Zbiór walidacyjny s³u¿y do wyboru jednego z wielu klasyfikatorów, tj.
tego, który ma najmniejszy b³¹d na niezale¿nym zbiorze walidacyjnym.
Niemniej jednak ocena takiego b³êdu na zbiorze walidacyjnym
mo¿e byæ zbyt optymistyczna, tj. za niska. St¹d potrzebujemy
jeszcze zbiór testowy do nieobci¹¿onej oceny b³êdu.
 
Dokonaæ eksperymentu symulacyjnyego, który poka¿e potrzebê
u¿ywania zbioru testowego na przyk³adzie predykcji binarnej.
Wygeneruj 256 obserwacji zmiennej binarnej, w której prawdopodobieñstwo
wyst¹pienia jedynki wynosi 50%. Nastêpnie stwórz 1024 niezale¿nych
klasfikatorów losowych, tj. takich których prawdopodobieñstwo
poprawnej klasyfikacji zmiennej binarnej wynosi 50%.
Podziel wygenerowany zbiór 256 obserwacji na równoliczne zbiory
walidacyjny i testowy. Przedstaw rozk³ad jakoœci klasyfikatorów
na zbiorze walidacyjnym oraz oceñ najlepszy z nich na zbiorze testowym.
 
set.seed(1)
observations <- 128
 
classifiers.number <- 1024
 
validation.data <- rep(c(0,1), observations / 2)
test.data <- rep(c(0,1), observations / 2)
 
 
 
classifier.accuracy <- function(i, ref.data) {
  prediction <- rbinom(length(ref.data), 1, 0.5)
  mean(prediction == ref.data)
}
 
accuracy.validation <- sapply(1:classifiers.number,
                              classifier.accuracy, ref.data = validation.data)
best.classifier <- which.max(accuracy.validation)
 
cat("Rozklad trafnosci na zbiorze walidacyjnym:\n")
print(summary(accuracy.validation))
cat("Trafnosc najlepszego klasyfikatora na zbiorze testowym:",
    classifier.accuracy(best.classifier, test.data),"\n")
 
 
 
 
 
 
 
dane <- read.csv(file.choose(),header = F)
 
names(dane)<-c("crime","zoned","industry","river", "nox","rooms","age","distances",
               "highways","tax","pt.ratio","blacks","lstat","medvalue")
 
names(dane)
 
attach(dane)
 
# Model pe³ny
mnk<-lm(medvalue~., dane)
 
summary(mnk)
 
 
model.AIC<-step(mnk, direction="both", k=2)
 
 
summary(model.AIC)
 
 
model.BIC<-step(lm(medvalue~., dane), direction="both", k=log(nrow(dane)))
 
 
 
library(leaps)
all<-regsubsets(medvalue~.,data=dane,nvmax=13)
summary(all)
 
 
par(mfrow=c(1,3))
plot(all, scale="adjr2", col=rainbow(256))
plot(all, scale="Cp", col=rainbow(256))
plot(all, scale="bic", col=rainbow(256))
 
 
 
 
 
#CV
 
model<-glm(medvalue~., data=dane)
 
library(boot)
dane.cv<-cv.glm(dane, model, K=2)
 
 
dane.cv5<-cv.glm(dane, model, K=5)
dane.cv5$delta
 
dane.cv10<-cv.glm(dane, model, K=10)
dane.cv10$delta
 
dane.cv_loo<-cv.glm(dane, model)
dane.cv_loo$delta
 
 
 
 
 
 
 
 
set.seed(1)
# autor: dr B. Kamiñski
cv.split <- function(dane, K) {
  cvf <- vector("list",K)
  # przemieszanie wierszy
  dane <- dane[sample(nrow(dane)),]
  # nadanie wierszom etykiet przynale¿noœci do k'tej podpróby
  fold.n <- rep(1:K,length.out=nrow(dane))
  for (i in 1:K) {
    cvf[[i]] <- list(train=dane[fold.n!=i,], valid=dane[fold.n==i,])
  }
  return(cvf)
}
 
K=5
zbior=cv.split(dane,K)
 
length(zbior)
 
 
names(zbior[[1]])
 
nrow(zbior[[1]]$valid)
nrow(zbior[[1]]$train)
 
nrow(zbior[[2]]$valid)
nrow(zbior[[3]]$valid)
nrow(zbior[[4]]$valid)
nrow(zbior[[5]]$valid)
 
mean(zbior[[1]]$train)
 
library(lasso2)
shrinkageFactor=seq(0.01,1,0.01)
 
pred.err=c()
for(b in shrinkageFactor){
  skr=0
  for(i in 1:K){
    lasso<-l1ce(medvalue~.,data=zbior[[i]]$train,bound=b)
    pred=predict(lasso,newdata=zbior[[i]]$valid)
    skr=skr+sum((pred-zbior[[i]]$valid$medvalue)^2)
  }
  pred.err=c(pred.err,skr)
}
 
length(pred.err)
length(shrinkageFactor)
par(mfrow=c(1,1))
plot(shrinkageFactor,pred.err)
 
plot(shrinkageFactor,pred.err,ylim=c(10000,15000))
grid()
 
plot(shrinkageFactor,pred.err,ylim=c(12000,13000))
grid()
 
plot(shrinkageFactor,pred.err,ylim=c(12400,12800))
grid()
 
shrinkageFactor[which.min(pred.err)]
 
 
 
lasso<-l1ce(medvalue~.,data=dane,bound=shrinkageFactor,standardize=T)
 
# ...na razie dla jednego parametru...
plot(shrinkageFactor,coef(lasso)[,2],
     "l",ylim=c(-18,4),lwd=2,col=2,xlim=c(0,1.2),
     xlab="Wzglêdne ograniczenie",
     ylab="Oszacowanie parametrów",
     main="Regresja LASSO")
grid()
text(1.1,coef(lasso)[length(shrinkageFactor),2],names(dane)[1],col=2)
 
# i teraz dla pozosta³ych...
for(i in 3:14){
  lines(shrinkageFactor,coef(lasso)[,i],col=i,lwd=2)
  text(1.1,coef(lasso)[length(shrinkageFactor),i],names(dane)[i-1],col=i)
}
 
 
# Pytanie: jaka jest optymalna
#	wartoϾ 'bound'?
abline(v=shrinkageFactor[which.min(pred.err)],lty=2,col="blue")
 
# na zbiorze walidacyjnym
lasso[which.min(pred.err)]
 
lasso[length(pred.err)]
 
 
# BOOTSTRAP
 
 
 
parametry<-function(X,d,y,i)
{Xt<-t(X)
 beta=(1/(Xt[,d]%*%X[d,]))%*%(Xt[,d]*y[d,])
 beta<-as.vector(beta)
 return(beta[i])
}
 
X1<-dane[-14]
X1<-as.matrix(X1)
y1<-dane[14]
y1<-as.vector(y1)
for (m in 1:13)
{print(boot(X1,parametry,1000, y=y1,i=m))
}
 
detach(dane, unload=TRUE)
 
 
# autor: dr B. Kamiñski
bootdane <- function(dane, K) {
  cvf <- vector("list",K)
  # przemieszanie wierszy
  dane <- dane[order(runif(nrow(dane))),]
  # nadanie wierszom etykiet przynale¿noœci do k'tej podpróby
  # fold.n <- rep(1:K,length.out=nrow(dane))
  for (i in 1:K) {
    drawing=sample(nrow(dane),replace=T)
    cvf[[i]] <- list(train=dane[drawing,], valid=dane[-unique(drawing),])
  }
  return(cvf)
}
 
zbior=bootdane(dane,100)
 
nrow(zbior[[1]]$train)
nrow(zbior[[1]]$valid)
 
nrow(zbior[[2]]$train)
nrow(zbior[[100]]$valid)
 
nrow(zbior[[3]]$train)
nrow(zbior[[3]]$valid)
 
K=50
 
pred.err=c()
for(b in shrinkageFactor){
  skr=0
  for(i in 1:K){
    lasso<-l1ce(medvalue~.,data=zbior[[i]]$train,bound=b)
    pred=predict(lasso,newdata=zbior[[i]]$valid)
    skr=skr+sum((pred-zbior[[i]]$valid$medvalue)^2)/length(pred)
  }
  skr=skr/K
  pred.err=c(pred.err,skr)
}
 
length(pred.err)
length(shrinkageFactor)
plot(shrinkageFactor,pred.err)
 
shrinkageFactor[which.min(pred.err)]
 
 
 
lasso<-l1ce(medvalue~.,data=dane,bound=shrinkageFactor,standardize=T)
 
# ...na razie dla jednego parametru...
plot(shrinkageFactor,coef(lasso)[,2],
     "l",ylim=c(-18,4),lwd=2,col=2,xlim=c(0,1.2),
     xlab="Wzglêdne ograniczenie",
     ylab="Oszacowanie parametrów",
     main="Regresja LASSO")
grid()
text(1.1,coef(lasso)[length(shrinkageFactor),2],names(dane)[1],col=2)
 
# i teraz dla pozosta³ych...
for(i in 3:14){
  lines(shrinkageFactor,coef(lasso)[,i],col=i,lwd=2)
  text(1.1,coef(lasso)[length(shrinkageFactor),i],names(dane)[i-1],col=i)
}
 
 
# Pytanie: jaka jest optymalna
#	wartoϾ 'bound'?
abline(v=shrinkageFactor[which.min(pred.err)],lty=2,col="blue")
 
# na zbiorze walidacyjnym
lasso[which.min(pred.err)]
 
lasso[length(pred.err)]