# This R environment comes with all of CRAN preinstalled, as well as many other helpful packages # The environment is defined by the kaggle/rstats docker image: https://github.com/kaggle/docker-rstats # For example, here's several helpful packages to load in library(ggplot2) # Data visualization library(readr) # CSV file I/O, e.g. the read_csv function # Input data files are available in the "../input/" directory. # For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory system("ls ../input") # Any results you write to the current directory are saved as output. library(data.table) library(tcltk2) library(randomForest) library(xgboost) library(Matrix) library(devtools) train <- fread('../input/train.csv', select = c( 'Producto_ID',"Agencia_ID",'Cliente_ID', 'Demanda_uni_equil',"Canal_ID","Ruta_SAK"), colClasses=c(Producto_ID="numeric", Agencia_ID = "numeric", Cliente_ID="numeric",Demanda_uni_equil="numeric",Canal_ID = "numeric ",Ruta_SAK = "numeric")) colClasses=c(Producto_ID="numeric",Agencia_ID = "numeric",Cliente_ID="numeric",Demanda_uni_equil="numeric" , Canal_ID = "numeric",Ruta_SAK = "numeric") setkey(train, Producto_ID, Agencia_ID, Cliente_ID) #calculate the overall median median <- train[, median(log(Demanda_uni_equil))] #calculate the product overall mean; call it M2 mean_Prod <- train[, mean(log(Demanda_uni_equil+1)), by = Producto_ID] setnames(mean_Prod,"V1","M2") #calculate the agent and product mean; call it M3 mean_Agent_Prod <- train[, mean(log(Demanda_uni_equil+1)),by = .(Producto_ID,Agencia_ID,Cliente_ID)] setnames(mean_Agent_Prod,"V1","M3") setkey(train, Producto_ID, Ruta_SAK) mean_Prod_Ruta <- train[, mean(log(Demanda_uni_equil+1)),by = .(Producto_ID,Ruta_SAK)] setnames(mean_Prod_Ruta,"V1","M4") max_Prod <- train[, max(log(Demanda_uni_equil+1)),by = .(Producto_ID)] setnames(max_Prod,"V1","M5") q_Prod <- train[, quantile(log(Demanda_uni_equil+1),.85),by = .(Producto_ID)] setnames(q_Prod,"V1","M6") setkey(train, Cliente_ID) max_Cliente <- train[, max(log(Demanda_uni_equil+1)),by = .(Cliente_ID)] setnames(max_Cliente,"V1","M7") q_Cliente <- train[, quantile(log(Demanda_uni_equil+1),.85),by = .(Cliente_ID)] setnames(q_Cliente,"V1","M8") setkey(train, Agencia_ID) max_Agen <- train[, max(log(Demanda_uni_equil+1)),by = .(Agencia_ID)] setnames(max_Agen,"V1","M9") q_Agen <- train[, quantile(log(Demanda_uni_equil+1),.85),by = .(Agencia_ID)] setnames(q_Agen,"V1","M10") gc() # I have been using samples to run my models. People with access to a lot of computing power should be able to use the whole thing train1 <- train[sample(nrow(train), 500000),] setkey(train1, Producto_ID, Agencia_ID, Cliente_ID) train2 <- merge(train1, mean_Agent_Prod, all.x = TRUE) train2$M2 <- merge(train1, mean_Prod, all.x = TRUE)$M2 setkey(train1, Producto_ID, Ruta_SAK) train2$M4 <- merge(train1, mean_Prod_Ruta, all.x = TRUE)$M4 setkey(train1, Producto_ID) train2$M5 <- merge(train1,max_Prod, all.x = TRUE)$M5 train2$M6 <- merge(train1,q_Prod,all.x = TRUE)$M6 setkey(train1, Cliente_ID) train2$M7 <- merge(train1,max_Cliente,all.x = TRUE)$M7 train2$M8 <- merge(train1,q_Cliente,all.x = TRUE)$M8 setkey(train1, Agencia_ID) train2$M9 <- merge(train1,max_Agen,all.x = TRUE)$M9 train2$M10 <- merge(train1,q_Agen,all.x = TRUE)$M10 train2$Canal_ID <- as.factor(train2$Canal_ID) train2 <- as.data.frame(train2) trainFea <- train2 for (rows in split(1:nrow(trainFea), ceiling((1:nrow(trainFea))/10000))) { ifelse(trainFea[rows,"Demanda_uni_equil"] < 0,.001,trainFea[rows,"Demanda_uni_equil"]) } trainFea$target <- log(trainFea$Demanda_uni_equil + 1) test <- fread('../input/test.csv', select = c("id", 'Producto_ID',"Agencia_ID",'Cliente_ID',"Canal_ID","Ruta_SAK"), colClasses=c(Producto_ID="numeric", Agencia_ID = "numeric", Cliente_ID="numeric",Canal_ID = "numeric",Ruta_SAK = "numeric")) test$Canal_ID = as.factor(test$Canal_ID) setkey(test, Producto_ID, Agencia_ID, Cliente_ID) test1 <- merge(test, mean_Agent_Prod,all.x = TRUE) test1$M2 <- merge(test, mean_Prod, all.x = TRUE)$M2 setkey(test, Producto_ID, Ruta_SAK) test1$M4 <- merge(test, mean_Prod_Ruta, all.x = TRUE)$M4 setkey(test, Producto_ID) rm(train, train1) gc() test1$M5 <- merge(test,max_Prod, all.x = TRUE)$M5 test1$M6 <- merge(test,q_Prod,all.x = TRUE)$M6 setkey(test, Cliente_ID) test1$M7 <- merge(test,max_Cliente,all.x = TRUE)$M7 test1$M8 <- merge(test,q_Cliente,all.x = TRUE)$M8 setkey(test, Agencia_ID) test1$M9 <- merge(test,max_Agen,all.x = TRUE)$M9 test1$M10 <- merge(test,q_Agen,all.x = TRUE)$M10 # this is done very crudely test1$M4 <- ifelse(is.na(test1$M4),2,test1$M4) test1$M5 <- ifelse(is.na(test1$M5),2,test1$M5) test1$M6 <- ifelse(is.na(test1$M6),2,test1$M6) test1$M7 <- ifelse(is.na(test1$M7),2,test1$M7) test1$M8 <- ifelse(is.na(test1$M8),2,test1$M8) test1$M2 <- ifelse(is.na(test1$M2),2,test1$M2) #Split the trianing set up into those with and without that one really good covariate test1G <- test1[!is.na(test1$M3),] test1B <- test1[is.na(test1$M3),] for(i in 1:1){ set.seed(i) tr <- trainFea h <- sample(nrow(tr), 400000) tr1 <-tr[-h,] feature.names <- names(tr)[-c(1,2,3,5,6,16,17)] gc() } dtrain <- xgb.DMatrix(data.matrix(tr1[,feature.names]), label=tr1$target) val<-tr[h,] gc() dval <- xgb.DMatrix(data.matrix(val[,feature.names]), label=val$target) watchlist <- list(eval = dval, train = dtrain) param <- list( objective = "reg:linear", # booster = "gblinear", min_child_weight = 101, eta = .01, max_depth = 13, # changed from default of 8 subsample = 0.75, colsample_bytree = 0.75, eval_metric = "rmse" # alpha = 0.0001, # lambda = 1 ) clf <- xgb.train( params = param, data = dtrain, nrounds = 500000, verbose = 1, early.stop.round = 100, print.every.n = 100, watchlist = watchlist, maximize = FALSE, nthread = 8) testFeaG <- as.data.frame(test1G) testFeaG$Demanda_uni_equil <- NA # use this model to make predictions on the test set where the best covariate isn't missing for (rows in split(1:nrow(testFeaG), ceiling((1:nrow(testFeaG))/10000))) { testFeaG[rows, "Demanda_uni_equil"] <- exp(predict(clf, data.matrix(testFeaG[rows,feature.names]),ntreelimit = clf$bestInd )) - 1 } ## now we make a model for when M3 is missing. Again, maybe we could impute this instead. feature.names1 <- names(tr)[-c(1,2,3,5,6,7,16,17)] gc() dtrain <- xgb.DMatrix(data.matrix(tr1[,feature.names1]), label=tr1$target) dval <- xgb.DMatrix(data.matrix(val[,feature.names1]), label=val$target) watchlist <- list(eval = dval, train = dtrain) param <- list( objective = "reg:linear", # booster = "gblinear", min_child_weight = 101, eta = .01, max_depth = 8, # changed from default of 8 subsample = 0.8, colsample_bytree = 0.8, eval_metric = "rmse" # alpha = 0.0001, # lambda = 1 ) clf1 <- xgb.train( params = param, data = dtrain, nrounds = 500000, verbose = 1, early.stop.round = 100, print.every.n = 50, watchlist = watchlist, maximize = FALSE, nthread = 8) # where median by product is null use overall median testFeaB <- as.data.frame(test1B) testFeaB$Demanda_uni_equil <- NA # this kind of loops solves some memory problems for (rows in split(1:nrow(testFeaB), ceiling((1:nrow(testFeaB))/10000))) { testFeaB[rows, "Demanda_uni_equil"] <- exp(predict(clf1, data.matrix(testFeaB[rows,feature.names1]),ntreelimit = clf1$bestInd )) - 1 } # bind the two test sets together testFea <- rbind(testFeaB,testFeaG) submit <- testFea[,c(4,16)] for (rows in split(1:nrow(submit), ceiling((1:nrow(submit))/10000))) { ifelse(submit[rows,"Demanda_uni_equil"] < 0,.001,submit[rows,"Demanda_uni_equil"]) }