R/logitnorm.R

Defines functions twSigmaLogitnorm modeLogitnorm ofModeLogitnorm momentsLogitnorm twCoefLogitnormE ofLogitnormE twCoefLogitnormMLEFlat twCoefLogitnormMLE ofLogitnormMLE twCoefLogitnormCi twCoefLogitnorm twCoefLogitnormN ofLogitnorm qlogitnorm dlogitnorm plogitnorm rlogitnorm invlogit logit

Documented in dlogitnorm invlogit logit modeLogitnorm momentsLogitnorm plogitnorm qlogitnorm rlogitnorm twCoefLogitnorm twCoefLogitnormCi twCoefLogitnormE twCoefLogitnormMLE twCoefLogitnormMLEFlat twCoefLogitnormN twSigmaLogitnorm

# 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)<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}}
){
	##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)))
}

Try the logitnorm package in your browser

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

logitnorm documentation built on May 31, 2017, 4:45 a.m.