lognorm: Functions for the Lognormal Distribution

Description

The lognormal distribution (Limpert et al. (2001) <doi:10.1641/0006-3568(2001)051 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>) and the approximation of the difference of two correlated lognormally distributed variables (Lo 2012 <doi:10.1155/2012/838397>).

Details

Essential functions are

• Compute summary statistics: getLognormMoments

• computing distribution parameters from summary statistics: getParmsLognormForModeAndUpper

• estimate distribution parameters from summary sample: estimateParmsLognormFromSample

• Approximate the sum of correlated lognormals: estimateSumLognormalSample

• Approximate the difference of two correlated lognormals: estimateDiffLognormal

Utilities for correlated data. These functions maybe moved to a separate package in future.

• Estimate summary statistics of autocorrelated data

  • standard error of the mean: seCor

  • effective number of observations computeEffectiveNumObs

  • variance: varCor

• Return the vector of effective components of the autocorrelation: computeEffectiveAutoCorr

Otherwise refer to the vignettes

See Also

Useful links:

• https://github.com/bgctw/lognorm

remove NA values at the start and end

Description

remove NA values at the start and end

Usage

.trimNA(x)

Arguments

Value

subset of x with leading and trailing NAs removed

Estimate vector of effective components of the autocorrelation

Description

Estimate vector of effective components of the autocorrelation

Usage

computeEffectiveAutoCorr(res, type = "correlation")

Arguments

Details

Returns all components before first negative autocorrelation

Value

numeric vector: strongest components of the autocorrelation function

References

Zieba 2011 Standard Deviation of the Mean of Autocorrelated Observations Estimated with the Use of the Autocorrelation Function Estimated From the Data

Examples

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

Compute the effective number of observations taking into account autocorrelation

Description

Compute the effective number of observations taking into account autocorrelation

Usage

computeEffectiveNumObs(
  res,
  effAcf = computeEffectiveAutoCorr(res),
  na.rm = FALSE,
  exact.na = TRUE
)

Arguments

Details

Assumes records of all times present. DO NOT REMOVE OR FILTER NA records before. The length of the time series is used.

Handling of NA values: The formula from Zieba 2011 is extended to subtract the number of missing pairs in the count of correlation terms. If 'exact.na=false' the original formula is used (after trimming edge-NAs).

Value

integer scalar: effective number of observations

References

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

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

Examples

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

Inference on the difference of two lognormals

Description

The distribution of y = a - b + s, where a and b are two lognormal random variables and s is a constant to be estimated, can be approximated by a lognormal distribution.

Usage

estimateDiffLognormal(mu_a, mu_b, sigma_a, sigma_b, corr = 0)

pDiffLognormalSample(
  mu_a,
  mu_b,
  sigma_a,
  sigma_b,
  corr = 0,
  q = 0,
  nSample = 1e+05
)

Arguments

Value

estimateDiffLognormal: numeric vector with components mu, sigma, and shift, the components of the shifted lognormal distribution.

pDiffLognormalSample: vector of probabilities

Functions

• estimateDiffLognormal: Estimate the shifted-lognormal approximation to difference of two lognormals

• pDiffLognormalSample: Distribution function for the difference of two lognormals based on sampling. Default provides the probability that the difference is significantly larger than zero.

Estimate lognormal distribution parameters from a sample

Description

Estimate lognormal distribution parameters from a sample

Usage

estimateParmsLognormFromSample(x, na.rm = FALSE)

estimateStdErrParms(x, na.rm = FALSE)

Arguments

Details

The expected value of a can be determined with higher accuracy the larger the sample. Here, the uncorrelated assumption is applied at the log scale and distribution parameters are returned with the same expected value as the sample, but with uncertainty (sigma) decreased by sqrt(nfin - 1).

Since with low relative error, the lognormal becomes very close to the normal distribution, the distribution of the mean can be well approximated by a normal with sd(mean(x)) ~ sd(x)/sqrt(n-1).

Value

numeric vector with components mu and sigma, i.e., the center parameter (mean at log scale, log(median)) and the scale parameter (standard deviation at log scale)

Functions

• estimateParmsLognormFromSample: Estimate lognormal distribution parameters from a sample

• estimateStdErrParms: Estimate parameters of the lognormal distribution of the mean from an uncorrelated sample

Examples

.mu <- log(1)
.sigma <- log(2)
n = 200
x <- exp(rnorm(n, mean = .mu, sd = .sigma))
exp(pL <- estimateParmsLognormFromSample(x)) # median and multiplicative stddev
c(mean(x), meanx <- getLognormMoments(pL["mu"],pL["sigma"])[,"mean"])
c(sd(x), sdx <- sqrt(getLognormMoments(pL["mu"],pL["sigma"])[,"var"]))

