convolve: Convolution of two probability distributions

In bayesmeta: Bayesian Random-Effects Meta-Analysis and Meta-Regression

Convolution of two probability distributions

Description

Compute the convolution of two probability distributions, specified through their densities or cumulative distribution functions (CDFs).

Usage

  convolve(dens1, dens2,
           cdf1=Vectorize(function(x){integrate(dens1,-Inf,x)$value}),
           cdf2=Vectorize(function(x){integrate(dens2,-Inf,x)$value}),
           delta=0.01, epsilon=0.0001)

Arguments

dens1, dens2	the two distributions' probability density functions.
cdf1, cdf2	the two distributions' cumulative distribution functions.
delta, epsilon	the parameters specifying the desired accuracy for approximation of the convolution, and with that determining the number of support points being used internally. Smaller values imply greater accuracy and greater computational burden (Roever and Friede, 2017).

Details

The distribution of the sum of two (independent) random variables technically results as a convolution of their probability distributions. In some cases, the calculation of convolutions may be done analytically; e.g., the sum of two normally distributed random variables again turns out as normally distributed (with mean and variance resulting as the sums of the original ones). In other cases, convolutions may need to be determined numerically. One way to achieve this is via the DIRECT algorithm; the present implementation is the one discussed by Roever and Friede (2017). Accuracy of the computations is determined by the delta (maximum divergence \delta) and epsilon (tail probability \epsilon) parameters.

Convolutions here are used within the funnel() function (to generate funnel plots), but are often useful more generally. The original probability distributions may be specified via their probability density functions or their cumulative distribution functions (CDFs). The convolve() function returns the convolution's density, CDF and quantile function (inverse CDF).

Value

A list with elements

delta	the \delta parameter.
epsilon	the \epsilon parameter.
binwidth	the bin width.
bins	the total number of bins.
support	a matrix containing the support points used internally for the convolution approximation.
density	the probability density function.
cdf	the cumulative distribution function (CDF).
quantile	the quantile function (inverse CDF).

Author(s)

Christian Roever [e-mail]

References

C. Roever, T. Friede. Discrete approximation of a mixture distribution via restricted divergence. Journal of Computational and Graphical Statistics, 26(1):217-222, 2017. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1080/10618600.2016.1276840")}.

See Also

bayesmeta, funnel.

Examples

## Not run: 
#  Skew-normal / logistic example:

dens1 <- function(x, shape=4)
# skew-normal distribution's density
# see also: http://azzalini.stat.unipd.it/SN/Intro
{
  return(2 * dnorm(x) * pnorm(shape * x))
}

dens2 <- function(x)
# logistic distribution's density
{
  return(dlogis(x, location=0, scale=1))
}

rskewnorm <- function(n, shape=4)
# skew-normal random number generation
# (according to http://azzalini.stat.unipd.it/SN/faq-r.html)
{
  delta <- shape / sqrt(shape^2+1)
  u1 <- rnorm(n); v <- rnorm(n)
  u2 <- delta * u1 + sqrt(1-delta^2) * v
  return(apply(cbind(u1,u2), 1, function(x){ifelse(x[1]>=0, x[2], -x[2])}))
}

# compute convolution:
conv <- convolve(dens1, dens2)

# illustrate convolution:
n <- 100000
x <- rskewnorm(n)
y <- rlogis(n)
z <- x + y

# determine empirical and theoretical quantiles:
p      <- c(0.001,0.01, 0.1, 0.5, 0.9, 0.99, 0.999)
equant <- quantile(z, prob=p)
tquant <- conv$quantile(p)

# show numbers:
print(cbind("p"=p, "empirical"=equant, "theoretical"=tquant))

# draw Q-Q plot:
rg <- range(c(equant, tquant))
plot(rg, rg, type="n", asp=1, main="Q-Q-plot",
     xlab="theoretical quantile", ylab="empirical quantile")
abline(0, 1, col="grey")
points(tquant, equant)

## End(Not run)

bayesmeta documentation built on July 9, 2023, 5:12 p.m.