## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----library------------------------------------------------------------------
library(algebraic.dist)

## ----mvn-construct------------------------------------------------------------
M <- mvn(mu = c(1, 2, 3),
         sigma = matrix(c(1.0, 0.5, 0.2,
                          0.5, 2.0, 0.3,
                          0.2, 0.3, 1.5), nrow = 3, byrow = TRUE))
M

## ----mvn-accessors------------------------------------------------------------
mean(M)
vcov(M)
dim(M)

## ----mvn-params---------------------------------------------------------------
params(M)

## ----mvn-sampling-------------------------------------------------------------
set.seed(42)
samp <- sampler(M)
X <- samp(100)
dim(X)
head(X, 5)

## ----mvn-marginal-------------------------------------------------------------
M13 <- marginal(M, c(1, 3))
M13
mean(M13)
vcov(M13)

## ----mvn-marginal-univariate--------------------------------------------------
M1 <- marginal(M, 1)
M1
mean(M1)
vcov(M1)

## ----mvn-conditional-closedform-----------------------------------------------
M2 <- mvn(mu = c(0, 0), sigma = matrix(c(1, 0.8, 0.8, 1), 2, 2))

# Condition on X2 = 1
M_cond <- conditional(M2, given_indices = 2, given_values = 1)
mean(M_cond)
vcov(M_cond)

## ----mvn-conditional-mc-------------------------------------------------------
set.seed(42)
M_cond_mc <- conditional(M2, P = function(x) abs(x[2] - 1) < 0.1)
mean(M_cond_mc)

## ----affine-portfolio---------------------------------------------------------
returns <- mvn(
  mu = c(0.05, 0.08),
  sigma = matrix(c(0.04, 0.01,
                   0.01, 0.09), 2, 2)
)

# 60/40 portfolio weights
w <- matrix(c(0.6, 0.4), nrow = 1)
portfolio <- affine_transform(returns, A = w)
mean(portfolio)
vcov(portfolio)

## ----affine-sumdiff-----------------------------------------------------------
X <- mvn(mu = c(3, 7), sigma = matrix(c(4, 1, 1, 9), 2, 2))

# A maps (X1, X2) -> (X1 + X2, X1 - X2)
A <- matrix(c(1, 1,
              1, -1), nrow = 2, byrow = TRUE)
Y <- affine_transform(X, A = A)
mean(Y)
vcov(Y)

## ----mixture-construct--------------------------------------------------------
m <- mixture(
  components = list(normal(-2, 1), normal(3, 0.5)),
  weights = c(0.4, 0.6)
)
m

## ----mixture-mean-------------------------------------------------------------
mean(m)

## ----mixture-vcov-------------------------------------------------------------
vcov(m)

## ----mixture-density, fig.height = 4, fig.width = 6---------------------------
f <- density(m)

x_vals <- seq(-6, 8, length.out = 200)
y_vals <- sapply(x_vals, f)
plot(x_vals, y_vals, type = "l", lwd = 2,
     main = "Bimodal mixture density",
     xlab = "x", ylab = "f(x)")

## ----mixture-cdf, fig.height = 4, fig.width = 6-------------------------------
F <- cdf(m)
F_vals <- sapply(x_vals, F)
plot(x_vals, F_vals, type = "l", lwd = 2,
     main = "Mixture CDF",
     xlab = "x", ylab = "F(x)")

## ----mixture-sampling, fig.height = 4, fig.width = 6--------------------------
set.seed(123)
s <- sampler(m)
samples <- s(2000)
hist(samples, breaks = 50, probability = TRUE, col = "lightblue",
     main = "Mixture samples with true density",
     xlab = "x")
lines(x_vals, y_vals, col = "red", lwd = 2)

## ----mixture-marginal---------------------------------------------------------
m2 <- mixture(
  list(mvn(c(0, 0), diag(2)), mvn(c(3, 3), diag(2))),
  c(0.5, 0.5)
)

m2_x1 <- marginal(m2, 1)
m2_x1
mean(m2_x1)

## ----mixture-conditional------------------------------------------------------
mc <- conditional(m2, given_indices = 2, given_values = 0)
mc
mean(mc)

## ----mixture-conditional-high-------------------------------------------------
mc_high <- conditional(m2, given_indices = 2, given_values = 3)
mc_high
mean(mc_high)

## ----gmm-setup----------------------------------------------------------------
class_a <- mvn(c(-2, -1), matrix(c(1.0, 0.3, 0.3, 0.8), 2, 2))
class_b <- mvn(c(2, 1.5), matrix(c(0.6, -0.2, -0.2, 1.2), 2, 2))

gmm <- mixture(
  components = list(class_a, class_b),
  weights = c(0.5, 0.5)
)
gmm

## ----gmm-scatter, fig.height = 5, fig.width = 6-------------------------------
set.seed(314)
samp_a <- sampler(class_a)(300)
samp_b <- sampler(class_b)(300)

plot(samp_a[, 1], samp_a[, 2], col = "steelblue", pch = 16, cex = 0.6,
     xlim = c(-6, 6), ylim = c(-5, 5),
     xlab = expression(X[1]), ylab = expression(X[2]),
     main = "Two-class GMM samples")
points(samp_b[, 1], samp_b[, 2], col = "tomato", pch = 16, cex = 0.6)
legend("topright", legend = c("Class A", "Class B"),
       col = c("steelblue", "tomato"), pch = 16)

## ----gmm-classify-neg---------------------------------------------------------
post_neg <- conditional(gmm, given_indices = 1, given_values = -1)
post_neg

## ----gmm-classify-neg-density, fig.height = 4, fig.width = 6------------------
f_post <- density(post_neg)
x2_grid <- seq(-5, 5, length.out = 200)
y_post <- sapply(x2_grid, f_post)
plot(x2_grid, y_post, type = "l", lwd = 2,
     xlab = expression(X[2]),
     ylab = expression(f(X[2] ~ "|" ~ X[1] == -1)),
     main = "Posterior density of X2 given X1 = -1")

## ----gmm-classify-pos---------------------------------------------------------
post_pos <- conditional(gmm, given_indices = 1, given_values = 2)
post_pos

## ----gmm-classify-pos-density, fig.height = 4, fig.width = 6------------------
f_post_pos <- density(post_pos)
y_post_pos <- sapply(x2_grid, f_post_pos)
plot(x2_grid, y_post_pos, type = "l", lwd = 2,
     xlab = expression(X[2]),
     ylab = expression(f(X[2] ~ "|" ~ X[1] == 2)),
     main = "Posterior density of X2 given X1 = 2")

## ----gmm-classify-mid---------------------------------------------------------
post_mid <- conditional(gmm, given_indices = 1, given_values = 0)
post_mid

## ----gmm-classify-mid-density, fig.height = 4, fig.width = 6------------------
f_post_mid <- density(post_mid)
y_post_mid <- sapply(x2_grid, f_post_mid)
plot(x2_grid, y_post_mid, type = "l", lwd = 2,
     xlab = expression(X[2]),
     ylab = expression(f(X[2] ~ "|" ~ X[1] == 0)),
     main = "Posterior density of X2 given X1 = 0 (ambiguous)")