inst/doc/joker.R

In joker: Probability Distributions and Parameter Estimation

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

## -----------------------------------------------------------------------------
library(knitr)

kable(
  data.frame(
    Distribution = c("Bernoulli", "Beta", "Binomial", "Categorical", "Cauchy",
                     "Chi-Square", "Dirichlet", "Fisher", "Gamma", "Geometric"),
    Class_Name = c("Bern", "Beta", "Binom", "Cat", "Cauchy", "Chisq", "Dir",
                   "Fisher", "Gam", "Geom"),
    Distribution2 = c("Laplace", "Log-Normal", "Multivariate Gamma",
                      "Multinomial", "Negative Binomial", "Normal", "Poisson",
                      "Student", "Uniform", "Weibull"),
    Class_Name2 = c("Laplace", "Lnorm", "Multigam", "Multinom", "Nbinom",
                    "Norm", "Pois", "Stud", "Unif", "Weib")
  ),
  col.names = c("Distribution", "Class Name", "Distribution", "Class Name"),
  caption = "Overview of the distributions implemented in the `joker` package, 
  along with their respective class names."
)

## -----------------------------------------------------------------------------
shape1 <- 1 
shape2 <- 2
D <- Beta(shape1, shape2)

## -----------------------------------------------------------------------------
d(D, 0.5)
dbeta(0.5, shape1, shape2)

p(D, 0.5)
pbeta(0.5, shape1, shape2)

qn(D, 0.75)
qbeta(0.75, shape1, shape2)

r(D, 2)
rbeta(2, shape1, shape2)

## -----------------------------------------------------------------------------
F1 <- p(D)
F1(0.5)

## -----------------------------------------------------------------------------
mean(D)
median(D)
mode(D)
var(D)
sd(D)
skew(D)
kurt(D)
entro(D)
finf(D)

## -----------------------------------------------------------------------------
mode(Beta(1, 1))

## -----------------------------------------------------------------------------
set.seed(1)
shape1 <- 1
shape2 <- 2
D <- Beta(shape1, shape2)
x <- r(D, 100)

## -----------------------------------------------------------------------------
ebeta(x, type = "mle")
ebeta(x, type = "me")
ebeta(x, type = "same")

## -----------------------------------------------------------------------------
e(D, x, type = "mle")

## -----------------------------------------------------------------------------
mle(D, x)
me(D, x)
same(D, x)

## ----eval=FALSE---------------------------------------------------------------
# mle("beta", x)
# mle("bEtA", x)
# e("Beta", x, type = "mle")

## -----------------------------------------------------------------------------
llbeta(x, shape1, shape2)
ll(D, x)

## -----------------------------------------------------------------------------
vbeta(shape1, shape2, type = "mle")
vbeta(shape1, shape2, type = "me")
vbeta(shape1, shape2, type = "same")

## ----eval=FALSE---------------------------------------------------------------
# avar(D, type = "mle")
# avar_mle(D)
# avar_me(D)
# avar_same(D)

## -----------------------------------------------------------------------------
D <- Beta(1, 2)

prm <- list(name = "shape1",
            val = seq(1, 5, by = 0.5))

x <- small_metrics(D, prm,
             obs = c(20, 50),
             est = c("mle", "same", "me"),
             sam = 1e3,
             seed = 1)

head(x@df)

## ----echo=FALSE, fig.height=8, fig.width=14, out.width="100%", fig.cap="Small-sample metrics comparison for MLE, ME, and SAME of the beta distribution shape1 parameter."----
plot(x)

## -----------------------------------------------------------------------------
D <- Dir(alpha = 1:4)

prm <- list(name = "alpha",
            pos = 1,
            val = seq(1, 5, by = 0.5))

x <- small_metrics(D, prm,
                   obs = c(20, 50),
                   est = c("mle", "same", "me"),
                   sam = 1e3,
                   seed = 1)

class(x)
head(x@df)

## -----------------------------------------------------------------------------
D <- Beta(1, 2)

prm <- list(name = "shape1",
            val = seq(1, 5, by = 0.1))

x <- large_metrics(D, prm,
                   est = c("mle", "same", "me"))

class(x)
head(x@df)

## ----echo=FALSE, fig.height=8, fig.width=14,  out.width="100%", fig.cap="Large-sample metrics comparison for MLE, ME, and SAME of the beta distribution shape1 parameter."----
plot(x)

## ----eval = FALSE-------------------------------------------------------------
# setClass("Beta",
#   contains = "Distribution",
#   slots = c(shape1 = "numeric", shape2 = "numeric"),
#   prototype = list(shape1 = 1, shape2 = 1))

## ----eval = FALSE-------------------------------------------------------------
# Beta <- function(shape1 = 1, shape2 = 1) {
#   new("Beta", shape1 = shape1, shape2 = shape2)
# }
# 
# D <- Beta(1, 2)
# D@shape1
# D@shape2

## ----eval = FALSE-------------------------------------------------------------
# setValidity("Beta", function(object) {
#   if(length(object@shape1) != 1) {
#     stop("shape1 has to be a numeric of length 1")
#   }
#   if(object@shape1 <= 0) {
#     stop("shape1 has to be positive")
#   }
#   if(length(object@shape2) != 1) {
#     stop("shape2 has to be a numeric of length 1")
#   }
#   if(object@shape2 <= 0) {
#     stop("shape2 has to be positive")
#   }
#   TRUE
# })

## -----------------------------------------------------------------------------
# probability density function
setMethod("d", signature = c(distr = "Beta", x = "numeric"),
          function(distr, x) {
            dbeta(x, shape1 = distr@shape1, shape2 = distr@shape2)
          })

# (theoretical) expectation
setMethod("mean",
          signature  = c(x = "Beta"),
          definition = function(x) {

  x@shape1 / (x@shape1 + x@shape2)

})

# moment estimator
setMethod("me",
          signature  = c(distr = "Beta", x = "numeric"),
          definition = function(distr, x) {

  m  <- mean(x)
  m2 <- mean(x ^ 2)
  d  <- (m - m2) / (m2 - m ^ 2)

  c(shape1 = d * m, shape2 = d * (1 - m))

})