Test: Normal - Likelihood • jeksterslabRdist

library(testthat)
library(microbenchmark)
library(jeksterslabRdist)
context("Test Normal - Likelihood.")

Parameters

n <- 50
mu <- runif(
  n = 1,
  min = 1,
  max = 2
)
sigma <- runif(
  n = 1,
  min = 1,
  max = 2
)
Variable <- c(
  "`n`",
  "`mu`",
  "`sigma`"
)
Description <- c(
  "Sample size ($n$).",
  "Population mean ($\\mu$).",
  "Population standard deviation ($\\sigma$)."
)
Value <- c(
  n,
  mu,
  sigma
)
knitr::kable(
  x = data.frame(
    Variable,
    Description,
    Value
  ),
  row.names = FALSE
)

Variable	Description	Value
`n`	Sample size (\(n\)).	50.000000
`mu`	Population mean (\(\mu\)).	1.305921
`sigma`	Population standard deviation (\(\sigma\)).	1.449006

Generate Data

x <- rnorm(
  n = n,
  mean = mu,
  sd = sigma
)
x_bar <- mean(x)
s <- sd(x)
Variable <- c(
  "`n`",
  "`x_bar`",
  "`s`"
)
Description <- c(
  "Sample size ($n$).",
  "Sample mean ($\\bar{x}$).",
  "Sample standard deviation ($s$)."
)
Value <- c(
  n,
  x_bar,
  s
)
knitr::kable(
  x = data.frame(
    Variable,
    Description,
    Value
  ),
  row.names = FALSE
)

Variable	Description	Value
`n`	Sample size (\(n\)).	50.000000
`x_bar`	Sample mean (\(\bar{x}\)).	1.236637
`s`	Sample standard deviation (\(s\)).	1.698222

Likelihood

results_dnorm_L <- prod(
  dnorm(
    x = x,
    mean = mu,
    sd = sigma,
    log = FALSE
  )
)
results_dnorm_ll <- -1 * sum(
  dnorm(
    x = x,
    mean = mu,
    sd = sigma,
    log = TRUE
  )
)
results_dnorm_2ll <- -2 * sum(
  dnorm(
    x = x,
    mean = mu,
    sd = sigma,
    log = TRUE
  )
)
results_normL <- normL(
  mu = mu,
  sigma = sigma,
  x = x
)
results_normll <- normll(
  mu = mu,
  sigma = sigma,
  x = x,
  neg = TRUE
)
results_normll_neg_false <- normll(
  mu = mu,
  sigma = sigma,
  x = x,
  neg = FALSE
)
results_normll_L <- -1 * log(results_normL)
results_norm2ll <- norm2ll(
  mu = mu,
  sigma = sigma,
  x = x,
  neg = TRUE
)
results_norm2ll_neg_false <- norm2ll(
  mu = mu,
  sigma = sigma,
  x = x,
  neg = FALSE
)
results_norm2ll_L <- -2 * log(results_normL)

Maximum Likelihood Estimation

foo <- function(theta, x) {
  -2 * sum(
    dnorm(
      x = x,
      mean = theta[1],
      sd = theta[2],
      log = TRUE
    )
  )
}
results_ml_dnorm <- opt(
  FUN = foo,
  start_values = c(mu, sigma),
  optim = TRUE,
  x = x
)
results_ml_norm2ll <- opt(
  FUN = normobj,
  start_values = c(mu, sigma),
  optim = TRUE,
  x = x,
  neg = TRUE
)
results_nlminb_ml_dnorm <- opt(
  FUN = foo,
  start_values = c(mu, sigma),
  optim = FALSE,
  x = x
)
results_nlminb_ml_norm2ll <- opt(
  FUN = normobj,
  start_values = c(mu, sigma),
  optim = FALSE,
  x = x,
  neg = TRUE
)

Summarize Results

knitr::kable(
  x = data.frame(
    dnorm_ll = results_dnorm_ll,
    normll = results_normll,
    normL = results_normll_L
  ),
  row.names = FALSE,
  caption = "Negative Log-Likelihood"
)

Negative Log-Likelihood
dnorm_ll	normll	normL
98.20025	98.20025	98.20025

knitr::kable(
  x = data.frame(
    dnorm_2ll = results_dnorm_2ll,
    norm2ll = results_norm2ll,
    normL = results_norm2ll_L
  ),
  row.names = FALSE,
  caption = "Negative Two Log-Likelihood"
)

Negative Two Log-Likelihood
dnorm_2ll	norm2ll	normL
196.4005	196.4005	196.4005

knitr::kable(
  x = data.frame(
    Item = c("$\\mu$", "$\\sigma$"),
    Parameter = c(mu, sigma),
    Sample = c(x_bar, s),
    dnorm_optim = results_ml_dnorm$par,
    dnorm_nlminb = results_nlminb_ml_dnorm$par,
    normal_negll_optim = results_ml_norm2ll$par,
    normal_negll_nlminb = results_nlminb_ml_norm2ll$par
  ),
  row.names = FALSE,
  col.names = c(
    "Item",
    "Parameter",
    "Closed-form solution",
    "MLE using dnorm (optim)",
    "MLE using dnorm (nlminb)",
    "MLE using normal_negll (optim)",
    "MLE using normal_negll (nlminb)"
  ),
  caption = "Maximum Likelihood Estimates (MLE)"
)

