Introduction to boostPM

The boostPM package estimates probability distributions with an ensemble of tree-based probability measures. The method is useful when a distribution contains nonlinear dependence, multimodality, or other structure that is difficult to describe with a single parametric family.

For each example, we fit an ensemble, evaluate its density on a grid, and draw new observations from the fitted distribution.

Simulated training data

Observations occupy rows and variables occupy columns. Both examples use 2,400 training observations and remain practical for routine package checks.

Scenario 1: independent beta variables

The first coordinate follows a beta distribution concentrated toward the lower boundary, while the second is concentrated toward the upper boundary. The coordinates are independent, giving a simple parametric reference case with an analytic joint density.

library(boostPM)

n_beta <- 2400L
n_sinusoidal <- 2400L
support <- cbind(c(0, 0), c(1, 1))

set.seed(11001)
beta_data <- cbind(
  x1 = rbeta(n_beta, shape1 = 2, shape2 = 5),
  x2 = rbeta(n_beta, shape1 = 5, shape2 = 2)
)

beta_density <- function(x1, x2) {
  dbeta(x1, shape1 = 2, shape2 = 5) *
    dbeta(x2, shape1 = 5, shape2 = 2)
}

Scenario 2: sinusoidal ridge

The first coordinate is uniform. The conditional mean of the second coordinate follows a sine wave. This produces nonlinear dependence that is not well represented by an elliptical distribution.

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

sinusoidal_density <- function(x1, x2) {
  mu <- 0.15 + 0.70 * (0.5 + 0.5 * sin(2 * pi * x1))
  dbeta(x2, mu * concentration, (1 - mu) * concentration)
}

old_par <- par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
draw_points(beta_data, "Independent beta variables")
draw_points(sinusoidal_data, "Sinusoidal ridge")

par(old_par)

Fit the ensembles

The tree-number, learning-rate, and stopping controls follow the accuracy-oriented experimental configuration. The marginal stage is capped at 100 trees per dimension, and the dependence stage is capped at 5,000 trees. A learning rate of 0.1 is combined with adaptive stopping based on the most recent 50 validation improvements. The split prior of 0.9 follows the archived implementation and public experiment code.

set.seed(31001)
fit_beta <- fit_boostpm(
  beta_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
)

set.seed(32001)
fit_sinusoidal <- fit_boostpm(
  sinusoidal_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
)

Inspect fitted objects

fit_boostpm() is silent by default. Set progress = "stage" when stage-level messages are useful. The fitted object supports standard S3 methods.

print(fit_beta)

#> boostPM fit
#>   Observations: 2400 
#>   Variables: 2 
#>   Trees: 217 
#>   Elapsed time: 0.188597 secs 
#>   Variable importance: x1 = 0.4482104, x2 = 0.4271163

summary(fit_sinusoidal)

#> Summary of boostPM fit
#> 
#> Call:
#> fit_boostpm(data = sinusoidal_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)
#> 
#> Fit dimensions:
#>   Observations: 2400 
#>   Variables: 2 
#>   Trees: 478 
#> 
#> Support:
#>    lower upper
#> x1     0     1
#> x2     0     1
#> 
#> Nodes per tree: min 1, median 1, max 373 
#> Maximum tree depth: min 0, median 0, max 16 
#>   Elapsed time: 0.728838 secs 
#>   Variable importance: x1 = 0.6305273, x2 = 0.6477278

Variable-importance, number-of-nodes, and maximum-depth plots are available through plot().

plot(fit_sinusoidal, type = "variable_importance")

Evaluate the fitted densities

predict() returns log density by default. Here both fitted densities are evaluated on an 80 by 80 midpoint grid.

grid_size <- 80L
x1_grid <- (seq_len(grid_size) - 0.5) / grid_size
x2_grid <- (seq_len(grid_size) - 0.5) / grid_size
evaluation_grid <- as.matrix(expand.grid(
  x1 = x1_grid,
  x2 = x2_grid
))

fitted_beta <- predict(
  fit_beta,
  evaluation_grid,
  type = "density"
)
fitted_sinusoidal <- predict(
  fit_sinusoidal,
  evaluation_grid,
  type = "density"
)

