Test: mc

Ivan Jacob Agaloos Pesigan

2020-07-16

library(testthat)
library(jeksterslabRboot)
context("Test mc.")

Length 1

Parameters

R <- 20000L
thetahat <- 100
vcovhat <- 1.5

Results

mc_star_length_1 <- mc(
  thetahat = thetahat,
  vcovhat = vcovhat,
  R = R
)
str(mc_star_length_1)
#>  num [1:20000] 101 100 100 100 101 ...
hist(
  mc_star_length_1,
  main = expression(
    paste(
      "Histogram of ",
      hat(theta),
      "*"
    )
  ),
  xlab = expression(
    paste(
      hat(theta),
      "*"
    )
  )
)

qqnorm(mc_star_length_1)
qqline(mc_star_length_1)

Length Greater than 1

Parameters

alphahat <- 0.3386
betahat <- 0.4510
alphahat_betahat <- alphahat * betahat
alphahat_betahat_ci_2.5 <- 0.0033
alphahat_betahat_ci_97.5 <- 0.2979
varhat_alphahat <- 0.1224^2
varhat_betahat <- 0.1460^2
thetahat <- c(
  alphahat,
  betahat
)
vcovhat <- matrix(
  data = c(
    varhat_alphahat,
    0.00,
    0.00,
    varhat_betahat
  ),
  ncol = 2
)

Results

mc_star_length_2 <- mc(
  thetahat = thetahat,
  vcovhat = vcovhat,
  R = R
)
str(mc_star_length_2)
#>  num [1:20000, 1:2] 0.0909 0.1302 0.4743 0.3352 0.2965 ...
#>  - attr(*, "dimnames")=List of 2
#>   ..$ : NULL
#>   ..$ : NULL
alphahat_betahat_star <- mc_star_length_2[, 1] * mc_star_length_2[, 2]
hist(
  alphahat_betahat_star,
  main = expression(
    paste(
      "Histogram of ",
      hat(alpha),
      hat(beta),
      "*"
    )
  ),
  xlab = expression(
    paste(
      hat(alpha),
      hat(beta),
      "*"
    )
  )
)

qqnorm(alphahat_betahat_star)
qqline(alphahat_betahat_star)

wald_out <- wald(
  thetahat = alphahat_betahat,
  sehat = sd(alphahat_betahat_star)
)
wald_out
#>    statistic            p           se      ci_0.05       ci_0.5       ci_2.5 
#>  1.997931065  0.045724134  0.076433368 -0.098797440 -0.044170709  0.002901952 
#>      ci_97.5      ci_99.5     ci_99.95 
#>  0.302515248  0.349587909  0.404214640

testthat

test_that("ci_2.5", {
  expect_equivalent(
    alphahat_betahat_ci_2.5,
    wald_out["ci_2.5"],
    tolerance = 0.05
  )
})

test_that("ci_97.5", {
  expect_equivalent(
    alphahat_betahat_ci_97.5,
    wald_out["ci_97.5"],
    tolerance = 0.05
  )
})