{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"suppressMessages({\n    library(arrow)\n    library(dplyr)\n    library(ggplot2)\n    library(ROCR)\n    library(duckdb)\n    library(Rcpp)\n    library(tsibble)\n    library(slider)\n})","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-11-10T01:57:06.410589Z","iopub.execute_input":"2023-11-10T01:57:06.414253Z","iopub.status.idle":"2023-11-10T01:57:09.224758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"source(\"../usr/lib/competition_r_code/competition_r_code.R\")","metadata":{"execution":{"iopub.status.busy":"2023-11-10T01:57:12.858310Z","iopub.execute_input":"2023-11-10T01:57:12.900392Z","iopub.status.idle":"2023-11-10T01:57:36.266507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load the data\n\nLoad previously pre-processed data with sleep and awake spans labeled as well as a few features added.","metadata":{}},{"cell_type":"code","source":"dir(\"../input/enrich-split-and-save-train-data-set.\")\ndtr1e <- read_parquet(\"../input/enrich-split-and-save-train-data-set/preproc_train_series_1.parquet\")\ndtr4e <- read_parquet(\"../input/enrich-split-and-save-train-data-set/preproc_train_series_4.parquet\")","metadata":{"execution":{"iopub.status.busy":"2023-11-10T01:57:39.150498Z","iopub.execute_input":"2023-11-10T01:57:39.152519Z","iopub.status.idle":"2023-11-10T01:58:26.259791Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary(dtr1e)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:00:41.377924Z","iopub.execute_input":"2023-11-05T21:00:41.379529Z","iopub.status.idle":"2023-11-05T21:01:43.508236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Select train and test subsets","metadata":{}},{"cell_type":"markdown","source":"# Linear models: SVM, logistic regression","metadata":{}},{"cell_type":"code","source":"aa_set <- which(dtr1e$status %in% c(\"asleep\",\"awake\") & \n                !is.na(dtr1e$enmo) & !is.na(dtr1e$anglez_diff) & \n                !is.na(dtr1e$enmo_sd10m) & !is.na(dtr1e$anglez_sd10m)\n               )\nlength(aa_set)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:02:29.173867Z","iopub.execute_input":"2023-11-05T21:02:29.175776Z","iopub.status.idle":"2023-11-05T21:02:33.865626Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_set <- sample(aa_set, 500000)\ntest_set <- sample(setdiff(aa_set, train_set), 500000)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:02:41.251090Z","iopub.execute_input":"2023-11-05T21:02:41.252354Z","iopub.status.idle":"2023-11-05T21:02:43.272846Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(LiblineaR)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:02:48.673557Z","iopub.execute_input":"2023-11-05T21:02:48.674680Z","iopub.status.idle":"2023-11-05T21:02:48.692121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All available features.","metadata":{}},{"cell_type":"code","source":"fml8 <- status~0+hour+enmo+enmo_diff+sqrt(enmo)+enmo_sd1m+enmo_sd2m+enmo_sd10m+enmo_sqrt_sd1m+enmo_sqrt_sd10m+\nanglez+anglez_diff+sqrt(abs(anglez))+anglez_sd1m+anglez_sd2m+anglez_sd10m+anglez_sqrt_sd1m+anglez_sqrt_sd10m\n\nx <- model.matrix(fml8, dtr1e[train_set,])\nx <- scale(x)\ny <- ifelse(dtr1e$status[train_set]==\"asleep\",0,1)","metadata":{"execution":{"iopub.status.busy":"2023-11-10T01:58:58.864611Z","iopub.execute_input":"2023-11-10T01:58:58.867642Z","iopub.status.idle":"2023-11-10T01:59:26.978931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"See cross-validated model accuracy for model types:\n\n- 0 – L2-regularized logistic regression (primal)\n- 2 – L2-regularized L2-loss support vector classification (primal)\n- 6 – L1-regularized logistic regression","metadata":{}},{"cell_type":"code","source":"cat(\"model type 0 accuracy\", LiblineaR(data=x,target=y,type=0, cost=1, cross=5), \"\\n\")\ncat(\"model type 2 accuracy\", LiblineaR(data=x,target=y,type=2, cost=1, cross=5), \"\\n\")\ncat(\"model type 6 accuracy\", LiblineaR(data=x,target=y,type=6, cost=1, cross=5), \"\\n\")","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:03:04.134875Z","iopub.execute_input":"2023-11-05T21:03:04.136198Z","iopub.status.idle":"2023-11-05T21:03:35.358738Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# XGBoost","metadata":{}},{"cell_type":"code","source":"suppressMessages(library(xgboost))","metadata":{"execution":{"iopub.status.busy":"2023-11-10T01:58:45.565017Z","iopub.execute_input":"2023-11-10T01:58:45.566794Z","iopub.status.idle":"2023-11-10T01:58:47.103399Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Allow NAs in the data this time","metadata":{}},{"cell_type":"code","source":"aa_set <- which(dtr1e$status %in% c(\"asleep\",\"awake\"))\ntrain_set <- sample(aa_set, 1000000)\ny <- ifelse(dtr1e$status[train_set]==\"asleep\",0,1)\nx <- model.matrix(fml8, model.frame(fml8,dtr1e[train_set,], na.action=na.pass))","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:27:01.072413Z","iopub.execute_input":"2023-11-10T02:27:01.075387Z","iopub.status.idle":"2023-11-10T02:27:06.905046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_set_series_ids <- unique(dtr1e$series_id[train_set])\nsave(train_set_series_ids, file = \"train_set_series_ids.RData\")","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:27:09.739503Z","iopub.execute_input":"2023-11-10T02:27:09.741772Z","iopub.status.idle":"2023-11-10T02:27:09.764761Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary(x)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:37:53.842385Z","iopub.execute_input":"2023-11-05T21:37:53.843695Z","iopub.status.idle":"2023-11-05T21:37:59.005065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cv_res <- xgb.cv(data=x, label=y, \n       nrounds=300, nfold=5, metrics=list(\"error\",\"auc\"), max_depth=5, eta=1, \n       objective = \"binary:logistic\", early_stopping_rounds=3,\n      verbose=FALSE)\ncv_res$evaluation_log[cv_res$best_iteration]","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:27:13.045878Z","iopub.execute_input":"2023-11-10T02:27:13.047740Z","iopub.status.idle":"2023-11-10T02:27:14.015311Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"param <- list(max_depth = 5, eta = 1, objective = \"binary:logistic\", eval_metric = \"error\")\nbst5 <- xgb.train(param, data=xgb.DMatrix(x, label=y), nrounds = cv_res$best_iteration)\nxgb.save(bst5,\"xgboost_sleep_state_5.model\")","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:27:18.199274Z","iopub.execute_input":"2023-11-10T02:27:18.201404Z","iopub.status.idle":"2023-11-10T02:27:18.388648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb.importance(model=bst5)","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:27:20.672217Z","iopub.execute_input":"2023-11-10T02:27:20.674244Z","iopub.status.idle":"2023-11-10T02:27:20.728925Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Verify against some sample series","metadata":{}},{"cell_type":"code","source":"sample_ser <- sample(levels(dtr4e$series_id)[which(table(dtr4e$series_id) > 0)],4)","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:28:11.989112Z","iopub.execute_input":"2023-11-10T02:28:11.991639Z","iopub.status.idle":"2023-11-10T02:28:13.896208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ser1_ix = dtr4e$series_id == sample_ser[1] &  !is.na(dtr4e$status)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:46.761005Z","iopub.execute_input":"2023-11-05T21:07:46.763013Z","iopub.status.idle":"2023-11-05T21:07:47.813232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xt <- model.matrix(fml8, model.frame(fml8, dtr4e[ser1_ix,], na.action=na.pass))\nyt <- ifelse(dtr4e$status[ser1_ix]==\"asleep\",0,1)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:47.816623Z","iopub.execute_input":"2023-11-05T21:07:47.818359Z","iopub.status.idle":"2023-11-05T21:07:51.537845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"yp <- predict(bst5, newdata=xt)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:51.541279Z","iopub.execute_input":"2023-11-05T21:07:51.543005Z","iopub.status.idle":"2023-11-05T21:07:51.637542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot(performance(prediction(yp, yt), \"tpr\", \"fpr\"))","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:51.640949Z","iopub.execute_input":"2023-11-05T21:07:51.642528Z","iopub.status.idle":"2023-11-05T21:07:52.261336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot(performance(prediction(yp, yt), \"acc\"))","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:52.264486Z","iopub.execute_input":"2023-11-05T21:07:52.266424Z","iopub.status.idle":"2023-11-05T21:07:52.752093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"performance(prediction(yp, yt), \"auc\")@y.values[[1]]","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:52.755738Z","iopub.execute_input":"2023-11-05T21:07:52.757753Z","iopub.status.idle":"2023-11-05T21:07:53.099086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prop.table(table(pred=yp>0.5,actual=yt), 2)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:53.102181Z","iopub.execute_input":"2023-11-05T21:07:53.103739Z","iopub.status.idle":"2023-11-05T21:07:53.231153Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean((yp>0.5) == (yt == 1))","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:53.233935Z","iopub.execute_input":"2023-11-05T21:07:53.235140Z","iopub.status.idle":"2023-11-05T21:07:53.248635Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Look at the xgb scores distribution over the timeline","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\n\nser = dtr4e[ser1_ix,]\nser$pred_score <- yp\nser_change <- ser[!is.na(ser$event),]\nser %>% filter(step<100000) %>% \nggplot() + geom_point(aes(timestamp, pred_score, col=status)) + \ngeom_vline(aes(xintercept = timestamp, colour = event), ser_change)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:07:56.565905Z","iopub.execute_input":"2023-11-05T21:07:56.567338Z","iopub.status.idle":"2023-11-05T21:07:59.341992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Verify against all other series","metadata":{}},{"cell_type":"code","source":"model_perf <- function(dtr) {\n  sers <- levels(dtr$series_id)[which(table(dtr$series_id) > 0)]\n  sapply(sers, function(s) {\n      dds <- dtr[dtr$series_id==s,]\n      ser_ix = !is.na(dds$status)\n      xt <- model.matrix(fml8, model.frame(fml8, dds[ser_ix,], na.action=na.pass))\n      yt <- ifelse(dds$status[ser_ix]==\"asleep\",0,1)\n      yp <- predict(bst5, newdata=xt)\n      auc <- performance(prediction(yp, yt), \"auc\")@y.values[[1]]\n      acc <- mean((yp>0.5) == (yt == 1))\n      c(series_id=s, accuracy=acc, AUC=auc)\n  }) %>% t()\n}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr4_perf <- as.data.frame(model_perf(dtr4e))","metadata":{"execution":{"iopub.status.idle":"2023-10-31T05:28:59.664214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dtr4_perf$accuracy <- as.numeric(dtr4_perf$accuracy)\ndtr4_perf$AUC <- as.numeric(dtr4_perf$AUC)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T05:28:59.669022Z","iopub.execute_input":"2023-10-31T05:28:59.671560Z","iopub.status.idle":"2023-10-31T05:28:59.690362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"summary(dtr4_perf)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T05:28:59.694621Z","iopub.execute_input":"2023-10-31T05:28:59.696314Z","iopub.status.idle":"2023-10-31T05:28:59.717118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot(dtr4_perf) + geom_boxplot(aes(accuracy)) ","metadata":{"execution":{"iopub.status.busy":"2023-10-31T05:28:59.721323Z","iopub.execute_input":"2023-10-31T05:28:59.722966Z","iopub.status.idle":"2023-10-31T05:29:00.015274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Look at the worst accuracy series","metadata":{}},{"cell_type":"code","source":"bad_ser <- dtr4_perf$series_id[which.min(dtr4_perf$accuracy)]\nser_ix = dtr4e$series_id == bad_ser & !is.na(dtr4e$status)\nxt <- model.matrix(fml8, model.frame(fml8,dtr4e[ser_ix,],na.action=na.pass))\nyt <- ifelse(dtr4e$status[ser_ix]==\"asleep\",0,1)\nyp <- predict(bst5, newdata=xt)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T05:29:00.018737Z","iopub.execute_input":"2023-10-31T05:29:00.020582Z","iopub.status.idle":"2023-10-31T05:29:02.008498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prop.table(table(pred=yp>0.5,actual=yt), 2)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T05:29:02.012720Z","iopub.execute_input":"2023-10-31T05:29:02.014946Z","iopub.status.idle":"2023-10-31T05:29:02.213193Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ser = dtr4e[ser_ix,]\nser$pred_score <- yp\nser_change <- ser[!is.na(ser$event),]\noptions(repr.plot.width=12, repr.plot.height=6)\nser %>% \nggplot() + geom_point(aes(timestamp, pred_score, col=status)) + \ngeom_vline(aes(xintercept = timestamp, colour = event), ser_change)","metadata":{"execution":{"iopub.status.busy":"2023-10-31T05:29:02.216559Z","iopub.execute_input":"2023-10-31T05:29:02.218356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# State undefined prediction","metadata":{}},{"cell_type":"code","source":"summary(dtr1e$status)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:10:22.667220Z","iopub.execute_input":"2023-11-05T21:10:22.668627Z","iopub.status.idle":"2023-11-05T21:10:23.835487Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"state_known <- which(!is.na(dtr1e$status))\nstate_unknown <- which(is.na(dtr1e$status))\nstate_model_train_set <- c(sample(state_known, 500000), \n                           sample(state_unknown,500000))","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:18:46.757946Z","iopub.execute_input":"2023-11-05T21:18:46.759864Z","iopub.status.idle":"2023-11-05T21:18:48.328912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"state_fml82 <- ~0+enmo+enmo_diff+sqrt(enmo)+enmo_sd1m+enmo_sd2m+enmo_sd10m+enmo_sqrt_sd1m+enmo_sqrt_sd10m+anglez+anglez_diff+sqrt(abs(anglez))+anglez_sd1m+anglez_sd2m+anglez_sd10m+anglez_sqrt_sd1m+anglez_sqrt_sd10m\ny <- is.na(dtr1e$status[state_model_train_set])\nx <- model.matrix(state_fml82, \n                  model.frame(state_fml82,dtr1e[state_model_train_set,], \n                              na.action=na.pass)\n                 )\ncv_res <- xgb.cv(data=x, label=y, \n       nrounds=300, nfold=5, metrics=list(\"error\",\"auc\"), max_depth=5, eta=1, \n       objective = \"binary:logistic\", early_stopping_rounds=3,\n      verbose=FALSE)\ncv_res$evaluation_log[cv_res$best_iteration]","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:18:55.813555Z","iopub.execute_input":"2023-11-05T21:18:55.814947Z","iopub.status.idle":"2023-11-05T21:22:16.951601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"param <- list(max_depth = 5, eta = 1, objective = \"binary:logistic\", eval_metric = \"error\")\nstate_bst5 <- xgb.train(param, data=xgb.DMatrix(x, label=y), nrounds = cv_res$best_iteration)\nxgb.save(state_bst5,\"xgboost_unknown_state_5.model\")","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:22:48.133787Z","iopub.execute_input":"2023-11-05T21:22:48.135036Z","iopub.status.idle":"2023-11-05T21:23:26.230010Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"xgb.importance(model=state_bst5)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:23:26.233196Z","iopub.execute_input":"2023-11-05T21:23:26.234822Z","iopub.status.idle":"2023-11-05T21:23:26.263443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_ser_id <- sample(unique(dtr4e$series_id[is.na(dtr4e$status)]),1)\nxt <- model.matrix(state_fml82, model.frame(state_fml82, \n                                            dtr4e[dtr4e$series_id == sample_ser_id,], \n                                            na.action=na.pass))\nyt <- is.na(dtr4e$status[dtr4e$series_id == sample_ser_id])\nyp <- predict(state_bst5, newdata=xt)","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:27:48.961158Z","iopub.execute_input":"2023-11-10T02:27:48.963068Z","iopub.status.idle":"2023-11-10T02:27:51.763011Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot(performance(prediction(yp, yt), \"tpr\", \"fpr\"))","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:29:35.246988Z","iopub.execute_input":"2023-11-05T21:29:35.248311Z","iopub.status.idle":"2023-11-05T21:29:35.651251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot(performance(prediction(yp, yt), \"prec\", \"rec\"))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot(performance(prediction(yp, yt), \"f\"))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot(performance(prediction(yp, yt), \"acc\"))","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:30:20.457940Z","iopub.execute_input":"2023-11-05T21:30:20.459257Z","iopub.status.idle":"2023-11-05T21:30:21.015532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prop.table(table(pred=yp>0.5,actual=yt), 2)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:31:31.760689Z","iopub.execute_input":"2023-11-05T21:31:31.762023Z","iopub.status.idle":"2023-11-05T21:31:31.898593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean((yp>0.5) == yt)","metadata":{"execution":{"iopub.status.busy":"2023-11-05T21:30:56.642879Z","iopub.execute_input":"2023-11-05T21:30:56.644235Z","iopub.status.idle":"2023-11-05T21:30:56.661709Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Find change point\n\nUse adjacent time windows to detect if that time was in sleep, awaken or a change from one state to another. Use sum of xgboost model scores to decide. Use one randomly selected series to work on as a sample.","metadata":{}},{"cell_type":"code","source":"dd = dtr4e[dtr4e$series_id == sample_ser[1],]\nfml82 <- ~0+hour+enmo+enmo_diff+sqrt(enmo)+enmo_sd1m+enmo_sd2m+enmo_sd10m+enmo_sqrt_sd1m+enmo_sqrt_sd10m+anglez+anglez_diff+sqrt(abs(anglez))+anglez_sd1m+anglez_sd2m+anglez_sd10m+anglez_sqrt_sd1m+anglez_sqrt_sd10m\nxt <- model.matrix(fml82, model.frame(fml82,dd,na.action=na.pass))\nyt <- ifelse(dd$status==\"awake\",1,0)\nyp <- predict(bst5, newdata=xt)\ndd$pred_score <- yp","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:28:28.699019Z","iopub.execute_input":"2023-11-10T02:28:28.700996Z","iopub.status.idle":"2023-11-10T02:28:38.017583Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ps <- dd %>% arrange(step) %>% pull(pred_score)\nstep <- dd %>% arrange(step) %>% pull(step)\npred_change <- change_prob_refine(ps, step, 2000)","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:28:41.720876Z","iopub.execute_input":"2023-11-10T02:28:41.722991Z","iopub.status.idle":"2023-11-10T02:28:43.106398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dd <- merge(dd, pred_change, by=\"step\", all=TRUE)\nnames(dd)[names(dd)==\"event.y\"] <- \"pred_change\"\nnames(dd)[names(dd)==\"event.x\"] <- \"event\"\ndd$pred_event <- as.factor(dd$pred_change)","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:28:46.782679Z","iopub.execute_input":"2023-11-10T02:28:46.784721Z","iopub.status.idle":"2023-11-10T02:28:48.368162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"See if predicted and actual change points are close","metadata":{}},{"cell_type":"code","source":"merge(pred_change, dd[!is.na(dd$event), c(\"step\",\"event\")],all=TRUE,by=\"step\") %>% \n  head(20)","metadata":{"execution":{"iopub.status.busy":"2023-11-10T02:28:51.450896Z","iopub.execute_input":"2023-11-10T02:28:51.452760Z","iopub.status.idle":"2023-11-10T02:28:51.499540Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"See if visually change points seems to be at the right place","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width=12, repr.plot.height=6)\n\ndd %>% filter(step<100000) %>% \n  ggplot() + geom_point(aes(timestamp, pred_score, col=status)) + \n  geom_vline(aes(xintercept = timestamp, colour=pred_change), dd[!is.na(dd$pred_change),])","metadata":{"execution":{"iopub.status.busy":"2023-11-04T18:48:55.235737Z","iopub.execute_input":"2023-11-04T18:48:55.237855Z","iopub.status.idle":"2023-11-04T18:48:59.196236Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}