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

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

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

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

## ----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------------------------------------------------------------------
print(fit_beta)
summary(fit_sinusoidal)

## ----variable-importance, fig.width=6.5, fig.height=4.5-----------------------
plot(fit_sinusoidal, type = "variable_importance")

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

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

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

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

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

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

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