Maximum Likelihood Estimates (MLE)
Item	Parameter	Closed-form solution	MLE using dnorm (optim)	MLE using dnorm (nlminb)	MLE using normal_negll (optim)	MLE using normal_negll (nlminb)
\(\mu\)	1.305921	1.236637	1.236839	1.236637	1.236839	1.236637
\(\sigma\)	1.449006	1.698222	1.681174	1.681154	1.681174	1.681154

Benchmarking

microbenchmark(
  dnorm_ll = -1 * sum(dnorm(x = x, mean = mu, sd = sigma, log = TRUE)),
  normll = normll(mu = mu, sigma = sigma, x = x, neg = TRUE),
  normll_neg_false = normll(mu = mu, sigma = sigma, x = x, neg = FALSE),
  dnorm_2ll = -2 * sum(dnorm(x = x, mean = mu, sd = sigma, log = TRUE)),
  norm2ll = norm2ll(mu = mu, sigma = sigma, x = x, neg = TRUE),
  norm2ll_neg_false = norm2ll(mu = mu, sigma = sigma, x = x, neg = FALSE)
)
#> Unit: microseconds
#>               expr   min     lq    mean median     uq    max neval
#>           dnorm_ll 1.629 1.7430 1.90184 1.8705 1.9815  2.765   100
#>             normll 1.043 1.1980 1.36585 1.3175 1.4650  2.963   100
#>   normll_neg_false 1.067 1.1880 1.50028 1.3425 1.4710 16.832   100
#>          dnorm_2ll 1.611 1.7655 2.30403 1.8960 2.0470 32.517   100
#>            norm2ll 1.054 1.1865 1.48943 1.3040 1.4855 12.895   100
#>  norm2ll_neg_false 1.065 1.1600 1.30007 1.2515 1.4100  2.060   100
microbenchmark(
  ml_dnorm_optim = opt(FUN = foo, start_values = c(mu, sigma), optim = TRUE, x = x),
  ml_dnorm_nlminb = opt(FUN = foo, start_values = c(mu, sigma), optim = FALSE, x = x),
  ml_normobj_optim = opt(FUN = normobj, start_values = c(mu, sigma), optim = TRUE, x = x, neg = TRUE),
  ml_normobj_nlminb = opt(FUN = normobj, start_values = c(mu, sigma), optim = FALSE, x = x, neg = TRUE)
)
#> Unit: microseconds
#>               expr     min       lq     mean   median       uq      max neval
#>     ml_dnorm_optim 129.422 144.3785 177.3860 169.5965 183.1555  534.880   100
#>    ml_dnorm_nlminb 127.350 154.1750 199.9750 162.7035 174.8530 2936.478   100
#>   ml_normobj_optim 118.128 149.8655 172.3927 163.0760 174.9120  345.270   100
#>  ml_normobj_nlminb 116.825 154.4895 177.6536 158.1360 179.7355  358.439   100

testthat

test_that("normL returns the same value as dnorm", {
  expect_equivalent(
    round(
      x = results_dnorm_L,
      digits = 2
    ),
    round(
      x = results_normL,
      digits = 2
    )
  )
})

test_that("normll returns the same value as dnorm", {
  expect_equivalent(
    round(
      x = results_dnorm_ll,
      digits = 2
    ),
    round(
      x = results_normll,
      digits = 2
    ),
    round(
      x = results_normll_neg_false * -1,
      digits = 2
    ),
    round(
      x = results_normll_L,
      digits = 2
    )
  )
})

test_that("norm2ll returns the same value as dnorm", {
  expect_equivalent(
    round(
      x = results_dnorm_2ll,
      digits = 2
    ),
    round(
      x = results_norm2ll,
      digits = 2
    ),
    round(
      x = results_norm2ll_neg_false * -1,
      digits = 2
    ),
    round(
      x = results_norm2ll_L,
      digits = 2
    )
  )
})

test_that("normobj maximizes to the same estimates as dnorm", {
  expect_equivalent(
    round(
      x = results_ml_dnorm$par,
      digits = 2
    ),
    round(
      x = results_nlminb_ml_dnorm$par,
      digits = 2
    ),
    round(
      x = results_ml_norm2ll$par,
      digits = 2
    ),
    round(
      x = results_nlminb_ml_norm2ll$par,
      digits = 2
    )
  )
})

test_that("NA output", {
  expect_true(
    is.na(normobj(
      theta = c(-2, -2),
      x = rnorm(10),
      neg = TRUE
    ))
  )
})

Test: Normal - Likelihood

Ivan Jacob Agaloos Pesigan

2020-08-06

Parameters

Generate Data

Likelihood

Maximum Likelihood Estimation

Summarize Results

Benchmarking

testthat