reworke twCoefLogitnormMLEFlat, mirgrate to testthat

DESCRIPTION

@@ -4,9 +4,13 @@ Version: 0.8.38
 Author: Thomas Wutzler
 Maintainer: Thomas Wutzler <[email protected]>
 Description: Density, distribution, quantile and random generation function for the logitnormal distribution. Estimation of the mode and the first two moments. Estimation of distribution parameters.
-Depends: 
-Suggests: RUnit, knitr, markdown, ggplot2, reshape2
+Suggests: 
+    knitr,
+    markdown,
+    ggplot2,
+    testthat (>= 3.0.0)
 VignetteBuilder: knitr
 License: GPL-2
 LazyData: true
 RoxygenNote: 7.1.1
+Config/testthat/edition: 3

NEWS.md

@@ -1,3 +1,7 @@
+# logitnorm 0.8.38
+- Reworked twCoefLogitnormMLEFlat
+- Migrated from RUnit to testthat
+
 # logitnorm 0.8.37
 
 - Density returns 0 instead of NAN outside (0,1)

R/logitnorm.R

@@ -97,6 +97,8 @@ qlogitnorm <- function(
   plogis(qn)
 }
 
+# derivative of logit(x)
+.dlogit <- function(x) 1/x + 1/(1-x)
 
 #------------------ estimating parameters from fitting to distribution 
 .ofLogitnorm <- function(
@@ -307,8 +309,6 @@ twCoefLogitnormMLE <- function(
 attr(twCoefLogitnormMLE,"ex") <- function(){
   # estimate the parameters, with mode 0.7 and upper quantile 0.9
   mode = 0.7; upper = 0.9
-  mode = 0.2; upper = 0.7
-  #mode = 0.5; upper = 0.9
   (theta <- twCoefLogitnormMLE(mode,upper))
   x <- seq(0,1,length.out = 41)[-c(1,41)]	# plotting grid
   px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function
@@ -321,18 +321,56 @@ attr(twCoefLogitnormMLE,"ex") <- function(){
   (theta <- twCoefLogitnormMLEFlat(mode))
 }
 
+.tmp.f <- function(){
+  (theta = twCoefLogitnormMLEFlat(0.9))
+  (theta = twCoefLogitnormMLEFlat(0.1))
+  (theta = twCoefLogitnormMLEFlat(0.5))
+  mle = c(0.7,0.2)
+  mle = seq(0.55,0.95,by=0.05)
+  mle = seq(0.95,0.99,by=0.01)
+  theta = twCoefLogitnormMLEFlat(mle)
+  sigma = theta[,2]
+  plot(sigma ~ mle)
+  (sigmac = 1/(2*mle*(1-mle)))
+  sigmac/sigma
+  plot(sigmac ~ mle)
+}
+
 twCoefLogitnormMLEFlat <- function(
-  ### Estimating coefficients of a maximally flat unimodal logitnormal distribution from mode 	
-  mle						##<< numeric vector: the mode of the density function
+  ### Estimating coefficients of a maximally flat unimodal logitnormal distribution given the mode 	
+  mle		##<< numeric vector: the mode of the density function
 ){
-  ##details<<
-  ## When increasing the sigma parameter, the distribution becomes
-  ## eventually becomes bi-model, i.e. has two maxima.	
-  ## This function estimates parameters for given mode, so that the distribution assigns high  
-  ## density to a maximum range, i.e. is maximally flat, but still is unimodal.
-  twCoefLogitnormMLE(mle,0)
+  vectorList <- lapply(mle, .twCoefLogitnormMLEFlat_scalar)
+  do.call(rbind, vectorList)
+}
+
+
+.twCoefLogitnormMLEFlat_scalar <- function(mle){
+  if (mle == 0.5) {
+    #sigma =  sqrt(dlogit(mle)/2) # 1.414214
+    return(c(mu = 0.0, sigma = sqrt(2)))
+  } 
+  is_right <- mle > 0.5
+  mler = ifelse(is_right, mle, 1-mle)
+  xt <- optimize(.ofModeFlat, interval = c(0,0.5), m=mler, logitm=logit(mler))$minimum
+  sigma2 <- (1/xt + 1/(1-xt))/2
+  mur <- logit(mler) - sigma2*(2*mler - 1)
+  mu <- ifelse(is_right, mur, -mur)
+  c(mu = mu, sigma = sqrt(sigma2))
 }
 
+.ofModeFlat <- function(x, m, logitm = logit(m)){
+  # lhs = 1/x + 1/(1-x)
+  # rhs = (logitm - logit(x))/(m-x)
+  # lhs = (m-x)/x + (m-x)/(1-x)
+  # rhs = (logitm - logit(x))
+  lhs = m/x + (m-1)/(1-x)
+  rhs = logitm - logit(x)
+  d = lhs - rhs
+  d*d
+}
+
+
 .ofLogitnormE <- function(
   ### Objective function used by \code{\link{coefLogitnormE}}. 
   theta				##<< theta[1] is the mu, theta[2] is the sigma

man/twCoefLogitnormMLE.Rd

@@ -21,8 +21,6 @@ rows correspond to rows in \code{mle}, \code{quant}, and \code{perc}}
 \examples{
 # estimate the parameters, with mode 0.7 and upper quantile 0.9
 mode = 0.7; upper = 0.9
-mode = 0.2; upper = 0.7
-#mode = 0.5; upper = 0.9
 (theta <- twCoefLogitnormMLE(mode,upper))
 x <- seq(0,1,length.out = 41)[-c(1,41)]	# plotting grid
 px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function

man/twCoefLogitnormMLEFlat.Rd

@@ -1,15 +1,12 @@
 \name{twCoefLogitnormMLEFlat}
 \alias{twCoefLogitnormMLEFlat}
 \title{twCoefLogitnormMLEFlat}
-\description{Estimating coefficients of a maximally flat unimodal logitnormal distribution from mode 	}
+\description{Estimating coefficients of a maximally flat unimodal logitnormal distribution given the mode 	}
 \usage{twCoefLogitnormMLEFlat(mle)}
 \arguments{
   \item{mle}{numeric vector: the mode of the density function}
 }
-\details{When increasing the sigma parameter, the distribution becomes
-eventually becomes bi-model, i.e. has two maxima.	
-This function estimates parameters for given mode, so that the distribution assigns high  
-density to a maximum range, i.e. is maximally flat, but still is unimodal.}
+
 
 
 \author{Thomas Wutzler}

tests/testthat.R

@@ -0,0 +1,4 @@
+library(testthat)
+library(logitnorm)
+
+test_check("logitnorm")

inst/unitTests/runitdlogitnormLog.R → tests/testthat/test_dlogitnormLog.R

@@ -1,22 +1,22 @@
-#require(logitnorm)
-
-test.dlogitnorm <- function() {
-	# Check dlogitnorm with arbitrary parameters at arbitrary points.
-	m <- 1.4
-	s <- 0.8
-	x <- c(0.1, 0.8);
-	# taken from https://en.wikipedia.org/w/index.php?title=Logit-normal_distribution&oldid=708557844
-	x.density <- 1/s/sqrt(2*pi) * exp(-(log(x/(1-x)) - m)^2/2/s^2) / x / (1-x)
-	stopifnot(all.equal(x.density, dlogitnorm(x, m, s)))
-	stopifnot(all.equal(log(x.density), dlogitnorm(x, m, s, log=TRUE)))
-}
-
-test.plogitnorm <- function() {
-	# Check plogitnorm with arbitrary parameters at arbitrary points.
-	m <- 1.4
-	s <- 0.8
-	x <- c(0.1, 0.8);
-	stopifnot(all.equal(pnorm(log(x/(1-x)), m, s), plogitnorm(x, m, s)))
-	stopifnot(all.equal(pnorm(log(x/(1-x)), m, s, log=TRUE), log(plogitnorm(x, m, s))))
-}
-
+#require(logitnorm)
+
+test_that("dlogitnorm",  {
+	# Check dlogitnorm with arbitrary parameters at arbitrary points.
+	m <- 1.4
+	s <- 0.8
+	x <- c(0.1, 0.8);
+	# taken from https://en.wikipedia.org/w/index.php?title=Logit-normal_distribution&oldid=708557844
+	x.density <- 1/s/sqrt(2*pi) * exp(-(log(x/(1-x)) - m)^2/2/s^2) / x / (1-x)
+	expect_equal(x.density, dlogitnorm(x, m, s))
+	expect_equal(log(x.density), dlogitnorm(x, m, s, log=TRUE))
+})
+
+test_that("plogitnorm",  {
+	# Check plogitnorm with arbitrary parameters at arbitrary points.
+	m <- 1.4
+	s <- 0.8
+	x <- c(0.1, 0.8);
+	expect_equal(pnorm(log(x/(1-x)), m, s), plogitnorm(x, m, s))
+	expect_equal(pnorm(log(x/(1-x)), m, s, log=TRUE), log(plogitnorm(x, m, s)))
+})
+