# stddev decreases (each sample about 0.9) to about 0.07
# for the mean with n replicated samples
se <- estimateStdErrParms(x)
sqrt(getLognormMoments(se["mu"],se["sigma"])[,"var"])
sd(x)/sqrt(n-1) # well approximated by normal
# expected value stays the same
c(meanx, getLognormMoments(se["mu"],se["sigma"])[,"mean"])

Estimate the parameters of the lognormal approximation to the sum

Description

Estimate the parameters of the lognormal approximation to the sum

Estimate the parameters of the lognormal approximation to the sum

Usage

estimateSumLognormalSample(
  mu,
  sigma,
  resLog,
  effAcf = computeEffectiveAutoCorr(resLog),
  isGapFilled = logical(0),
  na.rm = TRUE
)

estimateSumLognormalSampleExpScale(mean, sigmaOrig, ...)

estimateSumLognormal(
  mu,
  sigma,
  effAcf = c(),
  corr = Diagonal(length(mu)),
  corrLength = if (inherits(corr, "ddiMatrix")) 0 else nTerm,
  sigmaSum = numeric(0),
  isStopOnNoTerm = FALSE,
  na.rm = isStopOnNoTerm
)

Arguments

Details

If there are no gap-filled values, i.e. all(!isGapFilled) or !length(isGapFilled) (the default), distribution parameters are estimated using all the samples. Otherwise, the scale parameter (uncertainty) is first estimated using only the non-gapfilled records.

Also use isGapFilled == TRUE for records, where sigma cannot be trusted. When setting sigma to missing, this is also affecting the expected value.

If there are only gap-filled records, assume uncertainty to be (before v0.1.5: the largest uncertainty of given gap-filled records.) the mean of the given multiplicative standard deviation

Value

numeric vector with components mu, sigma, and nEff, i.e. the parameters of the lognormal distribution at log scale and the number of effective observations.

Functions

• estimateSumLognormalSample: In addition to estimateSumLognormal take care of missing values and estimate correlation terms.

• estimateSumLognormalSampleExpScale: Before calling estimateSumLognormalSample estimate lognormal parameters from value and its uncertainty given on original scale.

• estimateSumLognormal: Estimate the parameters of the lognormal approximation to the sum

References

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

Examples

  # 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"])

Construct the full correlation matrix from autocorrelation components.

Description

Construct the full correlation matrix from autocorrelation components.

Usage

getCorrMatFromAcf(nRow, effAcf)

Arguments

Compute summary statistics of a log-normal distribution

Description

Compute summary statistics of a log-normal distribution

Usage

getLognormMoments(mu, sigma, m = exp(mu + sigma2/2) - shift, shift = 0)

getLognormMedian(mu, sigma, shift = 0)

getLognormMode(mu, sigma, shift = 0)

Arguments

Value

for getLognormMoments a numeric matrix with columns mean (expected value at original scale) , var (variance at original scale) , and cv (coefficient of variation: sqrt(var)/mean). For the other functions a numeric vector of the required summary.

Functions

• getLognormMoments: get the expected value, variance, and coefficient of variation

• getLognormMedian: get the median

• getLognormMode: get the mode

References

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

See Also

scaleLogToOrig

Examples

  # 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"])
  #
  getLognormMedian(mu = log(1), sigma = log(2))
  getLognormMode(mu = log(1), sigma = c(log(1.2),log(2)))

Calculate mu and sigma of lognormal from summary statistics.

Description

Calculate mu and sigma of lognormal from summary statistics.

Usage

getParmsLognormForMedianAndUpper(median, upper, sigmaFac = qnorm(0.99))

getParmsLognormForMeanAndUpper(mean, upper, sigmaFac = qnorm(0.99))

getParmsLognormForLowerAndUpper(lower, upper, sigmaFac = qnorm(0.99))

getParmsLognormForLowerAndUpperLog(lowerLog, upperLog, sigmaFac = qnorm(0.99))

getParmsLognormForModeAndUpper(mle, upper, sigmaFac = qnorm(0.99))

getParmsLognormForMoments(mean, var, sigmaOrig = sqrt(var))

getParmsLognormForExpval(mean, sigmaStar)

Arguments

Details

For getParmsLognormForMeanAndUpper there are two valid solutions, and the one with lower sigma , i.e. the not so strongly skewed solution is returned.

Value

numeric matrix with columns 'mu' and 'sigma', the parameter of the lognormal distribution. Rows correspond to rows of inputs.

Functions

• getParmsLognormForMedianAndUpper: Calculates mu and sigma of lognormal from median and upper quantile.

• getParmsLognormForMeanAndUpper: Calculates mu and sigma of lognormal from mean and upper quantile.

• getParmsLognormForLowerAndUpper: Calculates mu and sigma of lognormal from lower and upper quantile.

• getParmsLognormForLowerAndUpperLog: Calculates mu and sigma of lognormal from lower and upper quantile at log scale.

• getParmsLognormForModeAndUpper: Calculates mu and sigma of lognormal from mode and upper quantile.

• getParmsLognormForMoments: Calculate mu and sigma from moments (mean anc variance)

