{"cells":[{"metadata":{},"cell_type":"markdown","source":"This kernel extends [https://www.kaggle.com/azc2019/earthquake-forecast-with-random-forests](https://www.kaggle.com/azc2019/earthquake-forecast-with-random-forests) by using XGBoost instead of Random Forests. \n\nThanks to Michael Mayer for the original kernel - appreciate your work!\n\nI've seen plenty of good <1.5 scores in Python but I prefer R as my workhorse so sharing this kernel with other R users, as they shared theirs.\n\nPlease upvote Michael's original kernel if you find this useful. Upvoting mine would be appreciated as well.\n\n* Added additional features based on papers provided by LANL\n* Changed rolling windows\n* Removed disk saving and loading\n* Added tuned XGBoost parameters\n\nThere is scope to improve score specifically around choosing validation folds and if you do, please share your kernel as well."},{"metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"cell_type":"code","source":"## Importing packages\n\n# This R environment comes with all of CRAN and many other helpful packages preinstalled.\n# You can see which packages are installed by checking out the kaggle/rstats docker image: \n# https://github.com/kaggle/docker-rstats\n\noptions(digits = 15)\n\nlibrary(data.table)\nlibrary(xgboost)\nlibrary(tidyverse)\nlibrary(dplyr)\nlibrary(e1071)\nlibrary(roll)\n\nrateChange <- function(x) ((last(x)-first(x))/first(x))*100\n\n# Calculates a bunch of statistics from vector x\nunivariate_stats <- function(x, tag = NULL, p = c(0, 0.25, 0.75, 1)) {\n  x <- x[!is.na(x)]\n  \n  out <- c(\n    min = min(x),\n    max = max(x),\n    range = max(x) - min(x),\n    mean = mean(x),\n    sd = sd(x),\n    sk = skewness(x),\n    kr = kurtosis(x),\n    vr = var(x),\n    mad = mad(x),\n    chn = rateChange(x),\n    setNames(quantile(x, p = p, names = FALSE), paste0(\"q\", p))\n    )\n  \n  if (is.null(tag)) \n    return(out)\n  \n  names(out) <- paste(names(out), tag, sep = \"_\")\n  out\n}\n\n# Feature extraction\ncreate_X <- function(x, rolling_windows = c(10,100,1000,10000)) {\n  stats_full <- univariate_stats(x = x, \"full\", p = c(0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,\n                                                      0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99))\n\n  stats_abs <- univariate_stats(x = abs(x), \"abs\")\n  \n  # Rolling versions of x\n  x_mat <- as.matrix(x, ncol = 1)\n  \n  roll_sd_k <- lapply(rolling_windows, function(k) roll_sd(x_mat, width = k))\n  stats_roll_sd <- Map(univariate_stats, roll_sd_k, tag = paste(\"roll_sd\", rolling_windows, sep = \"_\"))\n  \n  roll_mean_k <- lapply(rolling_windows, function(k) roll_mean(x_mat, width = k))\n  stats_roll_mean <- Map(univariate_stats, roll_mean_k, tag = paste(\"roll_mean\", rolling_windows, sep = \"_\"))\n  \n  c(stats_full, stats_abs, unlist(stats_roll_sd), unlist(stats_roll_mean))\n}\n\n# Length of test data sets\nn_test <- 150000\n\n# By how much do we shift the time window of n_test rows within earthquake?\nstride <- 75000\n\n# Positions of earthquakes (see answers in https://www.kaggle.com/c/LANL-Earthquake-Prediction/discussion/77390). \n# Used to create contiguous folds for cross-validation and train/validation split\nearthquakes <- c(\n  5656573,\n  50085877,\n  104677355,\n  138772452,\n  187641819,\n  218652629,\n  245829584,\n  307838916,\n  338276286,\n  375377847,\n  419368879,\n  461811622,\n  495800224,\n  528777114,\n  585568143,\n  621985672) + 1\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This runs for quite a while."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Figure out colnames of input as well as number of features\nraw <- fread(file.path(\"..\", \"input\", \"train.csv\"), nrows = 0, data.table = FALSE)\nnames_input <- names(raw)\nnames_features <- names(create_X(raw[[1]]))\n\ntrain_x <- matrix(NA, nrow = 0, ncol = length(names_features), dimnames = list(NULL, names_features))\ntrain_y <- numeric()\n\nfor (s in seq_along(earthquakes)) {\n  cat(s, \"\\n\")\n  \n  # Read data between two earthquakes\n  raw <- fread(file.path(\"..\", \"input\", \"train.csv\"), \n               nrows = c(earthquakes[1], diff(earthquakes))[s], \n               skip = c(0, earthquakes)[s] + 1)\n  setnames(raw, names_input)\n  \n  # How many times do we calculate features for this data chunk?\n  n_steps <- (nrow(raw) - n_test) %/% stride\n  \n  # Init feature matrix and vector of response\n  y <- numeric(n_steps)\n  X <- matrix(NA, nrow = n_steps, ncol = length(names_features), dimnames = list(NULL, names_features))\n  \n  # Loop through chunk and build up y and X\n  pb <- txtProgressBar(0, n_steps, style = 3)\n  \n  for (i in seq_len(n_steps)) {\n    setTxtProgressBar(pb, i)\n    from <- 1 + stride * (i - 1)\n    to <- n_test + stride * (i - 1)\n    X[i, ] <- create_X(raw$acoustic_data[from:to])\n    y[i] <- raw$time_to_failure[to]\n  }\n  \n  train_x <- rbind(train_x, X)\n  train_y <- c(train_y, y)\n}\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"I've tried Catboost and LightGBM but decided to stick with XGBoost as the paper is based on XGBoost. Parameters was tuned with gridsearch and set in this markdown. Depth of 5 and 6 improves the accuracy during training but might overfit so sticking with 4."},{"metadata":{"trusted":true},"cell_type":"code","source":"# XGBoost\ndtrain <- xgb.DMatrix(data = train_x, label = train_y) \n\nparam <- list(eta = 0.003,\n              max_depth=4,  \n              subsample=0.5,\n              colsample_bytree=0.5,\n              alpha=0.1,\n              min_samples_split=10,\n              min_samples_leaf=10,\n              max_leaf_nodes=10,\n              objective='reg:linear',\n              random_state=42,\n              verbose = TRUE,\n              eval_metric = \"mae\",\n              booster = \"gbtree\"\n)\n\nseed.number  <-  sample.int(10000, 1)\nset.seed(seed.number)\ncv.nround = 50000\ncv.nfold = 5\nmdcv <- xgb.cv(data=dtrain, params = param,  \n               nfold=cv.nfold, nrounds=cv.nround,\n               verbose = T,early_stopping_rounds = 10, shuffle=TRUE, stratified = TRUE\n)\n\nset.seed(42)\nmodel.xgb <- xgboost(data = dtrain,\n                     nrounds = mdcv$best_iteration,\n                     print_every_n = 200,\n                     params = param\n)\n\n# Remove as required based on below\nimp<-xgb.importance(colnames(X), model = model.xgb)\nprint(imp)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Submission\nsubmission <- fread(file.path(\"..\", \"input\", \"sample_submission.csv\"))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# This also takes a long time to run\nload_and_prepare <- function(file) {\n  seg <- fread(file.path(\"..\", \"input\", \"test\", paste0(file, \".csv\")))\n  create_X(seg$acoustic_data)\n}\nall_test <- lapply(submission$seg_id, load_and_prepare)\nall_test2 <- do.call(rbind, all_test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pred.xgb <- predict(model.xgb, all_test2)\n\nsubmission$time_to_failure <- pred.xgb\nfwrite(submission, paste0(\"submission_xgb_stride_\", stride, \".csv\"))","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"R","language":"R","name":"ir"},"language_info":{"mimetype":"text/x-r-source","name":"R","pygments_lexer":"r","version":"3.4.2","file_extension":".r","codemirror_mode":"r"}},"nbformat":4,"nbformat_minor":1}