In PejLab/bln: Binomial Logit Normal Distribution

knitr::opts_chunk$set(
  collapse = FALSE,
  comment = '#>',
  echo = TRUE,
  tidy = TRUE
)

Beta-Logit-Normal Distribution

# Load the bln package
library(bln)
# Find the included simluations conducted in MATLAB and read them in
blnfile <- system.file('extdata', 'matlab_simluations.txt', package = 'bln')
matlab <- read.table(
  file = blnfile,
  header = TRUE, 
  comment.char = '#', 
  as.is = TRUE
)

# Using the same x, size (x + xc), mean, and standard deviations
# run the approximate BLN PDF R code
system.time(
  expr = r.approx <- dbln(
    x = matlab$x,
    size = matlab$x + matlab$xc,
    mean = matlab$mu,
    sd = sqrt(x = matlab$v),
    approximate = TRUE
  )
)
plot(
  x = log(x = r.approx), 
  y = log(x = matlab$px), 
  main = cor(x = r.approx, y = matlab$px)
)

# Using the same x, size (x + xc), mean, and standard deviations,
# run the approximate BLN PDF C++ code
system.time(
  expr = cpp.approx <- dblnpp(
    x = matlab$x, 
    size = matlab$x + matlab$xc, 
    mean = matlab$mu, 
    sd = sqrt(x = matlab$v)
  )
)
plot(
  x = cpp.approx, 
  y = matlab$px, 
  main = cor(x = cpp.approx, y = matlab$px)
)

# Using the same x, size (x + xc), mean, and standard deviations
# run the exact BLN PDF R code
system.time(
  expr = r.exact <- dbln(
    x = matlab$x,
    size = matlab$x + matlab$xc,
    mean = matlab$mu,
    sd = sqrt(x = matlab$v),
    approximate = FALSE
  )
)
plot(
  x = log(x = r.exact), 
  y = log(x = matlab$px_acc), 
  main = cor(x = r.exact, y = matlab$px_acc)
)