• getParmsLognormForExpval: Calculate mu and sigma from expected value and geometric standard deviation

References

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

Examples

# example 1: a distribution with mode 1 and upper bound 5
(thetaEst <- getParmsLognormForModeAndUpper(1,5))
mle <- exp(thetaEst[1] - thetaEst[2]^2)
all.equal(mle, 1, check.attributes = FALSE)

# plot the distributions
xGrid = seq(0,8, length.out = 81)[-1]
dxEst <- dlnorm(xGrid, meanlog = thetaEst[1], sdlog = thetaEst[2])
plot( dxEst~xGrid, type = "l",xlab = "x",ylab = "density")
abline(v = c(1,5),col = "gray")

# example 2: true parameters, which should be rediscovered
theta0 <- c(mu = 1, sigma = 0.4)
mle <- exp(theta0[1] - theta0[2]^2)
perc <- 0.975		# some upper percentile, proxy for an upper bound
upper <- qlnorm(perc, meanlog = theta0[1], sdlog = theta0[2])
(thetaEst <- getParmsLognormForModeAndUpper(
  mle,upper = upper,sigmaFac = qnorm(perc)) )

#plot the true and the rediscovered distributions
xGrid = seq(0,10, length.out = 81)[-1]
dx <- dlnorm(xGrid, meanlog = theta0[1], sdlog = theta0[2])
dxEst <- dlnorm(xGrid, meanlog = thetaEst[1], sdlog = thetaEst[2])
plot( dx~xGrid, type = "l")
#plot( dx~xGrid, type = "n")
#overplots the original, coincide
lines( dxEst ~ xGrid, col = "red", lty = "dashed")

# example 3: explore varying the uncertainty (the upper quantile)
x <- seq(0.01,1.2,by = 0.01)
mle = 0.2
dx <- sapply(mle*2:8,function(q99){
  theta = getParmsLognormForModeAndUpper(mle,q99,qnorm(0.99))
  #dx <- dDistr(x,theta[,"mu"],theta[,"sigma"],trans = "lognorm")
  dx <- dlnorm(x,theta[,"mu"],theta[,"sigma"])
})
  matplot(x,dx,type = "l")
# Calculate mu and sigma from expected value and geometric standard deviation
.mean <- 1
.sigmaStar <- c(1.3,2)
(parms <- getParmsLognormForExpval(.mean, .sigmaStar))
# multiplicative standard deviation must equal the specified value
cbind(exp(parms[,"sigma"]), .sigmaStar)

Scale standard deviation between log and original scale.

Description

When comparing values at log scale that have different sd at original scale, better compare log(mean) instead of mu.

Usage

scaleLogToOrig(logmean, sigma)

scaleOrigToLog(mean, sd)

Arguments

Value

numeric matrix with columns mean, and sd at original scale

Functions

• scaleLogToOrig: get logmean and sigma at log scale

• scaleOrigToLog: get mean and sd at original scale

Examples

  xLog <- data.frame(logmean = c(0.8, 0.8), sigma = c(0.2, 0.3))
  xOrig <- as.data.frame(scaleLogToOrig(xLog$logmean, xLog$sigma))
  xLog2 <- as.data.frame(scaleOrigToLog(xOrig$mean, xOrig$sd))
  all.equal(xLog, xLog2)
  xLog3 <- as.data.frame(getParmsLognormForMoments(xOrig$mean, xOrig$sd^2))
  all.equal(xLog$sigma, xLog3$sigma) # but mu  < logmean 

Compute the standard error accounting for empirical autocorrelations

Description

Compute the standard error accounting for empirical autocorrelations

Usage

seCor(
  x,
  effCor = if (missing(effCov)) computeEffectiveAutoCorr(x) else effCov/var(x, na.rm =
    TRUE),
  na.rm = FALSE,
  effCov,
  nEff = computeEffectiveNumObs(x, effCor, na.rm = na.rm)
)

Arguments

Details

The default uses empirical autocorrelation estimates from the supplied data up to first negative component. For short series of x it is strongly recommended to to provide effCov that was estimated on a longer time series.

Value

numeric scalar of standard error of the mean of x

set off-diagonal values of a matrix

Description

set off-diagonal values of a matrix

Usage

setMatrixOffDiagonals(x, diag = 1:length(value), value, isSymmetric = FALSE)

Arguments

Value

matrix with modified diagonal elements

Compute the unbiased variance accounting for empirical autocorrelations

Description

Compute the unbiased variance accounting for empirical autocorrelations

Usage

varCor(
  x,
  effCor = computeEffectiveAutoCorr(x),
  na.rm = FALSE,
  nEff = computeEffectiveNumObs(x, effAcf = effCor)
)

Arguments

Details

The default uses empirical autocorrelation estimates from the supplied data up to first negative component. For short series of x it is strongly recommended to to provide effCov that was estimated on a longer time series.

Value

numeric scalar of unbiased variation of x