# 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
){
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{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 is mu, theta 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 <= 0) return( Inf )
tmp.predp = plogitnorm(quant, mu = theta, sigma = theta )
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, sigma = theta )
#tmp.diff =  tmp.predq - quant
sum(tmp.diff^2)
}

twCoefLogitnormN <- function(
### Estimating coefficients from a vector of quantiles and perentiles (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,sigma = theta)
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
, ...
){
# 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,"; quant = ",qmatI,"; perc = ",pmatI,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,sigma = theta)	#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,sigma = theta)	#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
){
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)<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]
if (mleI == 0.5) mleI = 0.5-.Machine$double.eps # 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 # 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]), perc = (percI<-perc[1+i0%%nc]))
#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 else
mu <- optimize( .ofLogitnormMLE, interval = intv, mle = (mleI), logitMle = logitMle, quant = (quantI<-quant[1+i0%%nc]), perc = (percI<-perc[1+i0%%nc]))$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 (theta <- twCoefLogitnormMLE(0.7,0.9)) x <- seq(0,1,length.out = 41)[-c(1,41)] # plotting grid px <- plogitnorm(x,mu = theta,sigma = theta) #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,sigma = theta) #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 ## 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 is the mu, theta 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, sigma = theta ) tmp.diff = tmp.predp - perc .exp <- momentsLogitnorm(theta,theta,...)["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 , ... ){ ##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, quant = mqpi, perc = mqpi, 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,sigma = thetaE) #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,sigma = thetaE) #density function plot(dx~x); abline(v = c(0.7,0.9),col = "gray") z <- rlogitnorm(1e5, mu = thetaE,sigma = thetaE) 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
){
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,sigma = theta)	#percentiles function
plot(px~x); abline(v = c(mle),col = "gray")
dx <- dlogitnorm(x,mu = theta,sigma = theta)	#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 May 1, 2019, 11:01 p.m.