Nothing
### normexp.R ---
##----------------------------------------------------------------------
## Author: Brice Ozenne
## Created: maj 6 2020 (14:06)
## Version:
## Last-Updated: jun 28 2023 (14:14)
## By: Brice Ozenne
## Update #: 75
##----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
##----------------------------------------------------------------------
##
### Code:
## * exponential distribution
## cumulative distribution fuction for Z = X + \rho Y
## where X follows a standard normal distribution
## and Y an exponential distribution with rate parameter \lambda
## denoting \Phi the cumulative distribution function of the standard normal distribution:
## F_Z(z) = \Prob[X + \rho Y < z]
## = \int f(x,y) \Ind[x + \rho y < z] dx dy
########################### rho>0
## = \int_{x \in [-\inf;z]} f(x) \int_{y \in [0;(z-x)/\rho] f(y) dx dy
## = \int_{x \in [-\inf;z]} f(x) (1-\exp(-(z-x)*\lambda/\rho)) dx
## = \Phi(z) - \exp(-z\lambda/rho)/\sqrt{2\pi} \int_{x \in [-\inf;z]} \exp(-x^2/2)\exp(x*\lambda/\rho) dx
## = \Phi(z) - \exp(-z\lambda/rho+\lambda^2/(2\rho^2))/\sqrt{2\pi} \int_{x \in [-\inf;z]} \exp(-(x-\lambda/rho)^2/2) dx
## = \Phi(z) - \exp(-z\lambda/rho+\lambda^2/(2\rho^2)) \Phi(z-\lambda/\rho)
##' @title Cumulative Distribution Function of a Gaussian Variable Plus an Exponential Variable
##' @noRd
##'
##' @examples
##' \dontrun{
##' n <- 1e6
##'
##' ## rho > 0
##' mean(rnorm(n) + 1.5 * rexp(n, rate = 2) <= 0.1)
##' pnormexp(0.1, rate = 2, rho = 1.5)
##' mean(rnorm(n) + 1.5 * rexp(n, rate = 2) <= 0.9)
##' pnormexp(0.9, rate = 2, rho = 1.5)
##'
##' ## rho < 0
##' mean(rnorm(n) - 1.5 * rexp(n, rate = 2) <= 0.1)
##' pnormexp(0.1, rate = 2, rho = -1.5)
##' mean(rnorm(n) - 1.5 * rexp(n, rate = 2) <= 0.9)
##' pnormexp(0.9, rate = 2, rho = -1.5)
##' }
pnormexp <- function(q, rate, rho){
if(abs(rho)<1e-12){
out <- stats::pnorm(q)
}else if(rho>0){
out <- stats::pnorm(q) - exp(-(rate/rho)*q+(rate/rho)^2/2)*stats::pnorm(q, mean = rate/rho)
}else{
## mean(rnorm(1e5) + rho * rexp(1e5, rate = 2) <= q) ## MONTE CARLO
## cubature::adaptIntegrate(f = function(x){dnorm(x[1])*dexp(x[2], rate = rate)*(x[1]+rho*x[2]<q)}, lowerLimit = c(-10,-10), upperLimit = c(10,10)) ## NUMERIC INTEGRATION
## integrate(f = function(x){ dnorm(x) * sapply(x, function(iX){integrate(f = function(y){dexp(y, rate = rate)}, lower = (q-iX)/rho, upper = 10)$value})}, lower = -10, upper = 10)
## integrate(f = function(x){ dnorm(x) * (1-pexp(q = (q-x)/rho, rate = rate))}, lower = -10, upper = 10)
## 1 - integrate(f = function(x){ dnorm(x) * pexp(q = (q-x)/rho, rate = rate)}, lower = -10, upper = 10)$value
## 1 - integrate(f = function(x){ dnorm(x) * pexp(q = (q-x)/rho, rate = rate)}, lower = q, upper = 10)$value
## 1 - integrate(f = function(x){ dnorm(x) * (1 - exp(-(q-x)*rate/rho))}, lower = q, upper = 10)$value
## pnorm(q) + integrate(f = function(x){ dnorm(x) * exp(-(q-x)*rate/rho)}, lower = q, upper = 10)$value
out <- stats::pnorm(q) + exp(-(rate/rho)*q+(rate/rho)^2/2)*(1-stats::pnorm(q, mean = rate/rho))
}
return(out)
}
##' @title Density of a Gaussian Variable Plus an Exponential Variable
##' @noRd
##'
##' @examples
##' \dontrun{
##' n <- 1e6
##'
##' c(qnormexp(0.5, rate = 2, rho = 1.5), quantile(rnorm(n) + 1.5 * rexp(n, rate = 2), 0.5))
##' c(qnormexp(0.95, rate = 1/10, rho = 1.5), quantile(rnorm(n) + 1.5 * rexp(n, rate = 1/10), 0.95))
##'
##' c(qnormexp(0.5, rate = 2, rho = -1.5), quantile(rnorm(n) - 1.5 * rexp(n, rate = 2), 0.5))
##' c(qnormexp(0.95, rate = 1/10, rho = -1.5), quantile(rnorm(n) - 1.5 * rexp(n, rate = 1/10), 0.95))
##' }
qnormexp <- function(p, rate, rho){
if(abs(rho)<1e-12){
out <- stats::qnorm(p)
}else if(rho>0){
out <- sapply(p, function(iP){
stats::uniroot(function(x){pnormexp(x, rate = rate, rho = rho) - iP},
lower = stats::qnorm(iP),
upper = (stats::qnorm(iP)+3) + (stats::qexp(iP, rate = rate) + 5/rate))$root
})
## iP <- tail(p,1)
## pnormexp(stats::qnorm(iP), rate = rate, rho = rho) - iP
## pnormexp((stats::qnorm(iP)+3) + (stats::qexp(iP, rate = rate) + 5/rate), rate = rate, rho = rho) - iP
## hist(rnorm(1e4) + rho * rexp(1e4, rate = rate))
}else if(rho<0){
out <- sapply(p, function(iP){
stats::uniroot(function(x){pnormexp(x, rate = rate, rho = rho) - iP},
lower = (stats::qnorm(iP)-3) - (stats::qexp(iP, rate = rate) + 5/rate),
upper = stats::qnorm(iP))$root
})
}
return(out)
}
## * Weibull distribution
## cumulative distribution fuction for Z = X + \rho Y
## where X follows a standard normal distribution
## and Y a weibull distribution with scale parameter \lambda and shape parameter k
## denoting \Phi the cumulative distribution function of the standard normal distribution:
## F_Z(z) = \Prob[X + \rho Y < z]
## = \int f(x,y) \Ind[x + \rho y < z] dx dy
## = \int_{x \in [-\inf;z]} f(x) \int_{y \in [0;(z-x)/\rho] f(y) dx dy
## = \int_{x \in [-\inf;z]} f(x) (1-\exp(-(z-x)^k/(\lambda\rho)^k)) dx
## = \Phi(z) - \int_{x \in [-\inf;z]} \exp(-x^2/2-(z-x)^k/(\lambda\rho)^k)\sqrt{2\pi} dx
##' @title Cumulative Distribution Function of a Gaussian Variable Plus an Weibull Variable
##' @noRd
##'
##' @examples
##' \dontrun{
##' n <- 1e6
##'
##' pnormweibull(0.1, scale = 1/2, shape = 1, rho = 1.5)
##' pnormweibull(0.8, scale = 1/2, shape = 1, rho = 1.5)
##' mean(rnorm(n) + 1.5 * rweibull(n, scale = 1/2, shape = 1) <= 0.1)
##' mean(rnorm(n) + 1.5 * rweibull(n, scale = 1/2, shape = 1) <= 0.8)
##'
##' pnormweibull(0.1, scale = 1/2, shape = 1, rho = -1.5)
##' pnormweibull(0.8, scale = 1/2, shape = 1, rho = -1.5)
##' mean(rnorm(n) - 1.5 * rweibull(n, scale = 1/2, shape = 1) <= 0.1)
##' mean(rnorm(n) - 1.5 * rweibull(n, scale = 1/2, shape = 1) <= 0.8)
##'
##' pnormweibull(0.1, scale = 1/2, shape = 2, rho = -1.5)
##' pnormweibull(0.8, scale = 1/2, shape = 2, rho = -1.5)
##' mean(rnorm(n) - 1.5 * rweibull(n, scale = 1/2, shape = 2) <= 0.1)
##' mean(rnorm(n) - 1.5 * rweibull(n, scale = 1/2, shape = 2) <= 0.8)
##' }
pnormweibull <- function(q, scale, shape, rho){
if(abs(rho)<1e-12){
out <- stats::pnorm(q)
}else{
if(rho>0){
if(shape==1){
out <- stats::pnorm(q) - exp(-(1/(scale*rho))*q+(1/(scale*rho))^2/2)*stats::pnorm(q, mean = 1/(scale*rho))
}else{
I <- stats::integrate(f = function(x){exp(-x^2/2)/sqrt(2*pi)*exp(-((q-x)/(rho*scale))^shape)}, lower = min(-4,q - 7^(1/shape)*rho*scale), upper = q)
out <- stats::pnorm(q) - I$value
}
}else if(rho<0){
if(shape==1){
out <- stats::pnorm(q) + exp(-(1/(scale*rho))*q+(1/(scale*rho))^2/2)*(1-stats::pnorm(q, mean = 1/(scale*rho)))
}else{
I <- stats::integrate(f = function(x){exp(-x^2/2)/sqrt(2*pi)*exp(-((q-x)/(rho*scale))^shape)}, lower = q, upper = max(4,q - 7^(1/shape)*rho*scale))
out <- stats::pnorm(q) + I$value
}
}
}
return(out)
}
##' @title Density of a Gaussian Variable Plus an Weibull Variable
##' @noRd
##'
##' @examples
##' \dontrun{
##' n <- 5e6
##'
##' c(qnormweibull(0.5, scale = 1/2, shape = 1, rho = 1.5),
##' quantile(rnorm(n) + 1.5 * rweibull(n, scale = 1/2, shape = 1), 0.5))
##' c(qnormweibull(0.95, scale = 10, shape = 1, rho = 1.5),
##' quantile(rnorm(n) + 1.5 * rweibull(n, scale = 10, shape = 1), 0.95))
##'
##' c(qnormweibull(0.5, scale = 1/2, shape = 2, rho = 1.5),
##' quantile(rnorm(n) + 1.5 * rweibull(n, scale = 1/2, shape = 2), 0.5))
##' c(qnormweibull(0.95, scale = 10, shape = 2, rho = 1.5),
##' quantile(rnorm(n) + 1.5 * rweibull(n, scale = 10, shape = 2), 0.95))
##'
##' c(qnormweibull(0.5, scale = 1/2, shape = 1, rho = -1.5),
##' quantile(rnorm(n) - 1.5 * rweibull(n, scale = 1/2, shape = 1), 0.5))
##' c(qnormweibull(0.95, scale = 10, shape = 1, rho = -1.5),
##' quantile(rnorm(n) - 1.5 * rweibull(n, scale = 10, shape = 1), 0.95))
##'
##' c(qnormweibull(0.5, scale = 1/2, shape = 2, rho = -1.5),
##' quantile(rnorm(n) - 1.5 * rweibull(n, scale = 1/2, shape = 2), 0.5))
##' c(qnormweibull(0.95, scale = 10, shape = 2, rho = -1.5),
##' quantile(rnorm(n) - 1.5 * rweibull(n, scale = 10, shape = 2), 0.95))
##'
##' }
qnormweibull <- function(p, scale, shape, rho){
if(abs(rho)<1e-12){
out <- stats::qnorm(p)
}else if(rho>0){
out <- sapply(p, function(iP){
stats::uniroot(function(x){pnormweibull(x, scale = scale, shape = shape, rho = rho) - iP},
lower = stats::qnorm(iP),
upper = (stats::qnorm(iP)+3) + (stats::qweibull(iP, scale = scale, shape = shape) + 5*scale))$root
})
}else if(rho<0){
out <- sapply(p, function(iP){
stats::uniroot(function(x){pnormweibull(x, scale = scale, shape = shape, rho = rho) - iP},
lower = (stats::qnorm(iP)-3) - (stats::qweibull(iP, scale = scale, shape = shape) + 5*scale),
upper = stats::qnorm(iP))$root
})
}
return(out)
}
##----------------------------------------------------------------------
### normexp.R ends here
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.