## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
old_options <- options(digits = 3, scipen = 999)
print_num <- function(x) print(sprintf("%.3f", x))

## ----setup, message = FALSE, warning = FALSE----------------------------------
library(algebraic.dist)

## ----data-gen-----------------------------------------------------------------
# we define the parameters of the MVN
M <- mvn(mu = 1:5, sigma = diag(1:5))
print(M)

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

## ----data-gen2----------------------------------------------------------------
# we generate a sample of size n
n <- 10000
# we sample from the MVN
data <- sampler(M)(n)

## -----------------------------------------------------------------------------
head(data, n=6)

## ----emp-construction---------------------------------------------------------
emp <- empirical_dist(data)
print(emp)

## -----------------------------------------------------------------------------
params(emp)

## ----support------------------------------------------------------------------
# generate a data frame with the dimension, supremum,
# and infimum of the MVN and empirical distribution
df <- data.frame(supremum.mvn = supremum(sup(M)),
                 supremum.emp = supremum(sup(emp)),
                 infimum.mvn = infimum(sup(M)),
                 infimum.emp = infimum(sup(emp)))
row.names(df) <- paste0("param",1:dim(sup(M)))
print(df)

## ----compare-basic-stats------------------------------------------------------
mean(M); mean(emp)

## ----compare-basic-stats2-----------------------------------------------------
vcov(M); vcov(emp)

## -----------------------------------------------------------------------------
diag(vcov(M)); diag(vcov(emp))

## -----------------------------------------------------------------------------
mu.emp <- mean(emp)
expectation(emp, function(x) (x - mu.emp)^2)
expectation(M, function(x) (x - mean(M))^2, control = list(n = n))

## -----------------------------------------------------------------------------
expectation(emp, function(x) (x - mu.emp)^2, control = list(compute_stats = TRUE))
expectation(M, function(x) (x - mean(M))^2, control = list(n = n, compute_stats = TRUE))

## ----rmap---------------------------------------------------------------------
mu.emp <- mean(emp)
mean(rmap(emp, function(x) (x - mu.emp)^2))
mean(rmap(M, function(x) (x - mean(M))^2, n = n))

## -----------------------------------------------------------------------------
x <- sampler(emp)(1)
y <- sampler(M)(1)
data.frame(
  ob = c("empirical(x)", "MVN(y)"),
  pdf.emp = c(density(emp)(x), density(emp)(y)),
  pdf.mvn = c(density(M)(x), density(M)(y)))

## ----marginal-test, fig.height = 4, fig.width = 6-----------------------------
X1.emp <- marginal(emp, 1)
F1.emp <- cdf(X1.emp)
curve(F1.emp(x), from = infimum(sup(X1.emp)), to = supremum(sup(X1.emp)), col = "blue", lty = 1)
X1 <- marginal(M, 1)
F1 <- cdf(X1)
curve(F1(x), from = infimum(sup(X1.emp)), to = supremum(sup(X1.emp)), add = TRUE, col = "red", lty = 2)
legend("topleft", legend = c("empirical", "MVN"), col = c("blue", "red"), lty = c(1, 2))

## ----marginal-test-zoom, fig.height = 4, fig.width = 6------------------------
curve(F1.emp(x), from = 1, to = 1.0005, col = "blue", lty = 1, type = "s")

## ----expectation-test---------------------------------------------------------
cat("E(X1) = ", expectation(X1, function(x) x), "( expected ", mean(X1), ")\n",
    "Var(X1) = ", expectation(X1,
                              function(x) {
                                  (x - expectation(X1,
                                                   function(x) x)
                                  )^2 
                              }),
    "( expected ", vcov(X1), ")\n",
    "Skewness(X1) = ", expectation(X1,
                                   function(x) {
                                       (x - expectation(X1, function(x) x))^3 / 
                                       expectation(X1, function(x) x)^3 
                                   }),
    "( expected 0 )\n",
    "E(X^2) = ", expectation(X1, function(x) x^2),
    "( expected ", vcov(X1) + mean(X1)^2, ")")

## ----joint-marginal-test, fig.height = 4, fig.width = 6-----------------------
dataX2X4.emp <- sampler(marginal(emp, c(2, 4)))(10000)
dataX2X4.mvn <- sampler(marginal(M, c(2, 4)))(10000)

# scatter plot a 2d sample. use xlab and ylab to label the axes
plot(dataX2X4.emp[,1], dataX2X4.emp[,2], pch = 20, col = "blue", xlab = "X2", ylab = "X4")
# overlay in green the MVN data
points(dataX2X4.mvn[,1], dataX2X4.mvn[,2], pch = 20, col = "green")
legend("topright", legend = c("empirical", "MVN"), col = c("blue", "green"), pch = c(20, 20))
title("Scatter plot of X2 and X4")

## ----multivariate-cdf---------------------------------------------------------
cdf(M)(c(1,2,3,4,0))

## ----conditional-test---------------------------------------------------------
mean(conditional(emp, function(x) x[1] + x[2] < 0))
mean(conditional(M, function(x) x[1] + x[2] < 0))

## -----------------------------------------------------------------------------
e <- edist(expression(x + y),
           list("x" = normal(mu = 0, var = 1),
                "y" = exponential(rate = 1)))

## -----------------------------------------------------------------------------
print(e)

## -----------------------------------------------------------------------------
mean(e)
vcov(e)

## -----------------------------------------------------------------------------
params(e)

## -----------------------------------------------------------------------------
s <- sampler(e)
samp <- s(10)  # Generate a sample of size 10
print(samp)

## ----include = FALSE----------------------------------------------------------
options(old_options)