boostPM

boostPM fits multivariate probability distributions using unsupervised tree boosting. It supports model fitting, analytic density evaluation, random sampling, and diagnostic summaries.

The stable release is available from CRAN. This repository contains the development version of the package. The immutable research implementation used for the published work is retained under original/, with documented regression comparisons against the current implementation.

Installation

Install the stable release from CRAN:

install.packages("boostPM")
library(boostPM)

To install the development version from GitHub:

# install.packages("remotes")
remotes::install_github("nawaya040/boostPM-cran")

Quick start

The following example fits a nonlinear two-dimensional distribution. The first coordinate is uniform, while the conditional mean of the second coordinate follows a sinusoidal curve.

library(boostPM)

n <- 2400L

set.seed(12001)
x1 <- runif(n)
mu <- 0.15 + 0.70 * (0.5 + 0.5 * sin(2 * pi * x1))
concentration <- 40
x2 <- rbeta(
  n,
  shape1 = mu * concentration,
  shape2 = (1 - mu) * concentration
)

training_data <- cbind(x1 = x1, x2 = x2)
support <- cbind(
  lower = c(0, 0),
  upper = c(1, 1)
)

set.seed(32001)
fit <- fit_boostpm(
  training_data,
  Omega = support,
  max_marginal_trees = 100,
  max_dependence_trees = 5000,
  n_bins = 100,
  max_split_depth = 15,
  min_node_observations = 10,
  c0 = 0.1,
  gamma = 0.5,
  early_stop = c(0, 50),
  prior_split_prob = 0.9,
  add_noise = FALSE,
  progress = "stage"
)

The tree counts are upper bounds. With early_stop = c(0, 50), a fitting stage can stop before reaching its maximum when the mean of recent held-out log-density improvements falls below zero.

Evaluate the fitted density on a regular grid:

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

Draw a generated sample with the same size as the training data:

simulated <- simulate(fit, nsim = n, seed = 42001)

Compare the training data, fitted density, and generated data:

point_color <- grDevices::adjustcolor("#17324D", alpha.f = 0.3)
old_par <- par(mfrow = c(1, 3), mar = c(4.2, 4.2, 2.8, 1.2))

plot(
  training_data,
  xlim = c(0, 1),
  ylim = c(0, 1),
  xaxs = "i",
  yaxs = "i",
  asp = 1,
  pch = 16,
  cex = 0.35,
  col = point_color,
  xlab = "x1",
  ylab = "x2",
  main = "Training data"
)

image(
  grid,
  grid,
  fitted_density,
  xlim = c(0, 1),
  ylim = c(0, 1),
  xaxs = "i",
  yaxs = "i",
  asp = 1,
  col = grDevices::hcl.colors(64, "YlOrRd", rev = TRUE),
  xlab = "x1",
  ylab = "x2",
  main = "Fitted density"
)
contour(
  grid,
  grid,
  fitted_density,
  add = TRUE,
  drawlabels = FALSE,
  col = grDevices::adjustcolor("black", alpha.f = 0.55),
  nlevels = 8
)

plot(
  simulated,
  xlim = c(0, 1),
  ylim = c(0, 1),
  xaxs = "i",
  yaxs = "i",
  asp = 1,
  pch = 16,
  cex = 0.35,
  col = point_color,
  xlab = "x1",
  ylab = "x2",
  main = "Generated data"
)

par(old_par)

Training data, fitted density, and generated data

This example uses the accuracy-oriented configuration described in the introductory vignette.

Inspect a fitted model

print() provides a compact overview of the fitted model, including the number of observations, variables, accepted trees, elapsed time, and variable importance.

print(fit)

summary() additionally reports the fitted call, support limits, and ranges of tree sizes and depths.

fit_summary <- summary(fit)
print(fit_summary)

The fitted object provides three diagnostic plots:

old_par <- par(mfrow = c(1, 3))

plot(fit, type = "variable_importance")
plot(fit, type = "tree_node_counts")
plot(fit, type = "tree_depths")

par(old_par)

When adaptive stopping is enabled, held-out improvements and acceptance decisions are stored in heldout_diagnostics.

tail(fit$heldout_diagnostics)

Variable importance is an unnormalized within-model diagnostic. Larger values indicate greater accumulated splitting improvement, but the values should not be interpreted as probabilities.

Input and reproducibility notes

data must be a finite numeric matrix with observations in rows and variables in columns. Missing values, infinite values, and constant columns are not supported.
A supplied Omega must contain one lower and one upper support limit for each variable. All training observations must lie strictly inside those limits. If Omega = NULL, the processed range of each variable is expanded by 10 percent on both sides.
Density evaluation outside the fitted support returns zero. Equivalently, predict(fit, newdata, type = "log_density") returns -Inf outside the support. Samples returned by simulate() remain inside the fitted support.
max_marginal_trees and max_dependence_trees are upper bounds. When early_stop is supplied, fitting can stop before either bound is reached. The setting early_stop = c(0, 50) uses the paper-oriented held-out stopping rule.
Larger values of n_bins provide a finer search over candidate split locations, at greater computational cost.
With add_noise = TRUE, tied observations are randomly perturbed before fitting. Call set.seed() before fit_boostpm() when reproducibility is required. For continuously simulated data without ties, add_noise = FALSE avoids unnecessary perturbation.
Fitting and simulation use R's random-number generator. The seed argument to simulate() provides a reproducible generated sample and changes the subsequent RNG state.
progress = "stage" reports the marginal and dependence fitting stages. It does not change the fitted numerical result.

The introductory vignette is available after installation:

vignette("boostPM-introduction", package = "boostPM")

Citation

When using boostPM, please cite the published method paper:

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

@article{AwayaMa2024,
  author  = {Awaya, Naoki and Ma, Li},
  title   = {Unsupervised Tree Boosting for Learning Probability Distributions},
  journal = {Journal of Machine Learning Research},
  year    = {2024},
  volume  = {25},
  number  = {198},
  pages   = {1--52},
  url     = {https://www.jmlr.org/papers/v25/22-0980.html}
}

For the installed-package citation record, use:

citation("boostPM")

Development and validation

Routine checks run on Windows, macOS, and Ubuntu through GitHub Actions. Original-versus-package regression and longer statistical validation runners live under validation/. Their results and performance records are documented under docs/ and benchmarks/.