# In lognorm: Functions for the Lognormal Distribution

# twDev::genVigs()
#rmarkdown::render("lognormalDiff.Rmd","md_document")

library(knitr)
opts_chunkset(out.extra = 'style="display:block; margin: auto"' #, fig.align = "center" #, fig.width = 4.6, fig.height = 3.2 , fig.width = 6, fig.height = 3.75 #goldener Schnitt 1.6 , dev.args = list(pointsize = 10) , dev = c('png','pdf') ) knit_hooksset(spar = function(before, options, envir) {
if (before) {
par( las = 1 )                   #also y axis labels horizontal
par(mar = c(2.0,3.3,0,0) + 0.3 )  #margins
par(tck = 0.02 )                          #axe-tick length inside plots
par(mgp = c(1.1,0.2,0) )  #positioning of axis title, axis labels, axis
}
})
library(lognorm)
if (!require(ggplot2) || !require(dplyr) || !require(tidyr) || !require(purrr)) {
print("To generate this vignette, ggplot2, dplyr, tidyr, and purrr are required.")
knit_exit()
}
themeTw <- ggplot2::theme_bw(base_size = 10) +
theme(axis.title = element_text(size = 9))


# Approximating the difference of lognormal random variables

Lo 2012 reports an approximation for the difference of two random variables by a shifted lognormal distribution.

Rather than approximating the density of $y = a - b$, it approximates the density of $y_s = a - b + s$, where $s$ is the shift. Hence, one has to subtract $s$ from provided mean and quantiles. One can can use the variance, and relative error but has to recompute the relative error.

## Two uncorrelated random variables

# generate nSample values of two lognormal random variables
mu1 = log(110)
mu2 = log(100)
sigma1 = 0.25
sigma2 = 0.15
#(coefDiff <- estimateSumLognormal( c(mu1,mu2), c(sigma1,sigma2) ))
(coefDiff <- estimateDiffLognormal(mu1,mu2, sigma1,sigma2, 0))
(expDiff <- getLognormMoments(coefDiff["mu"], coefDiff["sigma"])[,"mean"] -
coefDiff["shift"])


Several functions accept the shift argument to handle this already.

getLognormMoments(coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])
getLognormMode(coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])
getLognormMedian(coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])


For the functions from the stats package, the shifting has to be done manually.

p <- seq(0,1,length.out = 100)[-c(1,100)]
dsPredY <- data.frame(
var = "y",
q_shifted = qlnorm(p, coefDiff["mu"], coefDiff["sigma"] )
) %>%
mutate(
q = q_shifted - coefDiff["shift"],
d = dlnorm(q_shifted, coefDiff["mu"], coefDiff["sigma"])
)


A check by random numbers (dotted lines) shows close correspondence.

{R plotUncorr, echo=FALSE, fig.height=2.04, fig.width=3.27} nSample = 2000

ds <- data.frame( a = rlnorm(nSample, mu1, sigma1) , b = rlnorm(nSample, mu2, sigma2) ) %>% mutate( y = a - b ) dsw <- gather(ds, "var", "value", a, b, y) p1 <- dsw %>% filter(var == "y") %>% ggplot(aes(value, color = var)) + geom_density(linetype = "dotted") + geom_vline(xintercept = mean(ds$y), linetype = "dotted") p1 + geom_line(data = dsPredY, aes(q, d, color = var)) + geom_vline(aes(xintercept = expDiff, color = var), data = data.frame(var = "y", q = expDiff)) + themeTw + theme(legend.position = c(0.98,0.98), legend.justification = c(1,1)) + theme(axis.title.x = element_blank()) ## Test if difference is significantly different from zero The probability of the zero quantile needs to be larger than a significance level. We can compute it based on the lognormal approximation. Since Lo12 is only accurate if the expected difference is small compared to the expected sum, the probability of the difference being larger than zero can be estimated by a sampling both terms. r mu1 = log(120) mu2 = log(60) sigma1 = 0.25 sigma2 = 0.15 coefDiff <- estimateDiffLognormal( mu1,mu2,sigma1,sigma2, corr = -0.8 ) pLo <- plnorm(0 + coefDiff["shift"], coefDiff["mu"], coefDiff["sigma"]) pSample <- pDiffLognormalSample(mu1,mu2,sigma1,sigma2, corr = -0.8) c(pLo = as.numeric(pLo), pSample = pSample)  In the example both approaches give a probability of less than 5% so that we conclude that the difference is significant. p <- seq(0,1,length.out = 100)[-c(1,100)] dsPredY <- tibble( var = "y", q_shifted = qlnorm(p, coefDiff["mu"], coefDiff["sigma"] ), q = q_shifted - coefDiff["shift"], d = dlnorm(q_shifted, coefDiff["mu"], coefDiff["sigma"]) ) dsPredA <- tibble( var = "a", q = qlnorm(p, mu1, sigma1), d = dlnorm(q, mu1, sigma1) ) dsPredB <- tibble( var = "b", q = qlnorm(p, mu2, sigma2), d = dlnorm(q, mu2, sigma2) ) bind_rows(select(dsPredY, -q_shifted), dsPredA, dsPredB) %>% #filter(var == "a") %>% ggplot(aes(q, d, color = var)) + geom_line() + geom_vline(xintercept = 0, linetype = "dashed", color = "darkgrey") + themeTw + theme(legend.position = c(0.98,0.98), legend.justification = c(1,1)) + theme(axis.title.x = element_blank())  ## Two positively correlated variables if (!requireNamespace("mvtnorm")) { warning("Remainder of the vignette required mvtnorm installed.") knitr::opts_chunk$set(error = TRUE)
}
corr = 0.8
(coefDiff <- estimateDiffLognormal(mu1,mu2, sigma1,sigma2, corr = corr))
(expDiff <- getLognormMoments(
coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])[,"mean"])


