# Custom logistic regression model

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.

    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

    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.

# 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)`:

    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):

    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:

    // [[Rcpp::export]]
    SEXP make_my_cost_pelt_xptr() {
      Rcpp::XPtr<CostPeltFnPtr> 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)`:

    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:

    R -e 'knitr::knit("vignettes/examples-custom-model.Rmd.original", output = "vignettes/examples-custom-model.Rmd")'

# Appendix: all code snippets

    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)