Package {boostPM}

boostPM: Unsupervised Tree Boosting for Learning Probability Distributions

Description

Implements the unsupervised tree boosting method of Awaya and Ma (2024) for fitting tree-ensemble probability distributions, evaluating fitted densities, and generating samples.

Author(s)

Maintainer: Naoki Awaya nawaya@waseda.jp

References

Awaya, N. and Ma, L. (2024). Unsupervised Tree Boosting for Learning Probability Distributions. Journal of Machine Learning Research, 25, 1–52.

See Also

fit_boostpm(), predict.boostPM_fit(), simulate.boostPM_fit(), and plot.boostPM_fit().

Fit a boostPM Distribution

Description

Fits the unsupervised tree boosting procedure of Awaya and Ma (2024). Marginal trees are fitted first, followed by dependence trees. The fitted ensemble represents a probability distribution on a rectangular support.

Usage

fit_boostpm(
  data,
  Omega = NULL,
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 8,
  max_split_depth = 15,
  min_node_observations = 5,
  c0 = 0.1,
  gamma = 0.1,
  early_stop = NULL,
  prior_split_prob = 0.9,
  add_noise = TRUE,
  progress = c("none", "stage")
)

Arguments

Value

An object of class boostPM_fit, represented as a list containing:

trees

Serialized fitted trees used by package methods. This is an internal representation and is not a stable interface for direct editing.

residual_coordinates

Final residual coordinates as a variables-by-observations numeric matrix. This advanced diagnostic is not intended for direct modification.

tree_diagnostics

A data frame with one row per accepted tree and columns tree_index, stage, variable, node_count, and max_depth. stage is "marginal" or "dependence"; variable identifies the marginal variable and is NA for dependence trees.

variable_importance

Unnormalized accumulated splitting improvement attributed to each variable. Larger values indicate greater contribution within the fitted ensemble. The vector is named when data has column names.

heldout_diagnostics

A data frame with one row per held-out candidate and columns candidate_index, stage, variable, mean_log_density_improvement, and accepted, or NULL when early_stop = NULL. The candidate that triggers stopping has accepted = FALSE.

support

A variables-by-two matrix containing named lower and upper original-scale support limits. Row names follow the training variable names when available.

elapsed_time

Elapsed fitting time as a difftime object.

call

The matched call used to fit the model.

control

A named list of the validated controls actually used: max_marginal_trees, max_dependence_trees, n_bins, max_split_depth, min_node_observations, c0, gamma, early_stop, prior_split_prob, add_noise, and progress.

Fitting stages

The procedure first fits up to max_marginal_trees univariate trees for each variable. It then fits up to max_dependence_trees multivariate trees to the remaining dependence structure. These arguments are upper bounds; adaptive stopping can end a stage earlier.

Scale-specific learning rate

For a node A in the normalized unit cube, the learning rate is

c(A) = c_0 (1 - \log_2[\mathrm{vol}(A)])^{-\gamma}.

Thus c0 determines the global update size. When gamma = 0, every node uses c0. When gamma > 0, smaller-volume nodes receive smaller updates, which imposes stronger shrinkage on local structure.

Adaptive stopping

If early_stop = NULL, all observations fit each candidate tree and no adaptive stopping is applied. With early_stop = c(threshold, window), each candidate tree is fitted to a random 90 percent subset of the current residuals and evaluated by its average log density on the held-out 10 percent. A stage ends when the mean of the most recent window held-out scores is strictly below threshold. The candidate that triggers stopping is recorded in heldout_diagnostics with accepted = FALSE but is not added to the ensemble. c(0, 50) is the paper-oriented setting; the implementation uses a strict comparison with zero, whereas the paper describes non-positive improvement.

Support and tied observations

Data are transformed from Omega to the unit cube before fitting. If Omega = NULL, the range after optional tie perturbation is expanded by 10 percent on each side. Set add_noise = TRUE when a continuous distribution is assumed and ties are technical artifacts such as rounding. Set it to FALSE when perturbing ties would be inappropriate.

Split prior