Check with sampled distribution.

nSample <- 1e5
sigma_vec = c(sigma1, sigma2)
corrM <- setMatrixOffDiagonals(
diag(nrow = 2), value = corr, isSymmetric = TRUE)
covM <- diag(sigma_vec) %*% corrM %*% diag(sigma_vec)
xObsN <- exp(mvtnorm::rmvnorm(nSample, mean =  c(mu1, mu2), sigma = covM))
y = xObsN[,1] - xObsN[,2]

p <- seq(0,1,length.out = 100)[-c(1,100)]
dsPredY <- data.frame(
var = "y",
q_shifted = qlnorm(p, coefDiff["mu"], coefDiff["sigma"] )
) %>%
mutate(
q = q_shifted - coefDiff["shift"],
d = dlnorm(q_shifted, coefDiff["mu"], coefDiff["sigma"])
)
# density plot of the random draws
ggplot(data.frame(y = y), aes(y, color = "random draws")) +
geom_density() +
# line plot of the lognorm density approximation
geom_line(data = dsPredY, aes(q, d, color = "computed diff")) +
# expected value
geom_vline(xintercept = expDiff, linetype = "dashed") +
themeTw +
theme(legend.position = c(0.98,0.98), legend.justification = c(1,1)) +
theme(axis.title.x = element_blank()) +
theme(legend.title = element_blank())


The approximation for the difference of positively correlated random numbers predicts a narrower distribution than with the uncorrelated or negatively correlated difference above. However, this case less accurate and shows some deviations from the sampled distribution around the mode.

## Subtracting a variable with larger variance

The method of Lo12 requires $\sigma_b < \sigma_a$ and otherwise gives an error.

# generate nSample values of two lognormal random variables
mu1 = log(110)
mu2 = log(100)
sigma1 = 0.15
sigma2 = 0.25
try(coefDiff <- estimateDiffLognormal(mu1,mu2, sigma1,sigma2, 0))


But one can compute the density of $y_r = -y = b - a$ and plot the density of the shifted and negated distribution.

# note the switch of positions of mu1 and mu2: mu2 - mu1
#(coefDiff <- estimateDiffLognormal(mu1,mu2, sigma1,sigma2, 0)
(coefDiff <- estimateDiffLognormal(mu2,mu1, sigma2,sigma1, 0))
(expDiff <- -(getLognormMoments(
coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])[,"mean"]))

p <- seq(0,1,length.out = 100)[-c(1,100)]
dsPredY <- data.frame(
var = "y",
q_shifted_neg = qlnorm(p, coefDiff["mu"], coefDiff["sigma"] )
) %>%
mutate(
q = -(q_shifted_neg - coefDiff["shift"]),
d = dlnorm(q_shifted_neg, coefDiff["mu"], coefDiff["sigma"])
)

nSample = 2000
ds <- data.frame(
a = rlnorm(nSample, mu1, sigma1)
, b = rlnorm(nSample, mu2, sigma2)
) %>%  mutate(
y = a - b
)
dsw <- gather(ds, "var", "value", a, b, y)
p1 <- dsw %>% filter(var == "y") %>%
ggplot(aes(value, color = var)) + geom_density(linetype = "dotted") +
geom_vline(xintercept = mean(ds\$y), linetype = "dotted")
#

p1 + geom_line(data = dsPredY, aes(q, d, color = var)) +
geom_vline(aes(xintercept = expDiff, color = var),
data = data.frame(var = "y", q = expDiff)) +
themeTw +
theme(legend.position = c(0.98,0.98), legend.justification = c(1,1)) +
theme(axis.title.x = element_blank())


Because we subtract a large-variance lognormal variable, the distribution becomes right-skewed.

## Try the lognorm package in your browser

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

lognorm documentation built on March 11, 2021, 1:08 a.m.