true_beta <- beta_density(
  evaluation_grid[, 1],
  evaluation_grid[, 2]
)
true_sinusoidal <- sinusoidal_density(
  evaluation_grid[, 1],
  evaluation_grid[, 2]
)

The true and fitted panels share a color scale within each scenario. Values above the scenario-specific 99th percentile are clipped for display only.

beta_cap <- as.numeric(stats::quantile(
  c(true_beta, fitted_beta),
  0.99,
  names = FALSE
))
sinusoidal_cap <- as.numeric(stats::quantile(
  c(true_sinusoidal, fitted_sinusoidal),
  0.99,
  names = FALSE
))

old_par <- par(
  mfrow = c(2, 3),
  mar = c(4.5, 4.2, 2.8, 0.8),
  mgp = c(2.6, 0.8, 0)
)
draw_points(beta_data, "Training: beta")
draw_density(
  x1_grid,
  x2_grid,
  true_beta,
  "True density",
  beta_cap
)
draw_density(
  x1_grid,
  x2_grid,
  fitted_beta,
  "Fitted density",
  beta_cap
)
draw_points(sinusoidal_data, "Training: sinusoidal")
draw_density(
  x1_grid,
  x2_grid,
  true_sinusoidal,
  "True density",
  sinusoidal_cap
)
draw_density(
  x1_grid,
  x2_grid,
  fitted_sinusoidal,
  "Fitted density",
  sinusoidal_cap
)

par(old_par)

The fitted densities visually capture the main features of both the beta-product density and the nonlinear sinusoidal ridge.

Points outside the support have log density -Inf and density zero.

outside <- rbind(c(-0.10, 0.50), c(0.50, 0.50), c(1.10, 0.50))
predict(fit_sinusoidal, outside, type = "log_density")

#> [1]     -Inf 1.485438     -Inf

Generate and compare samples

The fitted ensemble is also a generative model. The following calls draw 1,000 observations from each fit.

set.seed(41001)
generated_beta <- simulate(fit_beta, nsim = 1000)

set.seed(42001)
generated_sinusoidal <- simulate(fit_sinusoidal, nsim = 1000)

old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
draw_points(beta_data, "Training: beta")
draw_points(
  generated_beta,
  "Generated: beta",
  color = "#2166AC",
  cex = 0.35
)
draw_points(sinusoidal_data, "Training: sinusoidal")
draw_points(
  generated_sinusoidal,
  "Generated: sinusoidal",
  color = "#2166AC",
  cex = 0.35
)

par(old_par)

The generated observations reproduce the main visible features of the beta marginal shapes and nonlinear sinusoidal structure in the training samples. Marginal empirical distribution functions provide a second comparison.

data_pairs <- list(
  list("Independent beta", beta_data, generated_beta),
  list("Sinusoidal", sinusoidal_data, generated_sinusoidal)
)

old_par <- par(mfrow = c(2, 2), mar = c(4, 4, 3, 1))
for (data_pair in data_pairs) {
  for (j in seq_len(2L)) {
    plot(
      stats::ecdf(data_pair[[2]][, j]),
      do.points = FALSE,
      verticals = TRUE,
      xlim = support[j, ],
      lwd = 2,
      col = "black",
      xlab = "",
      ylab = "Empirical cumulative probability",
      main = paste(data_pair[[1]], paste0("x", j))
    )
    mtext(colnames(data_pair[[2]])[j], side = 1, line = 2.5)
    lines(
      stats::ecdf(data_pair[[3]][, j]),
      verticals = TRUE,
      lwd = 2,
      col = "#2166AC"
    )
    legend(
      "bottomright",
      legend = c("Training", "Generated"),
      col = c("black", "#2166AC"),
      lwd = 2,
      bty = "n"
    )
  }
}

par(old_par)

Practical notes

Use a finite numeric matrix. Constant columns are not supported.
Omega must have two columns and positive width in each dimension.
n_bins controls the uniform candidate split grid and must be at least two.
c0 is strictly between zero and one; gamma is non-negative.
More trees and deeper weak learners increase flexibility and computation.