--- title: "Custom logistic regression model" output: rmarkdown::html_vignette description: | Examples of using custom model to reproduce change point detection in logistic regression models. vignette: > %\VignetteIndexEntry{Custom logistic regression model} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- We try to reproduce the logistic regression models with custom cost functions and show that the results are similar to the built-in logistic regression models. The built-in logistic regression model in `fastcpd()` is implemented with the help of the **fastglm** package. The **fastglm** package utilizes the iteratively reweighted least squares with the step-halving approach to help safeguard against convergence issues. If a custom cost function is used with gradient descent, we should expect the results will be similar to the built-in logistic regression model. Specifying the `cost`, `cost_gradient` and `cost_hessian` parameters below, we can obtain similar results as the built-in logistic regression model. ``` r set.seed(1) x <- matrix(rnorm(2000, 0, 1), ncol = 5) theta <- rbind(rnorm(5, 0, 1), rnorm(5, 2, 1)) y <- c( rbinom(250, 1, 1 / (1 + exp(-x[1:250, ] %*% theta[1, ]))), rbinom(150, 1, 1 / (1 + exp(-x[251:400, ] %*% theta[2, ]))) ) binomial_data <- data.frame(y = y, x = x) result <- fastcpd.binomial(cbind(y, x), r.progress = FALSE, cost_adjustment = NULL) summary(result) #> #> Call: #> fastcpd.binomial(data = cbind(y, x), r.progress = FALSE, cost_adjustment = NULL) #> #> Change points: #> 250 #> #> Cost values: #> 99.40994 36.89336 #> #> Parameters: #> segment 1 segment 2 #> 1 -1.0096300 3.46131278 #> 2 -1.7740956 2.10636392 #> 3 1.2736663 0.74124627 #> 4 0.3955351 0.07297997 #> 5 -0.1047536 2.00553153 ``` ``` r logistic_loss <- function(data, theta) { x <- data[, -1] y <- data[, 1] u <- x %*% theta nll <- -y * u + log(1 + exp(u)) nll[u > 10] <- -y[u > 10] * u[u > 10] + u[u > 10] sum(nll) } logistic_gradient <- function(data, theta) { x <- data[nrow(data), -1] y <- data[nrow(data), 1] c(-(y - 1 / (1 + exp(-x %*% theta)))) * x } logistic_hessian <- function(data, theta) { x <- data[nrow(data), -1] prob <- 1 / (1 + exp(-x %*% theta)) (x %o% x) * c((1 - prob) * prob) } result <- fastcpd( y ~ . - 1, binomial_data, epsilon = 1e-5, cost = logistic_loss, cost_gradient = logistic_gradient, cost_hessian = logistic_hessian, r.progress = FALSE ) summary(result) #> #> Call: #> fastcpd(formula = y ~ . - 1, data = binomial_data, cost = logistic_loss, #> cost_gradient = logistic_gradient, cost_hessian = logistic_hessian, #> epsilon = 1e-05, r.progress = FALSE) #> #> Change points: #> 252 #> #> Parameters: #> segment 1 segment 2 #> 1 0 0 #> 2 0 0 #> 3 0 0 #> 4 0 0 #> 5 0 0 ``` Note that the result obtained through custom cost functions is inferior compared to the one obtained through built-in models. We remark that the results can be improved with extra parameters already provided in the package. The detailed discussion of several advanced usages of the package can be found in [Advanced examples](examples-advanced.html). # Compiled (C++) custom cost functions The custom cost functions above are R closures: `cost`, `cost_gradient` and `cost_hessian` are called from C++ once per candidate segment via `Rcpp::Function`, which incurs R call overhead. For performance-critical use cases, `fastcpd()` also accepts **pre-compiled C++ cost functions**, passed through the very same `cost` / `cost_gradient` / `cost_hessian` arguments as an `externalptr` built with `Rcpp::XPtr`. The compiled and R-closure paths compute the same mathematical cost and a `family = "custom"` run accepts either interchangeably -- `fastcpd()` detects which kind of SEXP it received (`TYPEOF(x) == EXTPTRSXP`) and dispatches accordingly, with no other change to the call. To build a compiled cost, write a C++ function matching one of the two expected signatures -- a PELT-style `cost(data)`: ```cpp double my_cost_pelt(arma::mat const& data); typedef double (*CostPeltFnPtr)(arma::mat const&); ``` or a SeGD-style `cost(data, theta)` (optionally paired with a gradient and a Hessian): ```cpp double my_cost_sen(arma::mat const& data, arma::colvec const& theta); typedef double (*CostSenFnPtr)(arma::mat const&, arma::colvec const&); arma::colvec my_cost_gradient(arma::mat const& data, arma::colvec const& theta); typedef arma::colvec (*CostGradientFnPtr)(arma::mat const&, arma::colvec const&); arma::mat my_cost_hessian(arma::mat const& data, arma::colvec const& theta); typedef arma::mat (*CostHessianFnPtr)(arma::mat const&, arma::colvec const&); ``` Wrap the function pointer in an `Rcpp::XPtr`, tagging it with the literal string `fastcpd` expects for that role (`"fastcpd_cost_pelt"`, `"fastcpd_cost_sen"`, `"fastcpd_cost_gradient"` or `"fastcpd_cost_hessian"` -- this tag is how `fastcpd` validates that the pointer was built for the right purpose), and export a function that returns it: ```cpp // [[Rcpp::export]] SEXP make_my_cost_pelt_xptr() { Rcpp::XPtr xptr( new CostPeltFnPtr(&my_cost_pelt), true, Rcpp::wrap("fastcpd_cost_pelt")); return xptr; } ``` Because external pointers carry no `formals()`, `fastcpd()` cannot infer whether a compiled `cost` is PELT-style (one argument) or SeGD-style (two arguments) the way it does for R closures. Instead, set an integer `fastcpd_cost_arity` attribute on the returned pointer -- `1L` for `cost(data)`, `2L` for `cost(data, theta)`: ``` r my_cost_pelt_xptr <- make_my_cost_pelt_xptr() attr(my_cost_pelt_xptr, "fastcpd_cost_arity") <- 1L result <- fastcpd( ~ . - 1, data.frame(x = data), family = "custom", cost = my_cost_pelt_xptr, r.progress = FALSE ) summary(result) ``` A compiled `cost` can be used on its own (covering both the PELT-style and SeGD-style cases above), or paired with compiled `cost_gradient` and `cost_hessian` `externalptr`s the same way -- including the gradient/Hessian- driven PELT warm start (`fastcpd_cost_arity = 1L`), which is now solved by a from-scratch Brent / bound-constrained BFGS optimiser working directly through the unified `cost`/`cost_gradient`/`cost_hessian` wrappers, with no R-level callback and therefore no restriction on whether `cost` is an R closure or a compiled `externalptr`. # Notes This document is generated by the following code: ```shell R -e 'knitr::knit("vignettes/examples-custom-model.Rmd.original", output = "vignettes/examples-custom-model.Rmd")' ``` # Appendix: all code snippets ``` r knitr::opts_chunk$set( collapse = TRUE, comment = "#>", eval = TRUE, warning = FALSE ) library(fastcpd) set.seed(1) x <- matrix(rnorm(2000, 0, 1), ncol = 5) theta <- rbind(rnorm(5, 0, 1), rnorm(5, 2, 1)) y <- c( rbinom(250, 1, 1 / (1 + exp(-x[1:250, ] %*% theta[1, ]))), rbinom(150, 1, 1 / (1 + exp(-x[251:400, ] %*% theta[2, ]))) ) binomial_data <- data.frame(y = y, x = x) result <- fastcpd.binomial(cbind(y, x), r.progress = FALSE, cost_adjustment = NULL) summary(result) logistic_loss <- function(data, theta) { x <- data[, -1] y <- data[, 1] u <- x %*% theta nll <- -y * u + log(1 + exp(u)) nll[u > 10] <- -y[u > 10] * u[u > 10] + u[u > 10] sum(nll) } logistic_gradient <- function(data, theta) { x <- data[nrow(data), -1] y <- data[nrow(data), 1] c(-(y - 1 / (1 + exp(-x %*% theta)))) * x } logistic_hessian <- function(data, theta) { x <- data[nrow(data), -1] prob <- 1 / (1 + exp(-x %*% theta)) (x %o% x) * c((1 - prob) * prob) } result <- fastcpd( y ~ . - 1, binomial_data, epsilon = 1e-5, cost = logistic_loss, cost_gradient = logistic_gradient, cost_hessian = logistic_hessian, r.progress = FALSE ) summary(result) my_cost_pelt_xptr <- make_my_cost_pelt_xptr() attr(my_cost_pelt_xptr, "fastcpd_cost_arity") <- 1L result <- fastcpd( ~ . - 1, data.frame(x = data), family = "custom", cost = my_cost_pelt_xptr, r.progress = FALSE ) summary(result) ```