mixdiff: Difference of mixture distributions
In RBesT: R Bayesian Evidence Synthesis Tools

Description Usage Arguments Details Value Examples

Density, cumulative distribution function, quantile function and random number generation for the difference of two mixture distributions.

dmixdiff(mix1, mix2, x)

pmixdiff(mix1, mix2, q, lower.tail = TRUE)

qmixdiff(mix1, mix2, p, lower.tail = TRUE)

rmixdiff(mix1, mix2, n)

`mix1`	first mixture density
`mix2`	second mixture density
`x`	vector of values for which density values are computed
`q`	vector of quantiles for which cumulative probabilities are computed
`lower.tail`	logical; if `TRUE` (default), probabilities are P[X <= x], otherwise P[X > x].
`p`	vector of cumulative probabilities for which quantiles are computed
`n`	size of random sample

If x_1 ~ f_1(x) and x_2 ~ f_2(x), the density of the difference x = x_1 - x_2 is given by the convolution

f(x) = \int f_1(x) f_2(x - u) du = (f_1 * f_2)(x).

The cumulative distribution function equates to

F(x) = \int F_1(x+u) f_2(u) du.

Both integrals are performed over the full support of the densities and use the numerical integration function integrate.

Respective density, quantile, cumulative density or random numbers.

# 1. Difference between two beta distributions, i.e. Pr( mix1 - mix2 > 0)
mix1 <- mixbeta(c(1, 11, 4))
mix2 <- mixbeta(c(1, 8, 7))
pmixdiff(mix1, mix2, 0, FALSE)

# Interval probability, i.e. Pr( 0.3 > mix1 - mix2 > 0)
pmixdiff(mix1, mix2, 0.3) - pmixdiff(mix1, mix2, 0)

# 2. two distributions, one of them a mixture
m1 <- mixbeta( c(1,30,50))
m2 <- mixbeta( c(0.75,20,50),c(0.25,1,1))

# random sample of difference
set.seed(23434)
rM <- rmixdiff(m1, m2, 1E4)

# histogram of random numbers and exact density
hist(rM,prob=TRUE,new=TRUE,nclass=40)
curve(dmixdiff(m1,m2,x), add=TRUE, n=51)

# threshold probabilities for difference, at 0 and 0.2
pmixdiff(m1, m2, 0)
mean(rM<0)
pmixdiff(m1,m2,0.2)
mean(rM<0.2)

# median of difference
mdn <- qmixdiff(m1, m2, 0.5)
mean(rM<mdn)

# 95%-interval
qmixdiff(m1, m2, c(0.025,0.975))
quantile(rM, c(0.025,0.975))