Nothing
## ----eval=FALSE, include=FALSE------------------------------------------------
# # twDev::genVigs()
# #rmarkdown::render("lognormalDiff.Rmd","md_document")
## ----setup, include=FALSE-----------------------------------------------------
library(knitr)
opts_chunk$set(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_hooks$set(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))
## -----------------------------------------------------------------------------
# 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"])
## -----------------------------------------------------------------------------
getLognormMoments(coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])
getLognormMode(coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])
getLognormMedian(coefDiff["mu"], coefDiff["sigma"], shift = coefDiff["shift"])
## -----------------------------------------------------------------------------
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"])
)
## ----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())
## -----------------------------------------------------------------------------
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)
## ----plotSignifZero, echo=FALSE-----------------------------------------------
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())
## -----------------------------------------------------------------------------
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"])
## -----------------------------------------------------------------------------
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))
head(xObsN)
y = xObsN[,1] - xObsN[,2]
## ----pdfDiffCorr80, echo=FALSE, fig.height=2.04, fig.width=3.27---------------
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())
## -----------------------------------------------------------------------------
# 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))
## -----------------------------------------------------------------------------
# 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))
# note the minus sign in front
(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(
# note the minus sign in front
q = -(q_shifted_neg - coefDiff["shift"]),
d = dlnorm(q_shifted_neg, coefDiff["mu"], coefDiff["sigma"])
)
## ----plotDiffLargerVar, 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())
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.