From af6beb9112cdf60718d0861eacde3e4b2f1ff907 Mon Sep 17 00:00:00 2001 From: Thomas Wutzler Date: Fri, 8 Dec 2017 15:28:19 +0100 Subject: [PATCH] Avoid writing file reports during testing --- .Rbuildignore | 5 +- .gitignore | 12 +- .travis.yml | 16 +- DESCRIPTION | 2 +- NAMESPACE | 8 +- NEWS.md | 3 + R/logitnorm.R | 908 ++++++++++++++-------------- README.Rmd | 112 ++-- README.md | 112 ++-- cranComments.md | 15 +- inst/genData/_genPackage.R | 38 +- inst/genData/logitnorm-package.Rd | 110 ++-- inst/unitTests/runitLogitnorm.R | 298 ++++----- inst/unitTests/runitdlogitnormLog.R | 44 +- man/dlogitnorm.Rd | 58 +- man/logitnorm-package.Rd | 114 ++-- man/twCoefLogitnormMLE.Rd | 72 +-- man/twCoefLogitnormMLEFlat.Rd | 40 +- man/twSigmaLogitnorm.Rd | 60 +- tests/doRUnit.R | 53 +- vignettes/logitnorm.Rmd | 268 ++++---- 21 files changed, 1172 insertions(+), 1176 deletions(-) diff --git a/.Rbuildignore b/.Rbuildignore index 22416ae..fdd8e90 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,4 +1,4 @@ -^.*\.Rproj$ +^.*\.Rproj$ ^\.Rproj\.user$ ^README\.Rmd$ ^README\.md$ @@ -11,7 +11,6 @@ vignettes/cache ^tmp$ ^.git$ ^.gitignore$ -^.project$ # added by eclipse +^.project$ ^cranComments.md$ - ^\.travis\.yml$ diff --git a/.gitignore b/.gitignore index 09e0919..8846cf8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,36 +1,28 @@ # History files .Rhistory .Rapp.history - # Session Data files .RData - # Example code in package build process *-Ex.R - # Output files from R CMD build /*.tar.gz - # Output files from R CMD check /*.Rcheck/ - # RStudio files .Rproj.user/ - # produced vignettes vignettes/*.html vignettes/*.pdf - # OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3 .httr-oauth - # knitr and R markdown default cache directories /*_cache/ /cache/ - # Temporary files created by R markdown *.utf8.md *.knit.md /.project /.svn - +.Rproj.user +*.Rproj diff --git a/.travis.yml b/.travis.yml index d3e94d7..b4c1213 100644 --- a/.travis.yml +++ b/.travis.yml @@ -1,8 +1,8 @@ -# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r - -language: R -sudo: false -cache: packages - -git: - depth: 3 +# R for travis: see documentation at https://docs.travis-ci.com/user/languages/r + +language: R +sudo: false +cache: packages + +git: + depth: 3 diff --git a/DESCRIPTION b/DESCRIPTION index 5fa3317..8889354 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: logitnorm Title: Functions for the Logitnormal Distribution -Version: 0.8.34 +Version: 0.8.35 Date: 2017-03-14 Author: Thomas Wutzler Maintainer: Thomas Wutzler diff --git a/NAMESPACE b/NAMESPACE index 0ca60f4..718a381 100755 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,4 +1,4 @@ -exportPattern("^[^\\.]") #exports all variables that do not start with a period. - -importFrom("stats", "dnorm", "integrate", "optim", "optimize", - "plogis", "pnorm", "qlogis", "qnorm", "rnorm") +exportPattern("^[^\\.]") #exports all variables that do not start with a period. + +importFrom("stats", "dnorm", "integrate", "optim", "optimize", + "plogis", "pnorm", "qlogis", "qnorm", "rnorm") diff --git a/NEWS.md b/NEWS.md index 608cf0b..6d4bb04 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,6 @@ +# logitnorm 0.8.35 +Avoid writing file reports during testing and installation. + # logitnorm 0.8.33 ## New features diff --git a/R/logitnorm.R b/R/logitnorm.R index 793e868..093d373 100755 --- a/R/logitnorm.R +++ b/R/logitnorm.R @@ -1,454 +1,454 @@ -# random, percentile, density and quantile function of the logit-normal distribution -# and estimation of parameters from percentiles by Sum of Squares Newton optimization -# - -logit <- function( - ### Transforming (0,1) to normal scale (-Inf Inf) - p,... -){ - ##details<< - ## function \eqn{ logit(p)= log \left( \frac{p}{1-p} \right) = log(p) - log(1-p) } - ##seealso<< \code{\link{invlogit}} - ##seealso<< \code{\link{logitnorm}} - qlogis(p,...) -} - -invlogit <- function( - ### Transforming (-Inf,Inf) to original scale (0,1) - q,... -){ - ##details<< - ## function \eqn{f(z) = \frac{e^{z}}{e^{z} + 1} \! = \frac{1}{1 + e^{-z}} \!} - ##seealso<< \code{\link{logit}} - ##seealso<< \code{\link{logitnorm}} - plogis(q,...) -} - -# for normal-logit transformation do not use scale and location, better use it in generating rnorm etc - -rlogitnorm <- function( - ### Random number generation for logitnormal distribution - mu = 0, sigma = 1, ##<< distribution parameters - ... ##<< arguments to \code{\link{rnorm}} -){ - ##seealso<< \code{\link{logitnorm}} - plogis( rnorm(mean=mu,sd=sigma,...) ) -} - -plogitnorm <- function( - ### Distribution function for logitnormal distribution - q - ,mu = 0, sigma = 1 ##<< distribution parameters - ,... -){ - ##seealso<< \code{\link{logitnorm}} - ql <- qlogis(q) - pnorm(ql,mean=mu,sd=sigma,...) -} - - -dlogitnorm <- function( - ### Density function of logitnormal distribution - q ##<< quantiles - ,mu = 0, sigma = 1 ##<< distribution parameters - ,log = FALSE ##<< if TRUE, the log-density is returned - ,... ##<< further arguments passed to \code{\link{dnorm}}: \code{mean}, and \code{sd} for mu and sigma respectively. -){ - ##details<< \describe{\item{Logitnorm distribution}{ - ## \itemize{ - ## \item{density function: dlogitnorm } - ## \item{distribution function: \code{\link{plogitnorm}} } - ## \item{quantile function: \code{\link{qlogitnorm}} } - ## \item{random generation function: \code{\link{rlogitnorm}} } - ## } - ## }} - ##seealso<< \code{\link{logitnorm}} - - ql <- qlogis(q) - # multiply (or add in the log domain) by the Jacobian (derivative) of the back-Transformation (logit) - if (log) { - dnorm(ql,mean=mu,sd=sigma,log=TRUE,...) - log(q) - log1p(-q) - } else { - dnorm(ql,mean=mu,sd=sigma,...) /q/(1-q) - } -} - -qlogitnorm <- function( - ### Quantiles of logitnormal distribution. - p - ,mu = 0, sigma = 1 ##<< distribution parameters - ,... -){ - ##seealso<< \code{\link{logitnorm}} - qn <- qnorm(p,mean=mu,sd=sigma,...) - plogis(qn) -} - - -#------------------ estimating parameters from fitting to distribution -.ofLogitnorm <- function( - ### Objective function used by \code{\link{coefLogitnorm}}. - theta ##<< theta[1] is mu, theta[2] is the sigma - ,quant ##<< q are the quantiles for perc - ,perc=c(0.5,0.975) -){ - # there is an analytical expression for the gradient, but here just use numerical gradient - - # calculating perc should be faster than calculating quantiles of the normal distr. - if( theta[2] <= 0) return( Inf ) - tmp.predp = plogitnorm(quant, mu=theta[1], sigma=theta[2] ) - tmp.diff = tmp.predp - perc - - #however, q not so long valley and less calls to objective function - #but q has problems with high sigma - #tmp.predq = qlogitnorm(perc, mean=theta[1], sigma=theta[2] ) - #tmp.diff = tmp.predq - quant - sum(tmp.diff^2) -} - -twCoefLogitnormN <- function( - ### Estimating coefficients of logitnormal distribution from a vector of quantiles and perentiles (non-vectorized). - quant ##<< the quantile values - ,perc=c(0.5,0.975) ##<< the probabilites for which the quantiles were specified - ,method="BFGS" ##<< method of optimization (see \code{\link{optim}}) - ,theta0=c(mu=0,sigma=1) ##<< starting parameters - ,returnDetails=FALSE ##<< if TRUE, the full output of optim is returned instead of only entry par - , ... ##<< further parameters passed to optim, e.g. control=list(maxit=1000) -){ - ##seealso<< \code{\link{logitnorm}} - names(theta0) <- c("mu","sigma") - tmp <- optim( theta0, .ofLogitnorm, quant=quant,perc=perc, method=method, ...) - if( tmp$convergence != 0) - warning(paste("coefLogitnorm: optim did not converge. theta0=",paste(theta0,collapse=","))) - if( returnDetails ) - tmp - else - tmp$par - ### named numeric vector with estimated parameters of the logitnormal distrubtion. - ### names: \code{c("mu","sigma")} -} -#mtrace(coefLogitnorm) -attr(twCoefLogitnormN,"ex") <- function(){ - # experiment of re-estimation the parameters from generated observations - thetaTrue <- c(mu=0.8, sigma=0.7) - obsTrue <- rlogitnorm(thetaTrue["mu"],thetaTrue["sigma"], n=500) - obs <- obsTrue + rnorm(100, sd=0.05) # some observation uncertainty - plot(density(obsTrue),col="blue"); lines(density(obs)) - - # re-estimate parameters based on the quantiles of the observations - (theta <- twCoefLogitnorm( median(obs), quantile(obs,probs=0.9), perc = 0.9)) - - # add line of estimated distribution - x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid - dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) - lines( dx ~ x, col="orange") -} - -twCoefLogitnorm <- function( - ### Estimating coefficients of logitnormal distribution from median and upper quantile - median ##<< numeric vector: the median of the density function - ,quant ##<< numeric vector: the upper quantile value - ,perc=0.975 ##<< numeric vector: the probability for which the quantile was specified - ,method="BFGS" ##<< method of optimization (see \code{\link{optim}}) - ,theta0=c(mu=0,sigma=1) ##<< starting parameters - ,returnDetails=FALSE ##<< if TRUE, the full output of optim is attached as attributes resOptim - , ... -){ - ##seealso<< \code{\link{logitnorm}} - # twCoefLogitnorm - names(theta0) <- c("mu","sigma") - nc <- c(length(median),length(quant),length(perc)) - n <- max(nc) - res <- matrix( as.numeric(NA), n,2, dimnames=list(NULL,c("mu","sigma"))) - resOptim <- list() - qmat <- cbind(median,quant) - pmat <- cbind(0.5,perc) - for( i in 1:n ){ # for each row in (recycled) vector - i0 <- i-1 - tmp <- optim( theta0, .ofLogitnorm, perc=pmatI<-as.numeric(pmat[1+i0%%nrow(pmat),]), quant=(qmatI<-qmat[1+i0%%nrow(qmat),]), method=method, ...) - if( tmp$convergence == 0) - res[i,] <- tmp$par - else - warning(paste("coefLogitnorm: optim did not converge. theta0=",paste(theta0,collapse=","),"; median=",qmatI[1],"; quant=",qmatI[2],"; perc=",pmatI[2],sep="")) - resOptim[[i]] <- tmp - } - if( returnDetails ) attr(res,"resOptim") <- resOptim - res - ### numeric matrix with columns \code{c("mu","sigma")} - ### rows correspond to rows in median, quant, and perc -} -#mtrace(twCoefLogitnorm) -attr(twCoefLogitnorm,"ex") <- function(){ - # estimate the parameters, with median at 0.7 and upper quantile at 0.9 - (theta <- twCoefLogitnorm(0.7,0.9)) - - x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid - px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function - plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.5,0.975),col="gray") - - dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function - plot(dx~x); abline(v=c(0.7,0.9),col="gray") - - # vectorized - (theta <- twCoefLogitnorm(seq(0.4,0.8,by=0.1),0.9)) -} - -twCoefLogitnormCi <- function( - ### Calculates mu and sigma of the logitnormal distribution from lower and upper quantile, i.e. confidence interval. - lower ##<< value at the lower quantile, i.e. practical minimum - ,upper ##<< value at the upper quantile, i.e. practical maximum - ,perc=0.975 ##<< numeric vector: the probability for which the quantile was specified - ,sigmaFac=qnorm(perc) ##<< sigmaFac=2 is 95% sigmaFac=2.6 is 99% interval - ,isTransScale = FALSE ##<< if true lower and upper are already on logit scale -){ - ##seealso<< \code{\link{logitnorm}} - if( !isTRUE(isTransScale) ){ - lower <- logit(lower) - upper <- logit(upper) - } - halfWidth <- (upper-lower)/2 - sigma <- halfWidth/sigmaFac - cbind( mu=upper-halfWidth, sigma=sigma ) - ### named numeric vector: mu and sigma parameter of the logitnormal distribution. -} -attr(twCoefLogitnormCi,"ex") <- function(){ - mu=2 - sd=c(1,0.8) - p=0.99 - lower <- l <- qlogitnorm(1-p, mu, sd ) # p-confidence interval - upper <- u <- qlogitnorm(p, mu, sd ) # p-confidence interval - cf <- twCoefLogitnormCi(lower,upper) - all.equal( cf[,"mu"] , c(mu,mu) ) - all.equal( cf[,"sigma"] , sd ) -} - -.ofLogitnormMLE <- function( - ### Objective function used by \code{\link{coefLogitnormMLE}}. - mu ##<< numeric vector of proposed parameter - ## make sure that logit(mu)0.5 and logit(mu)>mle for mle<0.5 - ,mle ##<< the mode of the density distribution - ,logitMle=logit(mle)##<< may provide for performance reasons - ,quant ##<< q are the quantiles for perc - ,perc=c(0.975) -){ - # given mu and mle, we can calculate sigma - sigma2 = (logitMle-mu)/(2*mle-1) - ifelse( sigma2 <= 0, Inf, { - tmp.predp = plogitnorm(quant, mu=mu, sigma=sqrt(sigma2) ) - tmp.diff = tmp.predp - perc - tmp.diff^2 - }) -} - -twCoefLogitnormMLE <- function( - ### Estimating coefficients of logitnormal distribution from mode and upper quantile - mle ##<< numeric vector: the mode of the density function - ,quant ##<< numeric vector: the upper quantile value - ,perc=0.999 ##<< numeric vector: the probability for which the quantile was specified - #, ... ##<< Further parameters to \code{\link{optimize}} -){ - ##seealso<< \code{\link{logitnorm}} - # twCoefLogitnormMLE - nc <- c(length(mle),length(quant),length(perc)) - n <- max(nc) - res <- matrix( as.numeric(NA), n,2, dimnames=list(NULL,c("mu","sigma"))) - for( i in 1:n ){ - i0 <- i-1 - mleI<-mle[1+i0%%nc[1]] - if( mleI==0.5) mleI=0.5-.Machine$double.eps # in order to estimate sigma - # we now that mu is in (logit(mle),0) for mle < 0.5 and in (0,logit(mle)) for mle > 0.5 for unimodal distribution - # there might be a maximum in the middle and optimized misses the low part - # hence, first get near the global minimum by a evaluating the cost at a grid - # grid is spaced narrower at the edge - logitMle <- logit(mleI) - #intv <- if( logitMle < 0) c(logitMle+.Machine$double.eps,0) else c(0,logitMle-.Machine$double.eps) - upper <- abs(logitMle)-.Machine$double.eps - tmp <- sign(logitMle)*log(seq(1,exp(upper),length.out=40)) - oftmp<-.ofLogitnormMLE(tmp, mle=(mleI), logitMle=logitMle, quant=(quantI<-quant[1+i0%%nc[2]]), perc=(percI<-perc[1+i0%%nc[3]])) - #plot(tmp,oftmp) - # now optimize starting from the near minimum - imin <- which.min(oftmp) - intv <- tmp[ c(max(1,imin-1), min(length(tmp),imin+1)) ] - if( diff(intv)==0 ) mu <- intv[1] else - mu <- optimize( .ofLogitnormMLE, interval=intv, mle=(mleI), logitMle=logitMle, quant=(quantI<-quant[1+i0%%nc[2]]), perc=(percI<-perc[1+i0%%nc[3]]))$minimum - sigma <- sqrt((logitMle-mu)/(2*mleI-1)) - res[i,] <- c(mu,sigma) - } - res - ### numeric matrix with columns \code{c("mu","sigma")} - ### rows correspond to rows in mle, quant, and perc -} -#mtrace(coefLogitnorm) -attr(twCoefLogitnormMLE,"ex") <- function(){ - - # estimate the parameters, with mode 0.7 and upper quantile 0.9 - (theta <- twCoefLogitnormMLE(0.7,0.9)) - - x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid - px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function - plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.999),col="gray") - dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function - plot(dx~x); abline(v=c(0.7,0.9),col="gray") - - # vectorized - (theta <- twCoefLogitnormMLE(mle=seq(0.4,0.8,by=0.1),quant=0.9)) - - # flat - (theta <- twCoefLogitnormMLEFlat(0.7)) -} - -twCoefLogitnormMLEFlat <- function( - ### Estimating coefficients of a maximally flat unimodal logitnormal distribution from 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 - ## densitiy to a maximum range, i.e. is maximally flat, but still is unimodal. - twCoefLogitnormMLE(mle,0) -} - -.ofLogitnormE <- function( - ### Objective function used by \code{\link{coefLogitnormE}}. - theta ##<< theta[1] is the mu, theta[2] is the sigma - ,mean ##<< the expected value of the density distribution - ,quant ##<< q are the quantiles for perc - ,perc=c(0.975) - ,... ##<< further argument to \code{\link{momentsLogitnorm}}. -){ - tmp.predp = plogitnorm(quant, mu=theta[1], sigma=theta[2] ) - tmp.diff = tmp.predp - perc - .exp <- momentsLogitnorm(theta[1],theta[2],...)["mean"] - tmp.diff.e <- mean-.exp - sum(tmp.diff^2) + tmp.diff.e^2 -} - - -twCoefLogitnormE <- function( - ### Estimating coefficients of logitnormal distribution from expected value, i.e. mean, and upper quantile. - mean ##<< the expected value of the density function - ,quant ##<< the quantile values - ,perc=c(0.975) ##<< the probabilites for which the quantiles were specified - ,method="BFGS" ##<< method of optimization (see \code{\link{optim}}) - ,theta0=c(mu=0,sigma=1) ##<< starting parameters - ,returnDetails=FALSE ##<< if TRUE, the full output of optim is returned with attribut resOptim - , ... -){ - ##seealso<< \code{\link{logitnorm}} - # twCoefLogitnormE - names(theta0) <- c("mu","sigma") - mqp <- cbind(mean,quant,perc) - n <- nrow(mqp) - res <- matrix( as.numeric(NA), n,2, dimnames=list(NULL,c("mu","sigma"))) - resOptim <- list() - for( i in 1:n ){ - mqpi<- mqp[i,] - tmp <- optim( theta0, .ofLogitnormE, mean=mqpi[1], quant=mqpi[2], perc=mqpi[3], method=method, ...) - resOptim[[i]] <- tmp - if( tmp$convergence != 0) - warning(paste("coefLogitnorm: optim did not converge. theta0=",paste(theta0,collapse=","))) - else - res[i,] <- tmp$par - } - if( returnDetails ) - attr(res,"resOptim") <- resOptim - res - ### named numeric matrix with estimated parameters of the logitnormal distrubtion. - ### colnames: \code{c("mu","sigma")} -} -#mtrace(coefLogitnorm) -attr(twCoefLogitnormE,"ex") <- function(){ - # estimate the parameters - (thetaE <- twCoefLogitnormE(0.7,0.9)) - - x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid - px <- plogitnorm(x,mu=thetaE[1],sigma=thetaE[2]) #percentiles function - plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.5,0.975),col="gray") - dx <- dlogitnorm(x,mu=thetaE[1],sigma=thetaE[2]) #density function - plot(dx~x); abline(v=c(0.7,0.9),col="gray") - - z <- rlogitnorm(1e5, mu=thetaE[1],sigma=thetaE[2]) - mean(z) # about 0.7 - - # vectorized - (theta <- twCoefLogitnormE(mean=seq(0.4,0.8,by=0.1),quant=0.9)) -} - -momentsLogitnorm <- function( - ### First two moments of the logitnormal distribution by numerical integration - mu ##<< parameter mu - ,sigma ##<< parameter sigma - ,abs.tol=0 ##<< chaning default to \code{\link{integrate}} - ,... ##<< further parameters to the \code{\link{integrate}} function -){ - fExp <- function(x) plogis(x)*dnorm(x,mean=mu,sd=sigma) - .exp <- integrate(fExp,-Inf,Inf, abs.tol=abs.tol, ...)$value - fVar <- function(x) (plogis(x) - .exp)^2 * dnorm(x,mean=mu,sd=sigma) - .var <- integrate(fVar,-Inf,Inf, abs.tol=abs.tol, ...)$value - c( mean=.exp, var=.var ) - ### named numeric vector with components \itemize{ - ### \item{ \code{mean}: expected value, i.e. first moment} - ### \item{ \code{var}: variance, i.e. second moment } - ### } -} -attr(momentsLogitnorm,"ex") <- function(){ - (res <- momentsLogitnorm(4,1)) - (res <- momentsLogitnorm(5,0.1)) -} - -.ofModeLogitnorm <- function( - ### Objective function used by \code{\link{coefLogitnormMLE}}. - mle ##<< proposed mode - ,mu ##<< parameter mu - ,sd2 ##<< parameter sigma^2 -){ - if( mle <=0 | mle >=1) return(.Machine$double.xmax) - tmp.diff.mle <- mu-sd2*(1-2*mle) - logit(mle) - tmp.diff.mle^2 -} - -modeLogitnorm <- function( - ### Mode of the logitnormal distribution by numerical optimization - mu ##<< parameter mu - ,sigma ##<< parameter sigma - ,tol=invlogit(mu)/1000 ##<< precisions of the estimate -){ - ##seealso<< \code{\link{logitnorm}} - itv <- if(mu<0) c(0,invlogit(mu)) else c(invlogit(mu),1) - tmp <- optimize( .ofModeLogitnorm, interval=itv, mu=mu, sd2=sigma^2, tol=tol) - tmp$minimum -} - -twSigmaLogitnorm <- function( - ### Estimating coefficients of logitnormal distribution from mode and given mu - mle ##<< numeric vector: the mode of the density function - ,mu=0 ##<< for mu=0 the distribution will be the flattest case (maybe bimodal) -){ - ##details<< - ## For a mostly flat unimodal distribution use \code{\link{twCoefLogitnormMLE}(mle,0)} - sigma = sqrt( (logit(mle)-mu)/(2*mle-1) ) - ##seealso<< \code{\link{logitnorm}} - # twSigmaLogitnorm - cbind(mu=mu, sigma=sigma) - ### numeric matrix with columns \code{c("mu","sigma")} - ### rows correspond to rows in mle and mu -} -#mtrace(coefLogitnorm) -attr(twSigmaLogitnorm,"ex") <- function(){ - mle <- 0.8 - (theta <- twSigmaLogitnorm(mle)) - # - x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid - px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function - plot(px~x); abline(v=c(mle),col="gray") - dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function - plot(dx~x); abline(v=c(mle),col="gray") - # vectorized - (theta <- twSigmaLogitnorm(mle=seq(0.401,0.8,by=0.1))) -} - - - - - +# random, percentile, density and quantile function of the logit-normal distribution +# and estimation of parameters from percentiles by Sum of Squares Newton optimization +# + +logit <- function( + ### Transforming (0,1) to normal scale (-Inf Inf) + p,... +){ + ##details<< + ## function \eqn{ logit(p)= log \left( \frac{p}{1-p} \right) = log(p) - log(1-p) } + ##seealso<< \code{\link{invlogit}} + ##seealso<< \code{\link{logitnorm}} + qlogis(p,...) +} + +invlogit <- function( + ### Transforming (-Inf,Inf) to original scale (0,1) + q,... +){ + ##details<< + ## function \eqn{f(z) = \frac{e^{z}}{e^{z} + 1} \! = \frac{1}{1 + e^{-z}} \!} + ##seealso<< \code{\link{logit}} + ##seealso<< \code{\link{logitnorm}} + plogis(q,...) +} + +# for normal-logit transformation do not use scale and location, better use it in generating rnorm etc + +rlogitnorm <- function( + ### Random number generation for logitnormal distribution + mu = 0, sigma = 1, ##<< distribution parameters + ... ##<< arguments to \code{\link{rnorm}} +){ + ##seealso<< \code{\link{logitnorm}} + plogis( rnorm(mean=mu,sd=sigma,...) ) +} + +plogitnorm <- function( + ### Distribution function for logitnormal distribution + q + ,mu = 0, sigma = 1 ##<< distribution parameters + ,... +){ + ##seealso<< \code{\link{logitnorm}} + ql <- qlogis(q) + pnorm(ql,mean=mu,sd=sigma,...) +} + + +dlogitnorm <- function( + ### Density function of logitnormal distribution + q ##<< quantiles + ,mu = 0, sigma = 1 ##<< distribution parameters + ,log = FALSE ##<< if TRUE, the log-density is returned + ,... ##<< further arguments passed to \code{\link{dnorm}}: \code{mean}, and \code{sd} for mu and sigma respectively. +){ + ##details<< \describe{\item{Logitnorm distribution}{ + ## \itemize{ + ## \item{density function: dlogitnorm } + ## \item{distribution function: \code{\link{plogitnorm}} } + ## \item{quantile function: \code{\link{qlogitnorm}} } + ## \item{random generation function: \code{\link{rlogitnorm}} } + ## } + ## }} + ##seealso<< \code{\link{logitnorm}} + + ql <- qlogis(q) + # multiply (or add in the log domain) by the Jacobian (derivative) of the back-Transformation (logit) + if (log) { + dnorm(ql,mean=mu,sd=sigma,log=TRUE,...) - log(q) - log1p(-q) + } else { + dnorm(ql,mean=mu,sd=sigma,...) /q/(1-q) + } +} + +qlogitnorm <- function( + ### Quantiles of logitnormal distribution. + p + ,mu = 0, sigma = 1 ##<< distribution parameters + ,... +){ + ##seealso<< \code{\link{logitnorm}} + qn <- qnorm(p,mean=mu,sd=sigma,...) + plogis(qn) +} + + +#------------------ estimating parameters from fitting to distribution +.ofLogitnorm <- function( + ### Objective function used by \code{\link{coefLogitnorm}}. + theta ##<< theta[1] is mu, theta[2] is the sigma + ,quant ##<< q are the quantiles for perc + ,perc=c(0.5,0.975) +){ + # there is an analytical expression for the gradient, but here just use numerical gradient + + # calculating perc should be faster than calculating quantiles of the normal distr. + if( theta[2] <= 0) return( Inf ) + tmp.predp = plogitnorm(quant, mu=theta[1], sigma=theta[2] ) + tmp.diff = tmp.predp - perc + + #however, q not so long valley and less calls to objective function + #but q has problems with high sigma + #tmp.predq = qlogitnorm(perc, mean=theta[1], sigma=theta[2] ) + #tmp.diff = tmp.predq - quant + sum(tmp.diff^2) +} + +twCoefLogitnormN <- function( + ### Estimating coefficients of logitnormal distribution from a vector of quantiles and perentiles (non-vectorized). + quant ##<< the quantile values + ,perc=c(0.5,0.975) ##<< the probabilites for which the quantiles were specified + ,method="BFGS" ##<< method of optimization (see \code{\link{optim}}) + ,theta0=c(mu=0,sigma=1) ##<< starting parameters + ,returnDetails=FALSE ##<< if TRUE, the full output of optim is returned instead of only entry par + , ... ##<< further parameters passed to optim, e.g. control=list(maxit=1000) +){ + ##seealso<< \code{\link{logitnorm}} + names(theta0) <- c("mu","sigma") + tmp <- optim( theta0, .ofLogitnorm, quant=quant,perc=perc, method=method, ...) + if( tmp$convergence != 0) + warning(paste("coefLogitnorm: optim did not converge. theta0=",paste(theta0,collapse=","))) + if( returnDetails ) + tmp + else + tmp$par + ### named numeric vector with estimated parameters of the logitnormal distrubtion. + ### names: \code{c("mu","sigma")} +} +#mtrace(coefLogitnorm) +attr(twCoefLogitnormN,"ex") <- function(){ + # experiment of re-estimation the parameters from generated observations + thetaTrue <- c(mu=0.8, sigma=0.7) + obsTrue <- rlogitnorm(thetaTrue["mu"],thetaTrue["sigma"], n=500) + obs <- obsTrue + rnorm(100, sd=0.05) # some observation uncertainty + plot(density(obsTrue),col="blue"); lines(density(obs)) + + # re-estimate parameters based on the quantiles of the observations + (theta <- twCoefLogitnorm( median(obs), quantile(obs,probs=0.9), perc = 0.9)) + + # add line of estimated distribution + x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid + dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) + lines( dx ~ x, col="orange") +} + +twCoefLogitnorm <- function( + ### Estimating coefficients of logitnormal distribution from median and upper quantile + median ##<< numeric vector: the median of the density function + ,quant ##<< numeric vector: the upper quantile value + ,perc=0.975 ##<< numeric vector: the probability for which the quantile was specified + ,method="BFGS" ##<< method of optimization (see \code{\link{optim}}) + ,theta0=c(mu=0,sigma=1) ##<< starting parameters + ,returnDetails=FALSE ##<< if TRUE, the full output of optim is attached as attributes resOptim + , ... +){ + ##seealso<< \code{\link{logitnorm}} + # twCoefLogitnorm + names(theta0) <- c("mu","sigma") + nc <- c(length(median),length(quant),length(perc)) + n <- max(nc) + res <- matrix( as.numeric(NA), n,2, dimnames=list(NULL,c("mu","sigma"))) + resOptim <- list() + qmat <- cbind(median,quant) + pmat <- cbind(0.5,perc) + for( i in 1:n ){ # for each row in (recycled) vector + i0 <- i-1 + tmp <- optim( theta0, .ofLogitnorm, perc=pmatI<-as.numeric(pmat[1+i0%%nrow(pmat),]), quant=(qmatI<-qmat[1+i0%%nrow(qmat),]), method=method, ...) + if( tmp$convergence == 0) + res[i,] <- tmp$par + else + warning(paste("coefLogitnorm: optim did not converge. theta0=",paste(theta0,collapse=","),"; median=",qmatI[1],"; quant=",qmatI[2],"; perc=",pmatI[2],sep="")) + resOptim[[i]] <- tmp + } + if( returnDetails ) attr(res,"resOptim") <- resOptim + res + ### numeric matrix with columns \code{c("mu","sigma")} + ### rows correspond to rows in median, quant, and perc +} +#mtrace(twCoefLogitnorm) +attr(twCoefLogitnorm,"ex") <- function(){ + # estimate the parameters, with median at 0.7 and upper quantile at 0.9 + (theta <- twCoefLogitnorm(0.7,0.9)) + + x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid + px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function + plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.5,0.975),col="gray") + + dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function + plot(dx~x); abline(v=c(0.7,0.9),col="gray") + + # vectorized + (theta <- twCoefLogitnorm(seq(0.4,0.8,by=0.1),0.9)) +} + +twCoefLogitnormCi <- function( + ### Calculates mu and sigma of the logitnormal distribution from lower and upper quantile, i.e. confidence interval. + lower ##<< value at the lower quantile, i.e. practical minimum + ,upper ##<< value at the upper quantile, i.e. practical maximum + ,perc=0.975 ##<< numeric vector: the probability for which the quantile was specified + ,sigmaFac=qnorm(perc) ##<< sigmaFac=2 is 95% sigmaFac=2.6 is 99% interval + ,isTransScale = FALSE ##<< if true lower and upper are already on logit scale +){ + ##seealso<< \code{\link{logitnorm}} + if( !isTRUE(isTransScale) ){ + lower <- logit(lower) + upper <- logit(upper) + } + halfWidth <- (upper-lower)/2 + sigma <- halfWidth/sigmaFac + cbind( mu=upper-halfWidth, sigma=sigma ) + ### named numeric vector: mu and sigma parameter of the logitnormal distribution. +} +attr(twCoefLogitnormCi,"ex") <- function(){ + mu=2 + sd=c(1,0.8) + p=0.99 + lower <- l <- qlogitnorm(1-p, mu, sd ) # p-confidence interval + upper <- u <- qlogitnorm(p, mu, sd ) # p-confidence interval + cf <- twCoefLogitnormCi(lower,upper) + all.equal( cf[,"mu"] , c(mu,mu) ) + all.equal( cf[,"sigma"] , sd ) +} + +.ofLogitnormMLE <- function( + ### Objective function used by \code{\link{coefLogitnormMLE}}. + mu ##<< numeric vector of proposed parameter + ## make sure that logit(mu)0.5 and logit(mu)>mle for mle<0.5 + ,mle ##<< the mode of the density distribution + ,logitMle=logit(mle)##<< may provide for performance reasons + ,quant ##<< q are the quantiles for perc + ,perc=c(0.975) +){ + # given mu and mle, we can calculate sigma + sigma2 = (logitMle-mu)/(2*mle-1) + ifelse( sigma2 <= 0, Inf, { + tmp.predp = plogitnorm(quant, mu=mu, sigma=sqrt(sigma2) ) + tmp.diff = tmp.predp - perc + tmp.diff^2 + }) +} + +twCoefLogitnormMLE <- function( + ### Estimating coefficients of logitnormal distribution from mode and upper quantile + mle ##<< numeric vector: the mode of the density function + ,quant ##<< numeric vector: the upper quantile value + ,perc=0.999 ##<< numeric vector: the probability for which the quantile was specified + #, ... ##<< Further parameters to \code{\link{optimize}} +){ + ##seealso<< \code{\link{logitnorm}} + # twCoefLogitnormMLE + nc <- c(length(mle),length(quant),length(perc)) + n <- max(nc) + res <- matrix( as.numeric(NA), n,2, dimnames=list(NULL,c("mu","sigma"))) + for( i in 1:n ){ + i0 <- i-1 + mleI<-mle[1+i0%%nc[1]] + if( mleI==0.5) mleI=0.5-.Machine$double.eps # in order to estimate sigma + # we now that mu is in (logit(mle),0) for mle < 0.5 and in (0,logit(mle)) for mle > 0.5 for unimodal distribution + # there might be a maximum in the middle and optimized misses the low part + # hence, first get near the global minimum by a evaluating the cost at a grid + # grid is spaced narrower at the edge + logitMle <- logit(mleI) + #intv <- if( logitMle < 0) c(logitMle+.Machine$double.eps,0) else c(0,logitMle-.Machine$double.eps) + upper <- abs(logitMle)-.Machine$double.eps + tmp <- sign(logitMle)*log(seq(1,exp(upper),length.out=40)) + oftmp<-.ofLogitnormMLE(tmp, mle=(mleI), logitMle=logitMle, quant=(quantI<-quant[1+i0%%nc[2]]), perc=(percI<-perc[1+i0%%nc[3]])) + #plot(tmp,oftmp) + # now optimize starting from the near minimum + imin <- which.min(oftmp) + intv <- tmp[ c(max(1,imin-1), min(length(tmp),imin+1)) ] + if( diff(intv)==0 ) mu <- intv[1] else + mu <- optimize( .ofLogitnormMLE, interval=intv, mle=(mleI), logitMle=logitMle, quant=(quantI<-quant[1+i0%%nc[2]]), perc=(percI<-perc[1+i0%%nc[3]]))$minimum + sigma <- sqrt((logitMle-mu)/(2*mleI-1)) + res[i,] <- c(mu,sigma) + } + res + ### numeric matrix with columns \code{c("mu","sigma")} + ### rows correspond to rows in mle, quant, and perc +} +#mtrace(coefLogitnorm) +attr(twCoefLogitnormMLE,"ex") <- function(){ + + # estimate the parameters, with mode 0.7 and upper quantile 0.9 + (theta <- twCoefLogitnormMLE(0.7,0.9)) + + x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid + px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function + plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.999),col="gray") + dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function + plot(dx~x); abline(v=c(0.7,0.9),col="gray") + + # vectorized + (theta <- twCoefLogitnormMLE(mle=seq(0.4,0.8,by=0.1),quant=0.9)) + + # flat + (theta <- twCoefLogitnormMLEFlat(0.7)) +} + +twCoefLogitnormMLEFlat <- function( + ### Estimating coefficients of a maximally flat unimodal logitnormal distribution from 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 + ## densitiy to a maximum range, i.e. is maximally flat, but still is unimodal. + twCoefLogitnormMLE(mle,0) +} + +.ofLogitnormE <- function( + ### Objective function used by \code{\link{coefLogitnormE}}. + theta ##<< theta[1] is the mu, theta[2] is the sigma + ,mean ##<< the expected value of the density distribution + ,quant ##<< q are the quantiles for perc + ,perc=c(0.975) + ,... ##<< further argument to \code{\link{momentsLogitnorm}}. +){ + tmp.predp = plogitnorm(quant, mu=theta[1], sigma=theta[2] ) + tmp.diff = tmp.predp - perc + .exp <- momentsLogitnorm(theta[1],theta[2],...)["mean"] + tmp.diff.e <- mean-.exp + sum(tmp.diff^2) + tmp.diff.e^2 +} + + +twCoefLogitnormE <- function( + ### Estimating coefficients of logitnormal distribution from expected value, i.e. mean, and upper quantile. + mean ##<< the expected value of the density function + ,quant ##<< the quantile values + ,perc=c(0.975) ##<< the probabilites for which the quantiles were specified + ,method="BFGS" ##<< method of optimization (see \code{\link{optim}}) + ,theta0=c(mu=0,sigma=1) ##<< starting parameters + ,returnDetails=FALSE ##<< if TRUE, the full output of optim is returned with attribut resOptim + , ... +){ + ##seealso<< \code{\link{logitnorm}} + # twCoefLogitnormE + names(theta0) <- c("mu","sigma") + mqp <- cbind(mean,quant,perc) + n <- nrow(mqp) + res <- matrix( as.numeric(NA), n,2, dimnames=list(NULL,c("mu","sigma"))) + resOptim <- list() + for( i in 1:n ){ + mqpi<- mqp[i,] + tmp <- optim( theta0, .ofLogitnormE, mean=mqpi[1], quant=mqpi[2], perc=mqpi[3], method=method, ...) + resOptim[[i]] <- tmp + if( tmp$convergence != 0) + warning(paste("coefLogitnorm: optim did not converge. theta0=",paste(theta0,collapse=","))) + else + res[i,] <- tmp$par + } + if( returnDetails ) + attr(res,"resOptim") <- resOptim + res + ### named numeric matrix with estimated parameters of the logitnormal distrubtion. + ### colnames: \code{c("mu","sigma")} +} +#mtrace(coefLogitnorm) +attr(twCoefLogitnormE,"ex") <- function(){ + # estimate the parameters + (thetaE <- twCoefLogitnormE(0.7,0.9)) + + x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid + px <- plogitnorm(x,mu=thetaE[1],sigma=thetaE[2]) #percentiles function + plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.5,0.975),col="gray") + dx <- dlogitnorm(x,mu=thetaE[1],sigma=thetaE[2]) #density function + plot(dx~x); abline(v=c(0.7,0.9),col="gray") + + z <- rlogitnorm(1e5, mu=thetaE[1],sigma=thetaE[2]) + mean(z) # about 0.7 + + # vectorized + (theta <- twCoefLogitnormE(mean=seq(0.4,0.8,by=0.1),quant=0.9)) +} + +momentsLogitnorm <- function( + ### First two moments of the logitnormal distribution by numerical integration + mu ##<< parameter mu + ,sigma ##<< parameter sigma + ,abs.tol=0 ##<< chaning default to \code{\link{integrate}} + ,... ##<< further parameters to the \code{\link{integrate}} function +){ + fExp <- function(x) plogis(x)*dnorm(x,mean=mu,sd=sigma) + .exp <- integrate(fExp,-Inf,Inf, abs.tol=abs.tol, ...)$value + fVar <- function(x) (plogis(x) - .exp)^2 * dnorm(x,mean=mu,sd=sigma) + .var <- integrate(fVar,-Inf,Inf, abs.tol=abs.tol, ...)$value + c( mean=.exp, var=.var ) + ### named numeric vector with components \itemize{ + ### \item{ \code{mean}: expected value, i.e. first moment} + ### \item{ \code{var}: variance, i.e. second moment } + ### } +} +attr(momentsLogitnorm,"ex") <- function(){ + (res <- momentsLogitnorm(4,1)) + (res <- momentsLogitnorm(5,0.1)) +} + +.ofModeLogitnorm <- function( + ### Objective function used by \code{\link{coefLogitnormMLE}}. + mle ##<< proposed mode + ,mu ##<< parameter mu + ,sd2 ##<< parameter sigma^2 +){ + if( mle <=0 | mle >=1) return(.Machine$double.xmax) + tmp.diff.mle <- mu-sd2*(1-2*mle) - logit(mle) + tmp.diff.mle^2 +} + +modeLogitnorm <- function( + ### Mode of the logitnormal distribution by numerical optimization + mu ##<< parameter mu + ,sigma ##<< parameter sigma + ,tol=invlogit(mu)/1000 ##<< precisions of the estimate +){ + ##seealso<< \code{\link{logitnorm}} + itv <- if(mu<0) c(0,invlogit(mu)) else c(invlogit(mu),1) + tmp <- optimize( .ofModeLogitnorm, interval=itv, mu=mu, sd2=sigma^2, tol=tol) + tmp$minimum +} + +twSigmaLogitnorm <- function( + ### Estimating coefficients of logitnormal distribution from mode and given mu + mle ##<< numeric vector: the mode of the density function + ,mu=0 ##<< for mu=0 the distribution will be the flattest case (maybe bimodal) +){ + ##details<< + ## For a mostly flat unimodal distribution use \code{\link{twCoefLogitnormMLE}(mle,0)} + sigma = sqrt( (logit(mle)-mu)/(2*mle-1) ) + ##seealso<< \code{\link{logitnorm}} + # twSigmaLogitnorm + cbind(mu=mu, sigma=sigma) + ### numeric matrix with columns \code{c("mu","sigma")} + ### rows correspond to rows in mle and mu +} +#mtrace(coefLogitnorm) +attr(twSigmaLogitnorm,"ex") <- function(){ + mle <- 0.8 + (theta <- twSigmaLogitnorm(mle)) + # + x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid + px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function + plot(px~x); abline(v=c(mle),col="gray") + dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function + plot(dx~x); abline(v=c(mle),col="gray") + # vectorized + (theta <- twSigmaLogitnorm(mle=seq(0.401,0.8,by=0.1))) +} + + + + + diff --git a/README.Rmd b/README.Rmd index d350a1d..51b9c6d 100644 --- a/README.Rmd +++ b/README.Rmd @@ -1,56 +1,56 @@ ---- -output: github_document ---- - - - - -```{r, echo = FALSE} -knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>", - fig.path = "tools/README-" -) -``` - -[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/logitnorm)](http://cran.r-project.org/package=logitnorm) -[![Travis-CI Build Status](https://travis-ci.org/bgctw/logitnorm.svg?branch=master)](https://travis-ci.org/bgctw/logitnorm) - - -## Overview - -`logitnorm` package provides support for the univariate -[logit-normal -distribution](https://en.wikipedia.org/wiki/Logit-normal_distribution). In -addition to the usual random, density, percential, and quantile function, it -helps with estimating distribution parameters from observations statistics. - -## Installation - -```{r, eval = FALSE} -# From CRAN -install.packages("logitnorm") - -# Or the the development version from GitHub: -# install.packages("devtools") -devtools::install_github("bgctw/logitnorm") -``` - -## Usage - -See the package vignette for an introduction. - -A simple example estimates distribution parameters from observation -statistics of mode 0.7 and upper quantile 0.9. Next, the density is -computed and plotted across a range of quantiles. - -```{r example} -(theta <- twCoefLogitnormMLE(0.7,0.9)) -x <- seq(0,1, length.out=81) -d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) -plot(d~x,type="l") -abline(v=c(0.7,0.9), col="grey") -``` +--- +output: github_document +--- + + + + +```{r, echo = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.path = "tools/README-" +) +``` + +[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/logitnorm)](http://cran.r-project.org/package=logitnorm) +[![Travis-CI Build Status](https://travis-ci.org/bgctw/logitnorm.svg?branch=master)](https://travis-ci.org/bgctw/logitnorm) + + +## Overview + +`logitnorm` package provides support for the univariate +[logit-normal +distribution](https://en.wikipedia.org/wiki/Logit-normal_distribution). In +addition to the usual random, density, percential, and quantile function, it +helps with estimating distribution parameters from observations statistics. + +## Installation + +```{r, eval = FALSE} +# From CRAN +install.packages("logitnorm") + +# Or the the development version from GitHub: +# install.packages("devtools") +devtools::install_github("bgctw/logitnorm") +``` + +## Usage + +See the package vignette for an introduction. + +A simple example estimates distribution parameters from observation +statistics of mode 0.7 and upper quantile 0.9. Next, the density is +computed and plotted across a range of quantiles. + +```{r example} +(theta <- twCoefLogitnormMLE(0.7,0.9)) +x <- seq(0,1, length.out=81) +d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) +plot(d~x,type="l") +abline(v=c(0.7,0.9), col="grey") +``` diff --git a/README.md b/README.md index c98845f..e6dbb33 100644 --- a/README.md +++ b/README.md @@ -1,56 +1,56 @@ ---- -output: github_document ---- - - - - - - -[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/logitnorm)](http://cran.r-project.org/package=logitnorm) -[![Travis-CI Build Status](https://travis-ci.org/bgctw/logitnorm.svg?branch=master)](https://travis-ci.org/bgctw/logitnorm) - - -## Overview - -`logitnorm` package provides support for the univariate -[logit-normal -distribution](https://en.wikipedia.org/wiki/Logit-normal_distribution). In -addition to the usual random, density, percential, and quantile function, it -helps with estimating distribution parameters from observations statistics. - -## Installation - - -```r -# From CRAN -install.packages("logitnorm") - -# Or the the development version from GitHub: -# install.packages("devtools") -devtools::install_github("bgctw/logitnorm") -``` - -## Usage - -See the package vignette for an introduction. - -A simple example estimates distribution parameters from observation -statistics of mode 0.7 and upper quantile 0.9. Next, the density is -computed and plotted across a range of quantiles. - - -```r -(theta <- twCoefLogitnormMLE(0.7,0.9)) -#> mu sigma -#> [1,] 0.7608886 0.464783 -x <- seq(0,1, length.out=81) -d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) -plot(d~x,type="l") -abline(v=c(0.7,0.9), col="grey") -``` - -![plot of chunk example](tools/README-example-1.png) +--- +output: github_document +--- + + + + + + +[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/logitnorm)](http://cran.r-project.org/package=logitnorm) +[![Travis-CI Build Status](https://travis-ci.org/bgctw/logitnorm.svg?branch=master)](https://travis-ci.org/bgctw/logitnorm) + + +## Overview + +`logitnorm` package provides support for the univariate +[logit-normal +distribution](https://en.wikipedia.org/wiki/Logit-normal_distribution). In +addition to the usual random, density, percential, and quantile function, it +helps with estimating distribution parameters from observations statistics. + +## Installation + + +```r +# From CRAN +install.packages("logitnorm") + +# Or the the development version from GitHub: +# install.packages("devtools") +devtools::install_github("bgctw/logitnorm") +``` + +## Usage + +See the package vignette for an introduction. + +A simple example estimates distribution parameters from observation +statistics of mode 0.7 and upper quantile 0.9. Next, the density is +computed and plotted across a range of quantiles. + + +```r +(theta <- twCoefLogitnormMLE(0.7,0.9)) +#> mu sigma +#> [1,] 0.7608886 0.464783 +x <- seq(0,1, length.out=81) +d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) +plot(d~x,type="l") +abline(v=c(0.7,0.9), col="grey") +``` + +![plot of chunk example](tools/README-example-1.png) diff --git a/cranComments.md b/cranComments.md index 8f8e151..7762336 100644 --- a/cranComments.md +++ b/cranComments.md @@ -1,14 +1,13 @@ ## Test environments -* local Windows install, R 3.2.4 64bit -* all checks from r-forge (https://r-forge.r-project.org/R/?group_id=886) +* local R3.4.3 on Mint 17 64bit +* Ubuntu 14.04 (on travis-ci), R 3.4.1 +* win-builder (devel and release) ## R CMD check results -There were no ERRORs nor WARNINGs +There were no ERRORs nor WARNINGs nor NOTEs -On r-forge - Windows there was one Note: - - Files 'README.md' or 'NEWS.md' cannot be checked without 'pandoc' being installed. - - these files are listed in .Rbuildignore +## Notes +Avoid writing file reports during testing and installation, in response +to email from Kurt Hornik Dez, 2nd 2017 \ No newline at end of file diff --git a/inst/genData/_genPackage.R b/inst/genData/_genPackage.R index 4f2ff6e..613bce5 100755 --- a/inst/genData/_genPackage.R +++ b/inst/genData/_genPackage.R @@ -1,19 +1,19 @@ -#------------ generating package tw.DEMC -# inlcuding generation of Rd files from inline-docs -# install.packages("inlinedocs") - -.tmp.f <- function(){ - library(twMisc) - source("R/logitnorm.R") - twUtestF() -} - -#R CMD build logitnorm -#R CMD check --as-cran logitnorm_XXX | tee tmp.txt | grep -v OK$ -#R CMD check --no-vignettes --no-latex --no-codoc logitnorm -#R CMD INSTALL --html logitnorm - -.tmp.f <- function(){ - library(sos) - fres1 <- findFn("logitnormal") -} +#------------ generating package tw.DEMC +# inlcuding generation of Rd files from inline-docs +# install.packages("inlinedocs") + +.tmp.f <- function(){ + library(twMisc) + source("R/logitnorm.R") + twUtestF() +} + +#R CMD build logitnorm +#R CMD check --as-cran logitnorm_XXX | tee tmp.txt | grep -v OK$ +#R CMD check --no-vignettes --no-latex --no-codoc logitnorm +#R CMD INSTALL --html logitnorm + +.tmp.f <- function(){ + library(sos) + fres1 <- findFn("logitnormal") +} diff --git a/inst/genData/logitnorm-package.Rd b/inst/genData/logitnorm-package.Rd index d6d5cae..a941b73 100644 --- a/inst/genData/logitnorm-package.Rd +++ b/inst/genData/logitnorm-package.Rd @@ -1,55 +1,55 @@ -\name{logitnorm-package} -\alias{logitnorm-package} \alias{logitnorm} -\title{Utilities for the logitnormal distribution in R} -\description{ Utilities for the logitnormal distribution in R -\itemize{ -\item Density, distribution, quantile and random generation function. -\item Estimation of the mode and the first two moments. -\item Estimation of distribution parameters from observations. -} -}%description - - -\details{ -The package provides the main distribution functions: -\itemize{ -\item density \code{\link{dlogitnorm}}, -\item distribution \code{\link{plogitnorm}}, -\item quantile \code{\link{qlogitnorm}}, and -\item random generation function \code{\link{rlogitnorm}}. -}%itemize - -Transformation functions - \itemize{ - \item{ (0,1) -> (-Inf,Inf): \code{\link{logit}} } - \item{ (-Inf,Inf) -> (0,1): \code{\link{invlogit}} } - } - -Moments and mode - \itemize{ - \item{ Expected value and variance: \code{\link{momentsLogitnorm}} } - \item{ Mode: \code{\link{modeLogitnorm}} } - } - -Estimating parameters - \itemize{ - \item{from mode and upper quantile: \code{\link{twCoefLogitnormMLE}} } - \item{from mode and constraint to be unimodal and maximally flat: \code{\link{twCoefLogitnormMLEFlat}} } - \item{from median and upper quantile: \code{\link{twCoefLogitnorm}} } - \item{from expected value, i.e. mean and upper quantile: \code{\link{twCoefLogitnormE}} } - \item{from a confidence interval which is symmetric at normal scale: \code{\link{twCoefLogitnormCi}} } - \item{from prescribed quantiles: \code{\link{twCoefLogitnormN}} } - } - -Have a look at the package vignettes. - }%details - -\references{ -Frederic, P. & Lad, F. (2008) Two Moments of the Logitnormal Distribution. -Communications in Statistics-Simulation and Computation, 37, 1263-1269 -}%references -\author{Thomas Wutzler} - -\keyword{ package } - - +\name{logitnorm-package} +\alias{logitnorm-package} \alias{logitnorm} +\title{Utilities for the logitnormal distribution in R} +\description{ Utilities for the logitnormal distribution in R +\itemize{ +\item Density, distribution, quantile and random generation function. +\item Estimation of the mode and the first two moments. +\item Estimation of distribution parameters from observations. +} +}%description + + +\details{ +The package provides the main distribution functions: +\itemize{ +\item density \code{\link{dlogitnorm}}, +\item distribution \code{\link{plogitnorm}}, +\item quantile \code{\link{qlogitnorm}}, and +\item random generation function \code{\link{rlogitnorm}}. +}%itemize + +Transformation functions + \itemize{ + \item{ (0,1) -> (-Inf,Inf): \code{\link{logit}} } + \item{ (-Inf,Inf) -> (0,1): \code{\link{invlogit}} } + } + +Moments and mode + \itemize{ + \item{ Expected value and variance: \code{\link{momentsLogitnorm}} } + \item{ Mode: \code{\link{modeLogitnorm}} } + } + +Estimating parameters + \itemize{ + \item{from mode and upper quantile: \code{\link{twCoefLogitnormMLE}} } + \item{from mode and constraint to be unimodal and maximally flat: \code{\link{twCoefLogitnormMLEFlat}} } + \item{from median and upper quantile: \code{\link{twCoefLogitnorm}} } + \item{from expected value, i.e. mean and upper quantile: \code{\link{twCoefLogitnormE}} } + \item{from a confidence interval which is symmetric at normal scale: \code{\link{twCoefLogitnormCi}} } + \item{from prescribed quantiles: \code{\link{twCoefLogitnormN}} } + } + +Have a look at the package vignettes. + }%details + +\references{ +Frederic, P. & Lad, F. (2008) Two Moments of the Logitnormal Distribution. +Communications in Statistics-Simulation and Computation, 37, 1263-1269 +}%references +\author{Thomas Wutzler} + +\keyword{ package } + + diff --git a/inst/unitTests/runitLogitnorm.R b/inst/unitTests/runitLogitnorm.R index 2ce9401..f210e66 100755 --- a/inst/unitTests/runitLogitnorm.R +++ b/inst/unitTests/runitLogitnorm.R @@ -1,149 +1,149 @@ -.setUp <-function () { - #library(MASS) - .setUpDf <- within( list(),{ - x <- seq(0,1,length.out=41)[-c(1,41)]; #x[1] = x[1] + .Machine$double.eps; x[length(x)] <- x[length(x)]- .Machine$double.eps - lx <- logit(x) - }) - attach(.setUpDf) -} - -.tearDown <- function () { - #detach(.setUpDf) - detach() -} - - - -test.inverseSame <- function(){ - xnorm <- logit(x) - xinv <- invlogit(xnorm) - checkEquals(x, xinv) -} - -test.plogitnorm <- function(){ - px <- plogitnorm(x) #percentiles - checkEquals( pnorm(lx), px) - px2 <- plogitnorm(x,mu=2,sigma=1) - checkEquals( pnorm(lx,mean=2,sd=1), px2) - #plot( px ~ x) - #plot( px ~ logit(x)) - #lines( pnorm( logit(x)) ~ logit(x) ) -} - -test.twCoefLogitnorm <- function(){ - theta <- twCoefLogitnorm(0.7,0.9,perc=0.999) - px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function - dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function - #plot(px~x); abline(v=c(0.7,0.9)); abline(h=c(0.5,0.975)) - #plot(dx~x); abline(v=c(0.7,0.9)) - # upper percentile at 0.9 - checkEquals(which.min(abs(px-0.999)), which(x==0.9) ) - checkEquals(which.min(abs(px-0.5)), which.min(abs(x-0.7)) ) - # mode at 0.7 - #checkEquals(which(abs(x-0.7)<.Machine$double.eps), which.max(dx) ) -} - -test.twCoefLogitnormN <- function(){ - quant=c(0.7,0.8,0.9) - perc=c(0.5,0.75,0.975) - (theta <- twCoefLogitnormN( quant=quant, perc=perc )) - #px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function - #dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function - #plot(px~x); abline(v=quant,col="gray"); abline(h=perc,col="gray") -} - - -test.twCoefLogitnormMLE <- function(){ - theta <- twCoefLogitnormMLE(0.7,0.9,perc=0.975) - px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function - dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function - #plot(px~x); abline(v=c(0.7,0.9)); abline(h=c(0.5,0.975)) - #plot(dx~x); abline(v=c(0.7,0.9)) - # upper percentile at 0.9 - checkEquals(which.min(abs(px-0.975)), which(x==0.9) ) - # mode at 0.7 - checkEquals(which.min(abs(x-0.7)), which.max(dx) ) -} - - -test.twCoefLogitnormE <- function(){ - theta <- twCoefLogitnormE(0.7,0.9) - px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function - dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function - #plot(px~x); abline(v=c(0.7,0.9)); abline(h=c(0.5,0.975)) - #plot(dx~x); abline(v=c(0.7,0.9)) - # upper percentile at 0.9 - checkEquals(which.min(abs(px-0.975)), which(x==0.9) ) - # mean at 0.7 - checkEqualsNumeric( momentsLogitnorm(mu=theta[1],sigma=theta[2])["mean"], 0.7, tolerance=1e-3) - z <- rlogitnorm(1e5, mu=theta[1],sigma=theta[2]) - checkEqualsNumeric(0.7, mean(z), tolerance=5e-3 ) -} - -.tmp.f <- function(){ - px <- plogitnorm(x) #percentiles - plot( px ~ x ) - plot( qlogitnorm(px) ~ x ) #one to one line - plot( dlogitnorm(x,mu=0.9) ~ x, type="l" ) - abline( v=qlogitnorm(c(0.025,0.5,0.975), mu=0.9)) -} - -.tmp.f <- function(){ - library(MASS) - ?fitdistr #not implemented - quant = c(0.6,0.9) - perc = c(0.5,0.975) - theta0=c(mu=0,sigma=1) - method="BFGS" - #mtrace(ofLogitnorm) - #popt <- as.list(tmp$par) - popt <- as.list(coefLogitnorm(quant)) - popt2 <- as.list(coefLogitnorm(quant, perc=c(0.5,0.9995))) - ofLogitnorm(popt,quant,perc) - - plot( dlogitnorm(x,mu=popt$mu,sigma=popt$sigma) ~ x, type="l" ) - abline( v=qlogitnorm(c(0.025,0.5,0.975),mu=popt$mu,sigma=popt$sigma)) - - lines( dlogitnorm(x,mu=popt2$mu,sigma=popt2$sigma) ~ x, type="l", col="maroon" ) - abline( v=qlogitnorm(c(0.5,0.9995),mu=popt2$mu,sigma=popt2$sigma), col="maroon") - - popt <- as.list(coefLogitnorm(c(0.9, 0.9995), perc=c(0.5,0.9995))) - plot( dlogitnorm(x,mu=popt$mu,sigma=popt$sigma) ~ x, type="l" ) - abline( v=qlogitnorm(c(0.025,0.5,0.975),mu=popt$mu,sigma=popt$sigma)) - -} - -.tmp.f <- function(){ - #visualize the objective functions surface - quant = c(0.6,0.9) - perc = c(0.5,0.975) - perc = c(0.5,0.9995) - perc=c(0.5,upperBoundProb) - quant = parms.var[varDist=="logitnorm",c("qMedian","qUpper")] - quant = parms.var["epsF",c("qMedian","qUpper")] - quant = parms.var["epsG",c("qMedian","qUpper")] - quant = parms.var["epsP",c("qMedian","qUpper")] - - - tmp.n <- 80 - tmp.mu <- seq(0,1,length.out=tmp.n) - tmp.sigma <- seq(0.01,2.5,length.out=tmp.n) - tmp <- as.matrix(expand.grid( mu=tmp.mu, sigma=tmp.sigma)) - tmp.of <- apply(tmp,1,ofLogitnorm, quant=quant,perc=perc ) - tmp.ofm <- matrix(tmp.of, nrow=tmp.n ) - image(tmp.mu, tmp.sigma, tmp.ofm) - image(tmp.mu, tmp.sigma, exp(-0.5*tmp.ofm)) #very flat - image(tmp.mu, tmp.sigma, exp(-0.5*1/(1/800)*tmp.ofm)) #distort by decreasing Temp *1/T - tmp.o <- coefLogitnorm(quant=quant, perc=perc, returnDetails=TRUE) - tmp.o - popt <- as.list(tmp.o$par) - points( popt$mu, popt$sigma) - - windows() - unlist(popt) - plot( dlogitnorm(x,mu=popt$mu,sigma=popt$sigma) ~ x, type="l" ) - abline( v=qlogitnorm(perc,mu=popt$mu,sigma=popt$sigma)) - lines( dlogitnorm(x,mu=popt$mu,sigma=1) ~ x, type="l", col="maroon" ) - abline( v=qlogitnorm(c(0.025,0.5,0.975),mu=popt$mu,sigma=1), col="maroon") -} - +.setUp <-function () { + #library(MASS) + .setUpDf <- within( list(),{ + x <- seq(0,1,length.out=41)[-c(1,41)]; #x[1] = x[1] + .Machine$double.eps; x[length(x)] <- x[length(x)]- .Machine$double.eps + lx <- logit(x) + }) + attach(.setUpDf) +} + +.tearDown <- function () { + #detach(.setUpDf) + detach() +} + + + +test.inverseSame <- function(){ + xnorm <- logit(x) + xinv <- invlogit(xnorm) + checkEquals(x, xinv) +} + +test.plogitnorm <- function(){ + px <- plogitnorm(x) #percentiles + checkEquals( pnorm(lx), px) + px2 <- plogitnorm(x,mu=2,sigma=1) + checkEquals( pnorm(lx,mean=2,sd=1), px2) + #plot( px ~ x) + #plot( px ~ logit(x)) + #lines( pnorm( logit(x)) ~ logit(x) ) +} + +test.twCoefLogitnorm <- function(){ + theta <- twCoefLogitnorm(0.7,0.9,perc=0.999) + px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function + dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function + #plot(px~x); abline(v=c(0.7,0.9)); abline(h=c(0.5,0.975)) + #plot(dx~x); abline(v=c(0.7,0.9)) + # upper percentile at 0.9 + checkEquals(which.min(abs(px-0.999)), which(x==0.9) ) + checkEquals(which.min(abs(px-0.5)), which.min(abs(x-0.7)) ) + # mode at 0.7 + #checkEquals(which(abs(x-0.7)<.Machine$double.eps), which.max(dx) ) +} + +test.twCoefLogitnormN <- function(){ + quant=c(0.7,0.8,0.9) + perc=c(0.5,0.75,0.975) + (theta <- twCoefLogitnormN( quant=quant, perc=perc )) + #px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function + #dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function + #plot(px~x); abline(v=quant,col="gray"); abline(h=perc,col="gray") +} + + +test.twCoefLogitnormMLE <- function(){ + theta <- twCoefLogitnormMLE(0.7,0.9,perc=0.975) + px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function + dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function + #plot(px~x); abline(v=c(0.7,0.9)); abline(h=c(0.5,0.975)) + #plot(dx~x); abline(v=c(0.7,0.9)) + # upper percentile at 0.9 + checkEquals(which.min(abs(px-0.975)), which(x==0.9) ) + # mode at 0.7 + checkEquals(which.min(abs(x-0.7)), which.max(dx) ) +} + + +test.twCoefLogitnormE <- function(){ + theta <- twCoefLogitnormE(0.7,0.9) + px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function + dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function + #plot(px~x); abline(v=c(0.7,0.9)); abline(h=c(0.5,0.975)) + #plot(dx~x); abline(v=c(0.7,0.9)) + # upper percentile at 0.9 + checkEquals(which.min(abs(px-0.975)), which(x==0.9) ) + # mean at 0.7 + checkEqualsNumeric( momentsLogitnorm(mu=theta[1],sigma=theta[2])["mean"], 0.7, tolerance=1e-3) + z <- rlogitnorm(1e5, mu=theta[1],sigma=theta[2]) + checkEqualsNumeric(0.7, mean(z), tolerance=5e-3 ) +} + +.tmp.f <- function(){ + px <- plogitnorm(x) #percentiles + plot( px ~ x ) + plot( qlogitnorm(px) ~ x ) #one to one line + plot( dlogitnorm(x,mu=0.9) ~ x, type="l" ) + abline( v=qlogitnorm(c(0.025,0.5,0.975), mu=0.9)) +} + +.tmp.f <- function(){ + library(MASS) + ?fitdistr #not implemented + quant = c(0.6,0.9) + perc = c(0.5,0.975) + theta0=c(mu=0,sigma=1) + method="BFGS" + #mtrace(ofLogitnorm) + #popt <- as.list(tmp$par) + popt <- as.list(coefLogitnorm(quant)) + popt2 <- as.list(coefLogitnorm(quant, perc=c(0.5,0.9995))) + ofLogitnorm(popt,quant,perc) + + plot( dlogitnorm(x,mu=popt$mu,sigma=popt$sigma) ~ x, type="l" ) + abline( v=qlogitnorm(c(0.025,0.5,0.975),mu=popt$mu,sigma=popt$sigma)) + + lines( dlogitnorm(x,mu=popt2$mu,sigma=popt2$sigma) ~ x, type="l", col="maroon" ) + abline( v=qlogitnorm(c(0.5,0.9995),mu=popt2$mu,sigma=popt2$sigma), col="maroon") + + popt <- as.list(coefLogitnorm(c(0.9, 0.9995), perc=c(0.5,0.9995))) + plot( dlogitnorm(x,mu=popt$mu,sigma=popt$sigma) ~ x, type="l" ) + abline( v=qlogitnorm(c(0.025,0.5,0.975),mu=popt$mu,sigma=popt$sigma)) + +} + +.tmp.f <- function(){ + #visualize the objective functions surface + quant = c(0.6,0.9) + perc = c(0.5,0.975) + perc = c(0.5,0.9995) + perc=c(0.5,upperBoundProb) + quant = parms.var[varDist=="logitnorm",c("qMedian","qUpper")] + quant = parms.var["epsF",c("qMedian","qUpper")] + quant = parms.var["epsG",c("qMedian","qUpper")] + quant = parms.var["epsP",c("qMedian","qUpper")] + + + tmp.n <- 80 + tmp.mu <- seq(0,1,length.out=tmp.n) + tmp.sigma <- seq(0.01,2.5,length.out=tmp.n) + tmp <- as.matrix(expand.grid( mu=tmp.mu, sigma=tmp.sigma)) + tmp.of <- apply(tmp,1,ofLogitnorm, quant=quant,perc=perc ) + tmp.ofm <- matrix(tmp.of, nrow=tmp.n ) + image(tmp.mu, tmp.sigma, tmp.ofm) + image(tmp.mu, tmp.sigma, exp(-0.5*tmp.ofm)) #very flat + image(tmp.mu, tmp.sigma, exp(-0.5*1/(1/800)*tmp.ofm)) #distort by decreasing Temp *1/T + tmp.o <- coefLogitnorm(quant=quant, perc=perc, returnDetails=TRUE) + tmp.o + popt <- as.list(tmp.o$par) + points( popt$mu, popt$sigma) + + windows() + unlist(popt) + plot( dlogitnorm(x,mu=popt$mu,sigma=popt$sigma) ~ x, type="l" ) + abline( v=qlogitnorm(perc,mu=popt$mu,sigma=popt$sigma)) + lines( dlogitnorm(x,mu=popt$mu,sigma=1) ~ x, type="l", col="maroon" ) + abline( v=qlogitnorm(c(0.025,0.5,0.975),mu=popt$mu,sigma=1), col="maroon") +} + diff --git a/inst/unitTests/runitdlogitnormLog.R b/inst/unitTests/runitdlogitnormLog.R index caef2c9..e5ec086 100644 --- a/inst/unitTests/runitdlogitnormLog.R +++ b/inst/unitTests/runitdlogitnormLog.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.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)))) +} + diff --git a/man/dlogitnorm.Rd b/man/dlogitnorm.Rd index 1a3ceee..033bc1d 100755 --- a/man/dlogitnorm.Rd +++ b/man/dlogitnorm.Rd @@ -1,29 +1,29 @@ -\name{dlogitnorm} -\alias{dlogitnorm} -\title{dlogitnorm} -\description{Density function of logitnormal distribution } -\usage{dlogitnorm(q, mu = 0, sigma = 1, log = FALSE, ...)} -\arguments{ - \item{q}{quantiles} - \item{mu}{distribution parameters} - \item{sigma}{ -} - \item{log}{if TRUE, the log-density is returned} - \item{\dots}{further arguments passed to \code{\link{dnorm}}: \code{mean}, and \code{sd} for mu and sigma respectively.} -} -\details{\describe{\item{Logitnorm distribution}{ -\itemize{ -\item{density function: dlogitnorm } -\item{distribution function: \code{\link{plogitnorm}} } -\item{quantile function: \code{\link{qlogitnorm}} } -\item{random generation function: \code{\link{rlogitnorm}} } -} -}}} - - -\author{Thomas Wutzler} - - - -\seealso{\code{\link{logitnorm}}} - +\name{dlogitnorm} +\alias{dlogitnorm} +\title{dlogitnorm} +\description{Density function of logitnormal distribution } +\usage{dlogitnorm(q, mu = 0, sigma = 1, log = FALSE, ...)} +\arguments{ + \item{q}{quantiles} + \item{mu}{distribution parameters} + \item{sigma}{ +} + \item{log}{if TRUE, the log-density is returned} + \item{\dots}{further arguments passed to \code{\link{dnorm}}: \code{mean}, and \code{sd} for mu and sigma respectively.} +} +\details{\describe{\item{Logitnorm distribution}{ +\itemize{ +\item{density function: dlogitnorm } +\item{distribution function: \code{\link{plogitnorm}} } +\item{quantile function: \code{\link{qlogitnorm}} } +\item{random generation function: \code{\link{rlogitnorm}} } +} +}}} + + +\author{Thomas Wutzler} + + + +\seealso{\code{\link{logitnorm}}} + diff --git a/man/logitnorm-package.Rd b/man/logitnorm-package.Rd index f17c45e..684d840 100755 --- a/man/logitnorm-package.Rd +++ b/man/logitnorm-package.Rd @@ -1,57 +1,57 @@ -\name{logitnorm-package} -\alias{logitnorm-package} \alias{logitnorm} -\title{Utilities for the logitnormal distribution in R} -\description{ Utilities for the logitnormal distribution in R -\itemize{ -\item Density, distribution, quantile and random generation function. -\item Estimation of the mode and the first two moments. -\item Estimation of distribution parameters from observations. -} -}%description - - -\details{ -The logitnormal distribution is useful as a prior density for variables that are bounded between 0 and 1, -such as proportions. Fig. 1 displays its density for various combinations of parameters mu and sigma. - -The package provides the main distribution functions: -\itemize{ -\item density \code{\link{dlogitnorm}}, -\item distribution \code{\link{plogitnorm}}, -\item quantile \code{\link{qlogitnorm}}, and -\item random generation function \code{\link{rlogitnorm}}. -}%itemize - -Transformation functions - \itemize{ - \item{ (0,1) -> (-Inf,Inf): \code{\link{logit}} } - \item{ (-Inf,Inf) -> (0,1): \code{\link{invlogit}} } - } - -Moments and mode - \itemize{ - \item{ Expected value and variance: \code{\link{momentsLogitnorm}} } - \item{ Mode: \code{\link{modeLogitnorm}} } - } - -Estimating parameters - \itemize{ - \item{from mode and upper quantile: \code{\link{twCoefLogitnormMLE}} } - \item{from mode and constraint to be unimodal and maximally flat: \code{\link{twCoefLogitnormMLEFlat}} } - \item{from median and upper quantile: \code{\link{twCoefLogitnorm}} } - \item{from expected value, i.e. mean and upper quantile: \code{\link{twCoefLogitnormE}} } - \item{from a confidence interval which is symmetric at normal scale: \code{\link{twCoefLogitnormCi}} } - \item{from prescribed quantiles: \code{\link{twCoefLogitnormN}} } - } - - }%details - -\references{ -Frederic, P. & Lad, F. (2008) Two Moments of the Logitnormal Distribution. -Communications in Statistics-Simulation and Computation, 37, 1263-1269 -}%references -\author{Thomas Wutzler} - -\keyword{ package } - - +\name{logitnorm-package} +\alias{logitnorm-package} \alias{logitnorm} +\title{Utilities for the logitnormal distribution in R} +\description{ Utilities for the logitnormal distribution in R +\itemize{ +\item Density, distribution, quantile and random generation function. +\item Estimation of the mode and the first two moments. +\item Estimation of distribution parameters from observations. +} +}%description + + +\details{ +The logitnormal distribution is useful as a prior density for variables that are bounded between 0 and 1, +such as proportions. Fig. 1 displays its density for various combinations of parameters mu and sigma. + +The package provides the main distribution functions: +\itemize{ +\item density \code{\link{dlogitnorm}}, +\item distribution \code{\link{plogitnorm}}, +\item quantile \code{\link{qlogitnorm}}, and +\item random generation function \code{\link{rlogitnorm}}. +}%itemize + +Transformation functions + \itemize{ + \item{ (0,1) -> (-Inf,Inf): \code{\link{logit}} } + \item{ (-Inf,Inf) -> (0,1): \code{\link{invlogit}} } + } + +Moments and mode + \itemize{ + \item{ Expected value and variance: \code{\link{momentsLogitnorm}} } + \item{ Mode: \code{\link{modeLogitnorm}} } + } + +Estimating parameters + \itemize{ + \item{from mode and upper quantile: \code{\link{twCoefLogitnormMLE}} } + \item{from mode and constraint to be unimodal and maximally flat: \code{\link{twCoefLogitnormMLEFlat}} } + \item{from median and upper quantile: \code{\link{twCoefLogitnorm}} } + \item{from expected value, i.e. mean and upper quantile: \code{\link{twCoefLogitnormE}} } + \item{from a confidence interval which is symmetric at normal scale: \code{\link{twCoefLogitnormCi}} } + \item{from prescribed quantiles: \code{\link{twCoefLogitnormN}} } + } + + }%details + +\references{ +Frederic, P. & Lad, F. (2008) Two Moments of the Logitnormal Distribution. +Communications in Statistics-Simulation and Computation, 37, 1263-1269 +}%references +\author{Thomas Wutzler} + +\keyword{ package } + + diff --git a/man/twCoefLogitnormMLE.Rd b/man/twCoefLogitnormMLE.Rd index 547c5a4..1e41405 100755 --- a/man/twCoefLogitnormMLE.Rd +++ b/man/twCoefLogitnormMLE.Rd @@ -1,36 +1,36 @@ -\name{twCoefLogitnormMLE} -\alias{twCoefLogitnormMLE} -\title{twCoefLogitnormMLE} -\description{Estimating coefficients of logitnormal distribution from mode and upper quantile } -\usage{twCoefLogitnormMLE(mle, quant, perc = 0.999)} -\arguments{ - \item{mle}{numeric vector: the mode of the density function} - \item{quant}{numeric vector: the upper quantile value} - \item{perc}{numeric vector: the probability for which the quantile was specified} -} - -\value{numeric matrix with columns \code{c("mu","sigma")} -rows correspond to rows in mle, quant, and perc} - -\author{Thomas Wutzler} - - - -\seealso{\code{\link{logitnorm}}} -\examples{ - -# estimate the parameters, with mode 0.7 and upper quantile 0.9 -(theta <- twCoefLogitnormMLE(0.7,0.9)) - -x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid -px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function -plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.999),col="gray") -dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function -plot(dx~x); abline(v=c(0.7,0.9),col="gray") - -# vectorized -(theta <- twCoefLogitnormMLE(mle=seq(0.4,0.8,by=0.1),quant=0.9)) - - # flat - (theta <- twCoefLogitnormMLEFlat(0.7)) -} +\name{twCoefLogitnormMLE} +\alias{twCoefLogitnormMLE} +\title{twCoefLogitnormMLE} +\description{Estimating coefficients of logitnormal distribution from mode and upper quantile } +\usage{twCoefLogitnormMLE(mle, quant, perc = 0.999)} +\arguments{ + \item{mle}{numeric vector: the mode of the density function} + \item{quant}{numeric vector: the upper quantile value} + \item{perc}{numeric vector: the probability for which the quantile was specified} +} + +\value{numeric matrix with columns \code{c("mu","sigma")} +rows correspond to rows in mle, quant, and perc} + +\author{Thomas Wutzler} + + + +\seealso{\code{\link{logitnorm}}} +\examples{ + +# estimate the parameters, with mode 0.7 and upper quantile 0.9 +(theta <- twCoefLogitnormMLE(0.7,0.9)) + +x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid +px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function +plot(px~x); abline(v=c(0.7,0.9),col="gray"); abline(h=c(0.999),col="gray") +dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function +plot(dx~x); abline(v=c(0.7,0.9),col="gray") + +# vectorized +(theta <- twCoefLogitnormMLE(mle=seq(0.4,0.8,by=0.1),quant=0.9)) + + # flat + (theta <- twCoefLogitnormMLEFlat(0.7)) +} diff --git a/man/twCoefLogitnormMLEFlat.Rd b/man/twCoefLogitnormMLEFlat.Rd index 6e8fff3..b69135d 100644 --- a/man/twCoefLogitnormMLEFlat.Rd +++ b/man/twCoefLogitnormMLEFlat.Rd @@ -1,20 +1,20 @@ -\name{twCoefLogitnormMLEFlat} -\alias{twCoefLogitnormMLEFlat} -\title{twCoefLogitnormMLEFlat} -\description{Estimating coefficients of a maximally flat unimodal logitnormal distribution from 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 -densitiy to a maximum range, i.e. is maximally flat, but still is unimodal.} - - -\author{Thomas Wutzler} - - - - - +\name{twCoefLogitnormMLEFlat} +\alias{twCoefLogitnormMLEFlat} +\title{twCoefLogitnormMLEFlat} +\description{Estimating coefficients of a maximally flat unimodal logitnormal distribution from 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 +densitiy to a maximum range, i.e. is maximally flat, but still is unimodal.} + + +\author{Thomas Wutzler} + + + + + diff --git a/man/twSigmaLogitnorm.Rd b/man/twSigmaLogitnorm.Rd index d262b6f..9db89e0 100755 --- a/man/twSigmaLogitnorm.Rd +++ b/man/twSigmaLogitnorm.Rd @@ -1,30 +1,30 @@ -\name{twSigmaLogitnorm} -\alias{twSigmaLogitnorm} -\title{twSigmaLogitnorm} -\description{Estimating coefficients of logitnormal distribution from mode and given mu } -\usage{twSigmaLogitnorm(mle, mu = 0)} -\arguments{ - \item{mle}{numeric vector: the mode of the density function} - \item{mu}{for mu=0 the distribution will be the flattest case (maybe bimodal)} -} -\details{For a mostly flat unimodal distribution use \code{\link{twCoefLogitnormMLE}(mle,0)}} -\value{numeric matrix with columns \code{c("mu","sigma")} -rows correspond to rows in mle and mu} - -\author{Thomas Wutzler} - - - -\seealso{\code{\link{logitnorm}}} -\examples{ - mle <- 0.8 - (theta <- twSigmaLogitnorm(mle)) - # -x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid -px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function -plot(px~x); abline(v=c(mle),col="gray") -dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function -plot(dx~x); abline(v=c(mle),col="gray") -# vectorized -(theta <- twSigmaLogitnorm(mle=seq(0.401,0.8,by=0.1))) -} +\name{twSigmaLogitnorm} +\alias{twSigmaLogitnorm} +\title{twSigmaLogitnorm} +\description{Estimating coefficients of logitnormal distribution from mode and given mu } +\usage{twSigmaLogitnorm(mle, mu = 0)} +\arguments{ + \item{mle}{numeric vector: the mode of the density function} + \item{mu}{for mu=0 the distribution will be the flattest case (maybe bimodal)} +} +\details{For a mostly flat unimodal distribution use \code{\link{twCoefLogitnormMLE}(mle,0)}} +\value{numeric matrix with columns \code{c("mu","sigma")} +rows correspond to rows in mle and mu} + +\author{Thomas Wutzler} + + + +\seealso{\code{\link{logitnorm}}} +\examples{ + mle <- 0.8 + (theta <- twSigmaLogitnorm(mle)) + # +x <- seq(0,1,length.out=41)[-c(1,41)] # plotting grid +px <- plogitnorm(x,mu=theta[1],sigma=theta[2]) #percentiles function +plot(px~x); abline(v=c(mle),col="gray") +dx <- dlogitnorm(x,mu=theta[1],sigma=theta[2]) #density function +plot(dx~x); abline(v=c(mle),col="gray") +# vectorized +(theta <- twSigmaLogitnorm(mle=seq(0.401,0.8,by=0.1))) +} diff --git a/tests/doRUnit.R b/tests/doRUnit.R index 2690e6b..38d478c 100644 --- a/tests/doRUnit.R +++ b/tests/doRUnit.R @@ -1,60 +1,63 @@ ## unit tests will not be done if RUnit is not available -if(require("RUnit", quietly=TRUE)) { +if (require("RUnit", quietly = TRUE)) { ## --- Setup --- #pkg <- "PKG" # <-- Change to package name! pkg <- "logitnorm" # <-- Change to package name! - if(Sys.getenv("RCMDCHECK") == "FALSE") { + if (Sys.getenv("RCMDCHECK") == "FALSE") { ## Path to unit tests for standalone running under Makefile (not R CMD check) ## PKG/tests/../inst/unitTests path <- file.path(getwd(), "..", "inst", "unitTests") } else { ## Path to unit tests for R CMD check ## PKG.Rcheck/tests/../PKG/unitTests - path <- system.file(package=pkg, "unitTests") + path <- system.file(package = pkg, "unitTests") } cat("\nRunning unit tests\n") - print(list(pkg=pkg, getwd=getwd(), pathToUnitTests=path)) + print(list(pkg = pkg, getwd = getwd(), pathToUnitTests = path)) - library(package=pkg, character.only=TRUE) + library(package = pkg, character.only = TRUE) ## If desired, load the name space to allow testing of private functions ## if (is.element(pkg, loadedNamespaces())) - ## attach(loadNamespace(pkg), name=paste("namespace", pkg, sep=":"), pos=3) + ## attach(loadNamespace(pkg), name = paste("namespace", pkg, sep = ":"), pos = 3) ## ## or simply call PKG:::myPrivateFunction() in tests ## --- Testing --- ## Define tests - testSuite <- defineTestSuite(name=paste(pkg, "unit testing"), - dirs=path) + testSuite <- defineTestSuite(name = paste(pkg, "unit testing"), + dirs = path) ## Run tests <- runTestSuite(testSuite) - ## Default report name - pathReport <- file.path(path, "report") - - ## Report to stdout and text files - cat("------------------- UNIT TEST SUMMARY ---------------------\n\n") - printTextProtocol(tests, showDetails=FALSE) - printTextProtocol(tests, showDetails=FALSE, - fileName=paste(pathReport, "Summary.txt", sep="")) - printTextProtocol(tests, showDetails=TRUE, - fileName=paste(pathReport, ".txt", sep="")) - - ## Report to HTML file - printHTMLProtocol(tests, fileName=paste(pathReport, ".html", sep="")) - ## Return stop() to cause R CMD check stop in case of ## - failures i.e. FALSE to unit tests or ## - errors i.e. R errors tmp <- getErrors(tests) - if(tmp$nFail > 0 | tmp$nErr > 0) { + if (tmp$nFail > 0 | tmp$nErr > 0) { stop(paste("\n\nunit testing failed (#test failures: ", tmp$nFail, - ", #R errors: ", tmp$nErr, ")\n\n", sep="")) + ", #R errors: ", tmp$nErr, ")\n\n", sep = "")) } } else { warning("cannot run unit tests -- package RUnit is not available") -} \ No newline at end of file +} + +.tmp.f <- function(){ + #twutz 1712 removed writing test protocols, as its against CRAN policy + ## Default report name + pathReport <- file.path(path, "report") + + ## Report to stdout and text files + cat("------------------- UNIT TEST SUMMARY ---------------------\n\n") + printTextProtocol(tests, showDetails = FALSE) + printTextProtocol(tests, showDetails = FALSE, + fileName = paste(pathReport, "Summary.txt", sep = "")) + printTextProtocol(tests, showDetails = TRUE, + fileName = paste(pathReport, ".txt", sep = "")) + + ## Report to HTML file + printHTMLProtocol(tests, fileName = paste(pathReport, ".html", sep = "")) +} diff --git a/vignettes/logitnorm.Rmd b/vignettes/logitnorm.Rmd index 23b954e..cdfc869 100644 --- a/vignettes/logitnorm.Rmd +++ b/vignettes/logitnorm.Rmd @@ -1,134 +1,134 @@ - - -```{r setup, include=FALSE} -library(knitr) -opts_chunk$set(out.extra='style="display:block; margin: auto"' - #, fig.align="center" - #, fig.width=4.6, fig.height=3.2 - , fig.width=6, fig.height=3.75 #goldener Schnitt 1.6 - , dev.args=list(pointsize=10) - , dev=c('png','pdf') - ) -knit_hooks$set(spar = function(before, options, envir) { - if (before){ - par( las=1 ) #also y axis labels horizontal - par(mar=c(2.0,3.3,0,0)+0.3 ) #margins - par(tck=0.02 ) #axe-tick length inside plots - par(mgp=c(1.1,0.2,0) ) #positioning of axis title, axis labels, axis - } -}) -library(logitnorm) -if( !require(ggplot2) ){ - print("To generate this vignette, ggplo2 is required.") - exit(0) -} -library(reshape2) -# genVigs("logitnorm") -``` - -Basic usage -=============== - -Distribution --------------- - -The logitnormal distribution is useful as a prior density for variables that are -bounded between 0 and 1, such as proportions. The following figure -displays its density for various combinations of parameters mu (panels) and -sigma (lines). - - -```{r densityPlots, spar=TRUE, echo=FALSE, fig.width=6*2, fig.height=3.75*2} -xGrid =c(0+.Machine$double.eps,seq(0,1, length.out=81)[-c(1,81)],1-.Machine$double.eps) -#explore chanign sigma at mu=0 -theta0 <- expand.grid(mu=seq(0,2,length.out=9), -sigma=10^seq(-0.5,0.5,length.out=5)) -n <- nrow(theta0) - -.calcDensityGrid <- function( - ### Calculate lognormal desnity for given combinations - theta0 ##<< matrix with columns mu and sigma - ,xGrid = seq(0,1, length.out=81)[-c(1,81)] -){ - dx <- apply( theta0, 1, function(theta0i){ - dx <- dlogitnorm(xGrid, mu=theta0i[1], sigma=theta0i[2]) - }) - dimnames(dx) <- list(iX=NULL,iTheta=NULL) - ds <- melt(dx) # two variables over time - ds[1:10,] - ds$mu <- rep(as.factor(round(theta0[,1],2)), each=length(xGrid)) - ds$sigma <- rep(as.factor(round(theta0[,2],2)), each=length(xGrid)) - ds$x <- rep(xGrid, nrow(theta0)) - ds - ### data frame with columns value,x,mu and sigma -} - -ds <- .calcDensityGrid(theta0,xGrid=xGrid) -p1 <- ggplot(ds, aes(x=x, y=value, color=sigma, linetype=sigma)) + geom_line(size=1) + - facet_wrap(~mu,scales="free")+ - theme_bw() + - theme(axis.title.x = element_blank()) + ylab("density") + - theme() -print(p1) -``` - -Example: Plot the cumulative distribution -```{r cumDensityPlot, spar=TRUE} -x <- seq(0,1, length.out=81) -d <- plogitnorm(x, mu=0.5, sigma=0.5) -plot(d~x,type="l") -``` - -Mean and Variance -------------------- -The moments have no analytical solution. This package -estimates them by numerical integration: - -Example: estimate mean and standard deviation. -```{r momentsLogitnorm} -(theta <- momentsLogitnorm(mu=0.6,sigma=0.5)) -``` - -Mode ----- -The mode is found by setting derivatives to zero and optimizing -the resulting equation: $logit(x) = \sigma^2(2x-1)+\mu$. - -Example: estimate the mode -```{r modeLogitnorm} -(mle <- modeLogitnorm(mu=0.6,sigma=0.5)) -``` - -Parameter Estimation ------------------------ -from upper quantile and - - mode (Maximum Likelihood Estimate) - - mean (Expected value) - - median - -Example: estimate the parameters, with mode 0.7 and upper quantile 0.9 -```{r twCoefLogitnormMLE} -(theta <- twCoefLogitnormMLE(0.7,0.9)) -x <- seq(0,1, length.out=81) -d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) -plot(d~x,type="l") -abline(v=c(0.7,0.9), col="grey") -``` - -When increasing the $\sigma$ parameter, the distribution becomes -eventually becomes bi-model, i.e. has two maxima. The unimodal distribution for -a given mode with widest confidence intervals is obtained by -function `twCoefLogitnormMLEFlat`. - -```{r twCoefLogitnormMLEFlat} -(theta <- twCoefLogitnormMLEFlat(0.7)) -x <- seq(0,1, length.out=81) -d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) -plot(d~x,type="l") -abline(v=c(0.7), col="grey") -``` - - + + +```{r setup, include=FALSE} +library(knitr) +opts_chunk$set(out.extra='style="display:block; margin: auto"' + #, fig.align="center" + #, fig.width=4.6, fig.height=3.2 + , fig.width=6, fig.height=3.75 #goldener Schnitt 1.6 + , dev.args=list(pointsize=10) + , dev=c('png','pdf') + ) +knit_hooks$set(spar = function(before, options, envir) { + if (before){ + par( las=1 ) #also y axis labels horizontal + par(mar=c(2.0,3.3,0,0)+0.3 ) #margins + par(tck=0.02 ) #axe-tick length inside plots + par(mgp=c(1.1,0.2,0) ) #positioning of axis title, axis labels, axis + } +}) +library(logitnorm) +if( !require(ggplot2) ){ + print("To generate this vignette, ggplo2 is required.") + exit(0) +} +library(reshape2) +# genVigs("logitnorm") +``` + +Basic usage +=============== + +Distribution +-------------- + +The logitnormal distribution is useful as a prior density for variables that are +bounded between 0 and 1, such as proportions. The following figure +displays its density for various combinations of parameters mu (panels) and +sigma (lines). + + +```{r densityPlots, spar=TRUE, echo=FALSE, fig.width=6*2, fig.height=3.75*2} +xGrid =c(0+.Machine$double.eps,seq(0,1, length.out=81)[-c(1,81)],1-.Machine$double.eps) +#explore chanign sigma at mu=0 +theta0 <- expand.grid(mu=seq(0,2,length.out=9), +sigma=10^seq(-0.5,0.5,length.out=5)) +n <- nrow(theta0) + +.calcDensityGrid <- function( + ### Calculate lognormal desnity for given combinations + theta0 ##<< matrix with columns mu and sigma + ,xGrid = seq(0,1, length.out=81)[-c(1,81)] +){ + dx <- apply( theta0, 1, function(theta0i){ + dx <- dlogitnorm(xGrid, mu=theta0i[1], sigma=theta0i[2]) + }) + dimnames(dx) <- list(iX=NULL,iTheta=NULL) + ds <- melt(dx) # two variables over time + ds[1:10,] + ds$mu <- rep(as.factor(round(theta0[,1],2)), each=length(xGrid)) + ds$sigma <- rep(as.factor(round(theta0[,2],2)), each=length(xGrid)) + ds$x <- rep(xGrid, nrow(theta0)) + ds + ### data frame with columns value,x,mu and sigma +} + +ds <- .calcDensityGrid(theta0,xGrid=xGrid) +p1 <- ggplot(ds, aes(x=x, y=value, color=sigma, linetype=sigma)) + geom_line(size=1) + + facet_wrap(~mu,scales="free")+ + theme_bw() + + theme(axis.title.x = element_blank()) + ylab("density") + + theme() +print(p1) +``` + +Example: Plot the cumulative distribution +```{r cumDensityPlot, spar=TRUE} +x <- seq(0,1, length.out=81) +d <- plogitnorm(x, mu=0.5, sigma=0.5) +plot(d~x,type="l") +``` + +Mean and Variance +------------------- +The moments have no analytical solution. This package +estimates them by numerical integration: + +Example: estimate mean and standard deviation. +```{r momentsLogitnorm} +(theta <- momentsLogitnorm(mu=0.6,sigma=0.5)) +``` + +Mode +---- +The mode is found by setting derivatives to zero and optimizing +the resulting equation: $logit(x) = \sigma^2(2x-1)+\mu$. + +Example: estimate the mode +```{r modeLogitnorm} +(mle <- modeLogitnorm(mu=0.6,sigma=0.5)) +``` + +Parameter Estimation +----------------------- +from upper quantile and + - mode (Maximum Likelihood Estimate) + - mean (Expected value) + - median + +Example: estimate the parameters, with mode 0.7 and upper quantile 0.9 +```{r twCoefLogitnormMLE} +(theta <- twCoefLogitnormMLE(0.7,0.9)) +x <- seq(0,1, length.out=81) +d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) +plot(d~x,type="l") +abline(v=c(0.7,0.9), col="grey") +``` + +When increasing the $\sigma$ parameter, the distribution becomes +eventually becomes bi-model, i.e. has two maxima. The unimodal distribution for +a given mode with widest confidence intervals is obtained by +function `twCoefLogitnormMLEFlat`. + +```{r twCoefLogitnormMLEFlat} +(theta <- twCoefLogitnormMLEFlat(0.7)) +x <- seq(0,1, length.out=81) +d <- dlogitnorm(x, mu=theta[1,"mu"], sigma=theta[1,"sigma"]) +plot(d~x,type="l") +abline(v=c(0.7), col="grey") +``` + +