Skip to content

Commit

Permalink
spellcheck for CRAN
Browse files Browse the repository at this point in the history
  • Loading branch information
bgctw committed Oct 16, 2018
1 parent 12e1c41 commit 2113349
Show file tree
Hide file tree
Showing 18 changed files with 249 additions and 80 deletions.
9 changes: 7 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -3,8 +3,13 @@ Title: Functions for the Lognormal Distribution
Version: 0.1.3
Author: Thomas Wutzler
Maintainer: Thomas Wutzler <[email protected]>
Description: Provides Moments and other statistics of the lognormal distribution.
Estimation of the parameters of the sum of lognormals is implemented.
Description: The lognormal distribution
(Limpert et al (2001) <doi:10.1641/0006-3568(2001)051[0341:lndats]2.0.co;2>)
can characterize uncertainty that is bounded by zero.
This package provides estimation of distribution parameters, computation of
moments and other basic statistics, and an approximation of the distribution
of the sum of several correlated lognormally distributed variables
(Lo 2013 <doi:10.12988/ams.2013.39511>).
Imports: Matrix
Suggests: testthat, knitr, dplyr, ggplot2, mvtnorm, purrr, tidyr
VignetteBuilder: knitr
Expand Down
55 changes: 43 additions & 12 deletions R/autocorr.R
Original file line number Diff line number Diff line change
Expand Up @@ -10,22 +10,22 @@ computeEffectiveNumObs <- function(
## stripped before the computation proceeds.
){
##references<<
## Zieba & Ramza (2011)
## \code{Zieba & Ramza (2011)
## Standard Deviation of the Mean of Autocorrelated
## Observations Estimated with the Use of the Autocorrelation Function
## Estimated From the Data.
## Metrology and Measurement Systems,
## Walter de Gruyter GmbH, 18 10.2478/v10178-011-0052-x
## Walter de Gruyter GmbH, 18 10.2478/v10178-011-0052-x}
##
## Bayley & Hammersley (1946)
## \code{Bayley & Hammersley (1946)
## The "effective" number of independent observations in an autocorrelated
## time series.
## Supplement to the Journal of the Royal Statistical Society, JSTOR,8,184-197
## Supplement to the Journal of the Royal Statistical Society, JSTOR,8,184-197}
#
##details<< Handling of NA values: NAs at the beginning or end and are
## just trimmed before computation and pose no problem.
## However with NAs aside from edges, the return value isbiased low,
## because correlatations terms are subtracted for those positions.
## However with NAs aside from edges, the return value is biased low,
## because correlation terms are subtracted for those positions.
resTr <- .trimNA(res)
if (!isTRUE(na.rm) & any(is.na(resTr))) return(NA_integer_)
isFin <- is.finite(resTr)
Expand All @@ -44,6 +44,19 @@ computeEffectiveNumObs <- function(
##value<< integer scalar: effective number of observations
nEff
}
attr(computeEffectiveNumObs,"ex") <- function(){
# generate autocorrelated time series
res <- stats::filter(rnorm(1000), filter = rep(1,5), circular = TRUE)
res[100:120] <- NA
# plot the series of autocorrelated random variables
plot(res)
# plot their empirical autocorrelation function
acf(res, na.action = na.pass)
#effAcf <- computeEffectiveAutoCorr(res)
# the effective number of parameters is less than number of 1000 samples
(nEff <- computeEffectiveNumObs(res, na.rm = TRUE))
}


.trimNA <- function(
### remove NA values at the start and end
Expand Down Expand Up @@ -72,28 +85,46 @@ computeEffectiveAutoCorr <- function(
##details<<
## Returns all components before first negative autocorrelation
##references<<
## Zieba 2011 Standard Deviation of the Mean of Autocorrelated
## \code{Zieba 2011 Standard Deviation of the Mean of Autocorrelated
## Observations Estimated with the Use of the Autocorrelation Function Estimated
## From the Data
## From the Data}
n <- length(res)
if (!is.finite(nC)) return(c(1))
##value<< numeric vector: stongest compponents of the autocorrelation function
##value<< numeric vector: strongest components of the autocorrelation function
ans$acf[1:nC]
}
attr(computeEffectiveAutoCorr,"ex") <- function(){
# generate autocorrelated time series
res <- stats::filter(rnorm(1000), filter = rep(1,5), circular = TRUE)
res[100:120] <- NA
(effAcf <- computeEffectiveAutoCorr(res))
}


#' @export
varEffective <- function(
### estimate the variance of a correlated time series
res ##<< numeric of autocorrelated numbers, usually observation - model residuals
, nEff = computeEffectiveNumObs(res) ##<< effective number of observations
, nEff = computeEffectiveNumObs(res, na.rm = na.rm) ##<< effective
## number of observations
, na.rm = FALSE ##<< set to TRUE to remove NA cases before computation
, ... ##<< further arguments to \code{\link{var}}
) {
##details<< The BLUE is not anymore the usual variance, but a modified
## variance as given in Zieba 2011
## variance as given in \code{Zieba 2011}
a <- 1
var(res, ...) * nEff/(nEff - 1)
var(res, na.rm = na.rm, ...) * nEff/(nEff - 1)
### The estimated variance of the sample
}
attr(varEffective,"ex") <- function(){
# generate autocorrelated time series
res <- stats::filter(rnorm(1000), filter = rep(1,5), circular = TRUE)
res[100:120] <- NA
# if correlations are neglected, the estimate of the variance is biased low
(varNeglectCorr <- var(res, na.rm = TRUE))
(varCorr <- varEffective(res, na.rm = TRUE))
}


.tmp.f <- function(){
dss <- dsfP %>% mutate(
Expand Down
53 changes: 45 additions & 8 deletions R/lognorm.R
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
#' @export
getLognormMoments <- function(
### get the expected value and variance of a log-normal distribution
mu ##<< center parameter (mean at log scale, log(median))
, sigma ##<< scale parameter (standard deviation at log scale)
mu ##<< numeric vector of center parameter (mean at log scale, log(median))
, sigma ##<< numeric vector of scale parameter (standard deviation at log scale)
){
##references<< Limpert E, Stahel W & Abbt M (2001)
##references<< \code{Limpert E, Stahel W & Abbt M (2001)
## Log-normal Distributions across the Sciences: Keys and Clues.
## Oxford University Press (OUP) 51, 341,
## 10.1641/0006-3568(2001)051[0341:lndats]2.0.co;2
## 10.1641/0006-3568(2001)051[0341:lndats]2.0.co;2}
sigma2 <- sigma*sigma
m = exp(mu + sigma2/2)
v = (exp(sigma2) - 1)*exp(2*mu + sigma2)
Expand All @@ -18,6 +18,16 @@ getLognormMoments <- function(
, cv = as.vector(sqrt(v)/m) ##<< coefficient of variation: std/mean
)
}
attr(getLognormMoments,"ex") <- function(){
# start by estimating lognormal parameters from moments
.mean <- 1
.var <- c(1.3,2)^2
parms <- getParmsLognormForMoments(.mean, .var)
#
# computed moments must equal previous ones
(ans <- getLognormMoments(parms[,"mu"], parms[,"sigma"]))
cbind(.var, ans[,"var"])
}

#' @export
getLognormMedian <- function(
Expand All @@ -28,6 +38,9 @@ getLognormMedian <- function(
##value<< the median
exp(mu)
}
attr(getLognormMedian,"ex") <- function(){
getLognormMedian(mu = log(1), sigma = log(2))
}

#' @export
getLognormMode <- function(
Expand All @@ -38,19 +51,24 @@ getLognormMode <- function(
##value<< the mode
exp(mu - sigma*sigma)
}
attr(getLognormMode,"ex") <- function(){
# with larger sigma, the distribution is more skewed
# with mode further away from median = 1
getLognormMode(mu = log(1), sigma = c(log(1.2),log(2)))
}

#' @export
getParmsLognormForMoments <- function(
### get the mean and variance of a log-normal distribution
mean ##<< expected value at original scale
, var ##<< variance at original scale
, sigmaOrig = sqrt(var) ##<< alternatively to the variance,
## the standard devation at original scale can be given
, sigmaOrig = sqrt(var) ##<< standard deviation at original scale
##, can be specified alternatively to the variance
){
##references<< Limpert E, Stahel W & Abbt M (2001)
##references<< \code{Limpert E, Stahel W & Abbt M (2001)
## Log-normal Distributions across the Sciences: Keys and Clues.
## Oxford University Press (OUP) 51, 341,
## 10.1641/0006-3568(2001)051[0341:lndats]2.0.co;2
## 10.1641/0006-3568(2001)051[0341:lndats]2.0.co;2}
omega = 1 + (sigmaOrig/mean)^2
mu = log(mean / sqrt(omega))
sigma = sqrt(log(omega))
Expand All @@ -62,6 +80,12 @@ getParmsLognormForMoments <- function(
## (standard deviation at log scale)
)
}
attr(getParmsLognormForMoments,"ex") <- function(){
.mean <- 1
.var <- c(1.3,2)^2
getParmsLognormForMoments(.mean, .var)
}


#' @export
getParmsLognormForExpval <- function(
Expand All @@ -78,6 +102,13 @@ getParmsLognormForExpval <- function(
## (standard deviation at log scale)
)
}
attr(getParmsLognormForExpval,"ex") <- function(){
.mean <- 1
.sigmaStar <- c(1.3,2)
(parms <- getParmsLognormForExpval(.mean, .sigmaStar))
# multiplicative standard deviation must equal the specified value
cbind(exp(parms[,"sigma"]), .sigmaStar)
}

#' @export
estimateParmsLognormFromSample <- function(
Expand All @@ -94,3 +125,9 @@ estimateParmsLognormFromSample <- function(
## (standard deviation at log scale)
)
}
attr(estimateParmsLognormFromSample,"ex") <- function(){
.mu <- log(1)
.sigma <- log(2)
x <- exp(rnorm(50, mean = .mu, sd = .sigma))
estimateParmsLognormFromSample(x)
}
31 changes: 23 additions & 8 deletions R/lognormalSum.R
Original file line number Diff line number Diff line change
Expand Up @@ -7,10 +7,10 @@ estimateSumLognormalSample <- function(
## to estimate correlation
, effAcf = computeEffectiveAutoCorr(resLog) ##<< effective autocorrelation
## coefficients (may provide precomputed for efficiency or if the sample
## of resLog is too small) set to 1 to assume uncorrelated sample
## of \code{resLog} is too small) set to 1 to assume uncorrelated sample
, isGapFilled = logical(0) ##<< logical vector whether entry is gap-filled
## rather than an original measurement, see details
, na.rm = TRUE ##<< neglect terms with NA vlaues in mu or sigma
, na.rm = TRUE ##<< neglect terms with NA values in mu or sigma
){
# only one term, return the parameters
if (length(mu) == 1 ) return(
Expand Down Expand Up @@ -49,7 +49,8 @@ estimateSumLognormalSample <- function(
mu, sigma = sigma, sigmaSum = sigmaSum, effAcf = effAcf, na.rm = na.rm)
return(c(p , nEff = nEff))
}
##value<< numeric vector with components "mu", "sigma", and "nEff"
##value<< numeric vector with components \code{mu}, \code{sigma},
## and \code{nEff},
## the parameters of the lognormal distribution at log scale
## (Result of \code{link{estimateSumLognormal}})
## and the number of effective observations.
Expand All @@ -72,20 +73,20 @@ estimateSumLognormal <- function(
## to speed up computation by neglecting correlations among terms
## further apart.
## Set to zero to omit correlations.
, isStopOnNoTerm = FALSE ##<< if no finite estiamte is provide, by
## default, NA is returned for the sum.
, isStopOnNoTerm = FALSE ##<< if no finite estimate is provided then by
## default NA is returned for the sum.
## Set this to TRUE to issue an error instead.
, effAcf ##<< numeric vector of effective autocorrelation
## This overides arguments \code{corr} and \code{corrLength}
## This overrides arguments \code{corr} and \code{corrLength}
, na.rm = isStopOnNoTerm ##<< if there are terms with NA values in mu or sigma
## by default also the sum coefficients are NA. Set to TRUE to
## neglect such terms in the sum.
){
##references<<
## Lo C (2013) WKB approximation for the sum of two correlated lognormal
## \code{Lo C (2013) WKB approximation for the sum of two correlated lognormal
## random variables.
## Applied Mathematical Sciences, Hikari, Ltd., 7 , 6355-6367
## 10.12988/ams.2013.39511
## 10.12988/ams.2013.39511}
lengthMu <- length(mu)
if (length(sigma) == 1) sigma <- rep(sigma, lengthMu)
iFinite <- which( is.finite(mu) & is.finite(sigma))
Expand Down Expand Up @@ -136,6 +137,20 @@ estimateSumLognormal <- function(
## the parameters of the lognormal distribution at log scale
return(c(mu = as.vector(muSum), sigma = as.vector(sqrt(sigma2Eff))))
}
attr(estimateSumLognormal,"ex") <- function(){
# distribution of the sum of two lognormally distributed random variables
mu1 = log(110)
mu2 = log(100)
sigma1 = log(1.2)
sigma2 = log(1.6)
(coefSum <- estimateSumLognormal( c(mu1,mu2), c(sigma1,sigma2) ))
# repeat with correlation
(coefSumCor <- estimateSumLognormal( c(mu1,mu2), c(sigma1,sigma2), effAcf = c(1,0.9) ))
# expected value is equal, but variance with correlated variables is larger
getLognormMoments(coefSum["mu"],coefSum["sigma"])
getLognormMoments(coefSumCor["mu"],coefSumCor["sigma"])
}


estimateSumLognormalBenchmark <- function(
### Estimate the distribution parameters of the lognormal approximation to the sum
Expand Down
19 changes: 16 additions & 3 deletions cran-comments.md
Original file line number Diff line number Diff line change
@@ -1,13 +1,26 @@
## Note
Improved consistency with other R distribution methods.
Now with including examples in addition to vignettes.

Provide to the R community

* Approximating the sum of lognormals
* estimating distribution parameters from observations statistics
* computing moments and other statistics

## Test environments
* local R 3.4.4 on Mint17 64bit
* Ubuntu 14.04.5 LTS (on travis-ci), 3.5.0
* Ubuntu 14.04.5 LTS (on travis-ci), 14.04
* win-builder (devel and release)

## R CMD check results
There were no ERRORs nor WARNINGs nor NOTEs
There were no ERRORs nor WARNINGs
One note that this is a new submssion.

The first two possible misspellings in DESCRIPTION are caused by names in a
literature reference.
Third: Throughout this package I consistently write lognormal instead of log-normal.





8 changes: 4 additions & 4 deletions inst/genData/lognorm-package.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -41,14 +41,14 @@ Have a look at the package vignettes.
}%details

\references{
Limpert E, Stahel W & Abbt M (2001) Log-normal Distributions across the Sciences:
\code{Limpert E, Stahel W & Abbt M (2001) Log-normal Distributions across the Sciences:
Keys and Clues.BioScience, Oxford University Press (OUP), 51 , 341
10.1641/0006-3568(2001)051[0341:lndats]2.0.co;2
10.1641/0006-3568(2001)051[0341:lndats]2.0.co;2}

Lo C (2013) WKB approximation for the sum of two correlated lognormal
\code{Lo C (2013) WKB approximation for the sum of two correlated lognormal
random variables.
Applied Mathematical Sciences, Hikari, Ltd., 7 , 6355-6367
10.12988/ams.2013.39511
10.12988/ams.2013.39511}
}

\author{Thomas Wutzler}
Expand Down
13 changes: 9 additions & 4 deletions man/computeEffectiveAutoCorr.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -7,13 +7,18 @@
\item{res}{numeric of autocorrelated numbers, usually observation - model residuals}
}
\details{Returns all components before first negative autocorrelation}
\value{numeric vector: stongest compponents of the autocorrelation function}
\references{Zieba 2011 Standard Deviation of the Mean of Autocorrelated
\value{numeric vector: strongest components of the autocorrelation function}
\references{\code{Zieba 2011 Standard Deviation of the Mean of Autocorrelated
Observations Estimated with the Use of the Autocorrelation Function Estimated
From the Data}
From the Data}}
\author{Thomas Wutzler}





\examples{
# generate autocorrelated time series
res <- stats::filter(rnorm(1000), filter = rep(1,5), circular = TRUE)
res[100:120] <- NA
(effAcf <- computeEffectiveAutoCorr(res))
}
Loading

0 comments on commit 2113349

Please sign in to comment.