The default prior_split_prob = 0.9 follows the archived implementation and public experiment code. Setting prior_split_prob = 0.5 corresponds to the constant split prior implied by the 0.5 stopping probability reported in Appendix C of Awaya and Ma (2024). Larger values favor splitting before the data-dependent candidate likelihoods are incorporated.

Reproducibility and progress

Tree construction is stochastic even when add_noise = FALSE, and adaptive stopping also draws random fitting subsets. Call set.seed() before every fit that must be reproducible. Fitting is silent by default; set progress = "stage" for stage-level messages. Elapsed time is stored in the result and can be inspected with print() or summary(). Exact equality across operating systems has not been established.

References

Awaya, N. and Ma, L. (2024). Unsupervised Tree Boosting for Learning Probability Distributions. Journal of Machine Learning Research, 25, 1–52.

Examples

set.seed(42)
n <- 400L
x1 <- stats::rbeta(n, shape1 = 2, shape2 = 5)
x2 <- stats::rbeta(n, shape1 = 2 + 6 * x1, shape2 = 4)
x <- cbind(x1 = x1, x2 = x2)
support <- cbind(lower = c(0, 0), upper = c(1, 1))

set.seed(123)
fit <- fit_boostpm(
  x,
  Omega = support,
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  add_noise = FALSE
)
print(fit)

grid <- seq(0.01, 0.99, length.out = 50L)
evaluation_grid <- as.matrix(expand.grid(x1 = grid, x2 = grid))
fitted_density <- matrix(
  predict(fit, newdata = evaluation_grid, type = "density"),
  nrow = length(grid)
)
graphics::image(grid, grid, fitted_density, xlab = "x1", ylab = "x2")
graphics::contour(
  grid, grid, fitted_density, add = TRUE, drawlabels = FALSE
)

generated <- simulate(fit, nsim = 500, seed = 321)
graphics::plot(
  generated, xlab = "x1", ylab = "x2", main = "Generated observations"
)

Plot Diagnostics for a Fitted boostPM Distribution

Description

Displays a bar plot of variable importance, the number of nodes per tree, or maximum tree depths.

Usage

## S3 method for class 'boostPM_fit'
plot(
  x,
  type = c("variable_importance", "tree_node_counts", "tree_depths"),
  ...
)

Arguments

Value

The bar midpoints returned by graphics::barplot(), invisibly.

Examples

set.seed(42)
n <- 400L
x1 <- stats::rbeta(n, shape1 = 2, shape2 = 5)
x <- cbind(x1 = x1, x2 = stats::rbeta(n, 2 + 6 * x1, 4))
set.seed(123)
fit <- fit_boostpm(
  x,
  Omega = cbind(lower = c(0, 0), upper = c(1, 1)),
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  add_noise = FALSE
)
plot(fit, type = "variable_importance")
plot(fit, type = "tree_node_counts")
plot(fit, type = "tree_depths")

Evaluate a Fitted boostPM Density

Description

Evaluates the probability density represented by a fitted tree ensemble at supplied points.

Usage

## S3 method for class 'boostPM_fit'
predict(object, newdata, type = c("log_density", "density", "details"), ...)

Arguments

Details

Points outside the fitted rectangular support receive log density -Inf and density zero. At a split point, evaluation follows the left-child convention in the method specification. Because mean_log_density_path averages over every row of newdata, the entire path is -Inf when any row lies outside the support.

Value

For type = "log_density" or type = "density", a numeric vector with one entry per row of newdata, in the original row order. For type = "details", a list containing:

log_density

The same vector returned by type = "log_density".

mean_log_density_path

A numeric vector with one entry per fitted tree. Entry k is the mean log density over all rows of newdata after the first k trees have been applied.

Examples

set.seed(42)
n <- 400L
x1 <- stats::rbeta(n, shape1 = 2, shape2 = 5)
x <- cbind(x1 = x1, x2 = stats::rbeta(n, 2 + 6 * x1, 4))
set.seed(123)
fit <- fit_boostpm(
  x,
  Omega = cbind(lower = c(0, 0), upper = c(1, 1)),
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  add_noise = FALSE
)

