Density returns 0 instead of NAN outside (0,1) (#4)

closes #2

NEWS.md

@@ -1,3 +1,7 @@
+# logitnorm 0.8.37
+
+- Density returns 0 instead of NAN outside (0,1)
+
 # logitnorm 0.8.36
 Remove the library call to MASS.
 

R/logitnorm.R

@@ -66,13 +66,21 @@ dlogitnorm <- function(
   ## }
   ## }}
   ##seealso<< \code{\link{logitnorm}}
+  ##details<< The function is only defined in interval (0,1), but the density 
+  ## returns 0 outside the support region.
   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)
+    ifelse(
+      q <= 0 | q >= 1, 0, 
+      dnorm(ql, mean = mu, sd = sigma, log = TRUE, ...) - log(q) - log1p(-q)
+    )
   } else {
-    dnorm(ql,mean = mu,sd = sigma,...)/q/(1 - q)
+    ifelse(
+      q <= 0 | q >= 1, 0, 
+      dnorm(ql,mean = mu,sd = sigma,...)/q/(1 - q)
+    )
   }
 }
 

inst/unitTests/runitLogitnorm.R

@@ -1,3 +1,7 @@
+.tmp.f <- function(){
+  library(RUnit)
+}
+
 .setUp <- function() {
 	#library(MASS)
 	.setUpDf <- within( list(),{
@@ -32,6 +36,13 @@ test.plogitnorm <- function(){
 	#lines( pnorm( logit(x)) ~ logit(x) )
 }
 
+test.dlogitnorm <- function(){
+  q <- c(-1,0,0.5,1,2)
+  ans <- suppressWarnings(dlogitnorm(q))
+  checkEquals(c(0,0,1.595769,0,0), ans, tolerance = 1e-7)
+}
+
+
 test.twCoefLogitnorm <- function(){
 	theta <- twCoefLogitnorm(0.7,0.9,perc = 0.999)
 	px <- plogitnorm(x,mu = theta[1],sigma = theta[2])	#percentiles function

man/dlogitnorm.Rd

@@ -10,7 +10,8 @@
   \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.}
+  \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{
@@ -19,7 +20,10 @@
 \item{quantile function: \code{\link{qlogitnorm}} }
 \item{random generation function: \code{\link{rlogitnorm}} }
 }
-}}}
+}}
+
+The function is only defined in interval (0,1), but the density 
+returns 0 outside the support region.}
 
 
 \author{Thomas Wutzler}

man/logit.Rd

@@ -9,7 +9,7 @@
   \item{\dots}{
 }
 }
-\details{function \eqn{ logit(p)= log \left( \frac{p}{1-p} \right) = log(p) - log(1-p) }}
+\details{function \eqn{logit(p) = log \left( \frac{p}{1-p} \right) = log(p) - log(1-p)}}
 
 
 \author{Thomas Wutzler}

man/twCoefLogitnorm.Rd

@@ -28,13 +28,13 @@ rows correspond to rows in median, quant, and perc}
 # 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")
+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")
+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))
+(theta <- twCoefLogitnorm(seq(0.4,0.8,by = 0.1),0.9))
 }

man/twCoefLogitnormCi.Rd

@@ -8,7 +8,7 @@
   \item{lower}{value at the lower quantile, i.e. practical minimum}
   \item{upper}{value at the upper quantile, i.e. practical maximum}
   \item{perc}{numeric vector: the probability for which the quantile was specified}
-  \item{sigmaFac}{sigmaFac=2 is 95\% sigmaFac=2.6 is 99\% interval}
+  \item{sigmaFac}{sigmaFac = 2 is 95\% sigmaFac = 2.6 is 99\% interval}
   \item{isTransScale}{if true lower and upper are already on logit scale}
 }
 
@@ -20,9 +20,9 @@
 
 \seealso{\code{\link{logitnorm}}}
 \examples{
-mu=2
-sd=c(1,0.8)
-p=0.99
+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)	

man/twCoefLogitnormE.Rd

@@ -8,7 +8,7 @@
 \arguments{
   \item{mean}{the expected value of the density function}
   \item{quant}{the quantile values}
-  \item{perc}{the probabilites for which the quantiles were specified}
+  \item{perc}{the probabilities for which the quantiles were specified}
   \item{method}{method of optimization (see \code{\link{optim}})}
   \item{theta0}{starting parameters}
   \item{returnDetails}{if TRUE, the full output of optim is returned with attribute resOptim}
@@ -28,15 +28,15 @@ colnames: \code{c("mu","sigma")}}
 # 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")
+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])
+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))
+(theta <- twCoefLogitnormE(mean = seq(0.4,0.8,by = 0.1),quant = 0.9))
 }

man/twCoefLogitnormMLE.Rd

@@ -10,7 +10,7 @@
 }
 
 \value{numeric matrix with columns \code{c("mu","sigma")}
-rows correspond to rows in mle, quant, and perc}
+rows correspond to rows in \code{mle}, \code{quant}, and \code{perc}}
 
 \author{Thomas Wutzler}
 
@@ -22,15 +22,15 @@ rows correspond to rows in mle, quant, and perc}
 # 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")
+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))
+(theta <- twCoefLogitnormMLE(mle = seq(0.4,0.8,by = 0.1),quant = 0.9))
+
+# flat
+(theta <- twCoefLogitnormMLEFlat(0.7))
 }

man/twCoefLogitnormN.Rd

@@ -1,20 +1,20 @@
 \name{twCoefLogitnormN}
 \alias{twCoefLogitnormN}
 \title{twCoefLogitnormN}
-\description{Estimating coefficients of logitnormal distribution from a vector of quantiles and perentiles (non-vectorized).	}
+\description{Estimating coefficients from a vector of quantiles and perentiles (non-vectorized).	}
 \usage{twCoefLogitnormN(quant, perc = c(0.5, 0.975), 
     method = "BFGS", theta0 = c(mu = 0, sigma = 1), 
     returnDetails = FALSE, ...)}
 \arguments{
   \item{quant}{the quantile values}
-  \item{perc}{the probabilites for which the quantiles were specified}
+  \item{perc}{the probabilities for which the quantiles were specified}
   \item{method}{method of optimization (see \code{\link{optim}})}
   \item{theta0}{starting parameters}
   \item{returnDetails}{if TRUE, the full output of optim is returned instead of only entry par}
-  \item{\dots}{further parameters passed to optim, e.g. control=list(maxit=1000)}
+  \item{\dots}{further parameters passed to optim, e.g. \code{control = list(maxit = 1000)}}
 }
 
-\value{named numeric vector with estimated parameters of the logitnormal distrubtion.
+\value{named numeric vector with estimated parameters of the logitnormal distribution.
 names: \code{c("mu","sigma")}}
 
 \author{Thomas Wutzler}
@@ -23,17 +23,17 @@ names: \code{c("mu","sigma")}}
 
 \seealso{\code{\link{logitnorm}}}
 \examples{
-    # 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))
+# 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))
 
-    # 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")
+# 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")
 }

man/twSigmaLogitnorm.Rd

@@ -5,7 +5,7 @@
 \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)}
+  \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")}
@@ -17,14 +17,14 @@ rows correspond to rows in mle and mu}
 
 \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")
+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)))
+(theta <- twSigmaLogitnorm(mle = seq(0.401,0.8,by = 0.1)))
 }