# R/logitnorm.R In logitnorm: Functions for the Logitnormal Distribution

#### Documented in dlogitnorminvlogitlogitmodeLogitnormmomentsLogitnormplogitnormqlogitnormrlogitnormtwCoefLogitnormtwCoefLogitnormCitwCoefLogitnormEtwCoefLogitnormMLEtwCoefLogitnormMLEFlattwCoefLogitnormNtwSigmaLogitnorm

# 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)}
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}} \!}
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
n,          ##<< number of observations
mu = 0,     ##<< distribution parameter
sigma = 1,	##<< distribution parameter
...	##<< arguments to \code{\link{rnorm}}
){
plogis( rnorm(n = n, mean = mu, sd = sigma, ...) )
}

plogitnorm <- function(
### Distribution function for logitnormal distribution
q  ##<< vector of quantiles
, mu = 0, sigma = 1	##<< distribution parameters
, ...
){
ql <- qlogis(q)
pnorm(ql,mean = mu,sd = sigma,...)
}

dlogitnorm <- function(
### Density function of logitnormal distribution
x		##<< vector of 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}} }
## }
## }}
##details<< The function is only defined in interval (0,1), but the density
## returns 0 outside the support region.
ql <- qlogis(x)
# multiply (or add in the log domain) by the Jacobian (derivative) of the
# back-Transformation (logit)
if (log) {
ifelse(
x <= 0 | x >= 1, 0,
dnorm(ql, mean = mu, sd = sigma, log = TRUE, ...) - log(x) - log1p(-x)
)
} else {
ifelse(
x <= 0 | x >= 1, 0,
dnorm(ql,mean = mu,sd = sigma,...)/x/(1 - x)
)
}
}

qlogitnorm <- function(
### Quantiles of logitnormal distribution.
p ##<< vector of probabilities
, mu = 0, sigma = 1	##<< distribution parameters
, ...
){
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
#
# 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 from a vector of quantiles and percentiles (non-vectorized).
quant					##<< the quantile values
,perc = c(0.5,0.975)		##<< the probabilities 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. \code{control = list(maxit = 1000)}
){
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 distribution.
### 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
, ...
){
# twCoefLogitnorm
mu = logit(median)
upperLogit = logit(quant)
sigmaFac = qnorm(perc)
sigma = 1/sigmaFac*(upperLogit - mu)
c(mu = mu, sigma = sigma)
### 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
med = 0.7; upper = 0.9
med = 0.2; upper = 0.4
(theta <- twCoefLogitnorm(med,upper))

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(med,upper),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(med,upper),col = "gray")

# vectorized
(theta <- twCoefLogitnorm(seq(0.4,0.8,by = 0.1),0.9))

.tmp.f <- function(){
# xr = rlogitnorm(1e5, mu = theta["mu"], sigma = theta["sigma"])
# median(xr)
invlogit(theta["mu"])
qlogitnorm(0.975, mu = theta["mu"], sigma = theta["sigma"])
}
}

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
){
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, perc = p)
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)<mle for mle>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}}
){
# 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]]
quantI <- quant[1 + i0 %% nc[2]]
percI <- perc[1 + i0 %% nc[3]]
if (mleI == 0.5) {
# if mode is at 0.5 its symmetrical and corresponds to median
mu = 0
res[i,] = twCoefLogitnorm(0.5, quantI, percI)
break()
}
# in order to estimate sigma
# we now that mu is in (\code{logit(mle)},0) for \code{mle < 0.5} and in
# \code{(0,logit(mle))} for \code{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 that is spaced narrower at the edge
#
#intv <- if (logitMle < 0) c(logitMle+.Machine$double.eps,0) else c(0,logitMle-.Machine$double.eps)
logitMle <- logit(mleI)
upper <- abs(logitMle) - .Machine$double.eps muTry <- sign(logitMle)*log(seq(1,exp(upper),length.out = 40)) ofMuTry <- .ofLogitnormMLE( muTry, mle = (mleI), logitMle = logitMle , quant = (quantI <- quant[1 + i0 %% nc[2]]) , perc = (percI <- perc[1 + i0 %% nc[3]])) #plot(muTry,ofMuTry) # now optimize starting from the near minimum imin <- which.min(ofMuTry) intv <- muTry[ c(max(1,imin - 1), min(length(muTry),imin + 1)) ] if (diff(intv) == 0 ) mu <- intv[1] else mu <- optimize( .ofLogitnormMLE, interval = intv, mle = (mleI), logitMle = logitMle , quant = quantI , perc = percI)$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 \code{mle}, \code{quant}, and \code{perc}
}
#mtrace(coefLogitnorm)
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
plot(px~x); abline(v = c(mode,upper),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(mode,upper),col = "gray")
# vectorized
(theta <- twCoefLogitnormMLE(mle = seq(0.4,0.8,by = 0.1),quant = upper))
# flat
(theta <- twCoefLogitnormMLEFlat(mode))
}

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
## density 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 probabilities 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 attribute resOptim
, ...
){
# 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
### distribution.
### 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	##<< changing 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) )
# 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)))
}


## Try the logitnorm package in your browser

Any scripts or data that you put into this service are public.

logitnorm documentation built on Jan. 7, 2021, 9:09 a.m.