---
title: "Introduction to boostPM"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to boostPM}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  comment = "#>",
  fig.width = 7,
  fig.height = 4.8,
  fig.align = "center",
  warning = FALSE,
  message = FALSE
)

density_palette <- grDevices::hcl.colors(32, "YlOrRd", rev = TRUE)

draw_density <- function(x1, x2, values, title, display_cap) {
  surface <- matrix(values, nrow = length(x1), ncol = length(x2))
  surface <- pmin(surface, display_cap)
  image(
    x1,
    x2,
    surface,
    zlim = c(0, display_cap),
    col = density_palette,
    xlim = c(0, 1),
    ylim = c(0, 1),
    asp = 1,
    xlab = "x1",
    ylab = "x2",
    main = title
  )
  contour(
    x1,
    x2,
    surface,
    add = TRUE,
    drawlabels = FALSE,
    col = "grey30",
    nlevels = 7
  )
}

draw_points <- function(x, title, color = "black", cex = 0.35) {
  plot(
    x,
    xlim = c(0, 1),
    ylim = c(0, 1),
    xaxs = "i",
    yaxs = "i",
    xaxt = "n",
    yaxt = "n",
    asp = 1,
    pch = 16,
    cex = cex,
    col = grDevices::adjustcolor(color, alpha.f = 0.35),
    xlab = "x1",
    ylab = "x2",
    main = title
  )
  axis(1, at = seq(0, 1, by = 0.2))
  axis(2, at = seq(0, 1, by = 0.2))
}
```

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.

This vignette considers two distributions on the unit square:

1. a product of two beta distributions;
2. a sinusoidal conditional-density ridge.

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.

```{r beta-data}
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.

```{r sinusoidal-data}
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)
}
```

```{r training-data-plot, fig.width=7, fig.height=3.6}
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.

```{r fit-beta}
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
)
```

```{r fit-sinusoidal}
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.

```{r inspect}
print(fit_beta)
summary(fit_sinusoidal)
```

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

```{r variable-importance, fig.width=6.5, fig.height=4.5}
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.

```{r density-evaluation}
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.

```{r density-comparison, fig.width=9, fig.height=7}
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.

```{r outside-support}
outside <- rbind(c(-0.10, 0.50), c(0.50, 0.50), c(1.10, 0.50))
predict(fit_sinusoidal, outside, type = "log_density")
```

## Generate and compare samples

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

```{r simulate}
set.seed(41001)
generated_beta <- simulate(fit_beta, nsim = 1000)

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

```{r generated-joint-plot, fig.width=8, fig.height=7.2}
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.

```{r marginal-comparison, fig.width=8, fig.height=7.2}
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.

## References

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