grid <- seq(0.01, 0.99, length.out = 50L)
evaluation_grid <- as.matrix(expand.grid(x1 = grid, x2 = grid))
fitted_density <- matrix(
  predict(fit, evaluation_grid, type = "density"),
  nrow = length(grid)
)
graphics::image(grid, grid, fitted_density, xlab = "x1", ylab = "x2")
graphics::contour(
  grid, grid, fitted_density, add = TRUE, drawlabels = FALSE
)

Print a Fitted boostPM Distribution

Description

Prints a compact description of a boostPM_fit object.

Usage

## S3 method for class 'boostPM_fit'
print(x, ...)

Arguments

Value

The input object, invisibly.

Examples

set.seed(42)
n <- 400L
x1 <- stats::rbeta(n, shape1 = 2, shape2 = 5)
x <- cbind(x1 = x1, x2 = stats::rbeta(n, 2 + 6 * x1, 4))
set.seed(123)
fit <- fit_boostpm(
  x,
  Omega = cbind(lower = c(0, 0), upper = c(1, 1)),
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  add_noise = FALSE
)
print(fit)

Print a boostPM Fit Summary

Description

Prints the fitted call, data dimensions, support, structural tree diagnostics, variable importance, and elapsed fitting time.

Usage

## S3 method for class 'summary.boostPM_fit'
print(x, ...)

Arguments

Value

The input summary object, invisibly.

Examples

set.seed(42)
n <- 400L
x1 <- stats::rbeta(n, shape1 = 2, shape2 = 5)
x <- cbind(x1 = x1, x2 = stats::rbeta(n, 2 + 6 * x1, 4))
set.seed(123)
fit <- fit_boostpm(
  x,
  Omega = cbind(lower = c(0, 0), upper = c(1, 1)),
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  add_noise = FALSE
)
print(summary(fit))

Simulate from a Fitted boostPM Distribution

Description

Draws samples from the probability distribution represented by a fitted tree ensemble.

Usage

## S3 method for class 'boostPM_fit'
simulate(object, nsim = 1L, seed = NULL, ...)

Arguments

Details

With seed = NULL, randomness is controlled by calling set.seed() before simulate(). Supplying seed sets R's random-number generator and changes its subsequent state.

Value

A numeric matrix with nsim rows and one column per fitted variable, in the same order as the training data. Training variable names are preserved as column names when available. Every draw lies within the fitted rectangular support.

Examples

set.seed(42)
n <- 400L
x1 <- stats::rbeta(n, shape1 = 2, shape2 = 5)
x <- cbind(x1 = x1, x2 = stats::rbeta(n, 2 + 6 * x1, 4))
set.seed(123)
fit <- fit_boostpm(
  x,
  Omega = cbind(lower = c(0, 0), upper = c(1, 1)),
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  add_noise = FALSE
)
generated <- simulate(fit, nsim = 500, seed = 321)

old_par <- graphics::par(mfrow = c(1, 2))
graphics::plot(x, xlab = "x1", ylab = "x2", main = "Training data")
graphics::plot(
  generated,
  xlab = "x1",
  ylab = "x2",
  main = "Generated data"
)
graphics::par(old_par)

Summarize a Fitted boostPM Distribution

Description

Extracts compact structural diagnostics from a boostPM_fit object.

Usage

## S3 method for class 'boostPM_fit'
summary(object, ...)

Arguments

Value

An object of class summary.boostPM_fit, represented as a list containing:

call

The matched call used to fit the model.

n_observations

Number of training observations, or NA when it cannot be recovered from the fitted object.

n_variables

Number of fitted variables.

n_trees

Number of accepted trees.

support

The original-scale lower and upper support limits.

tree_diagnostics

The per-tree diagnostic data frame returned by fit_boostpm().

variable_importance

Accumulated improvement attributed to each variable.

elapsed_time

Elapsed fitting time as a difftime object.

Examples

set.seed(42)
n <- 400L
x1 <- stats::rbeta(n, shape1 = 2, shape2 = 5)
x <- cbind(x1 = x1, x2 = stats::rbeta(n, 2 + 6 * x1, 4))
set.seed(123)
fit <- fit_boostpm(
  x,
  Omega = cbind(lower = c(0, 0), upper = c(1, 1)),
  max_marginal_trees = 100,
  max_dependence_trees = 1000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  add_noise = FALSE
)
summary(fit)