Skip to content

convolution of the double exponential distribution

5 messages · Bickel, David, Duncan Murdoch, Matthias Kohl +2 more

#
Is there any R function that computes the convolution of the double
exponential distribution?

If not, is there a good way to integrate ((q+x)^n)*exp(-2x) over x from
0 to Inf for any value of q and for any positive integer n? I need to
perform the integration within a function with q and n as arguments. The
function integrate() is giving me this message:

"evaluation of function gave a result of wrong length"

David
_______________________________________
David R. Bickel  http://davidbickel.com
Research Scientist
Pioneer Hi-Bred International (DuPont)
Bioinformatics and Exploratory Research
7200 NW 62nd Ave.; PO Box 184
Johnston, IA 50131-0184
515-334-4739 Tel
515-334-4473 Fax
david.bickel at pioneer.com, bickel at prueba.info

This communication is for use by the intended recipient and ...{{dropped}}
#
On 12/22/2005 7:56 PM, Bickel, David wrote:
Under the substitution of y = q+x, that looks like a gamma integral. 
The x = 0 to Inf range translates into y = q to Inf, so you'll need an 
incomplete gamma function, such as pgamma.  Be careful to get the 
constant multiplier right.

Duncan Murdoch
#
Duncan Murdoch schrieb:
Hi,

you can use our package "distr".

require(distr)
## define double exponential distribution
loc <- 0 # location parameter
sca <- 1 # scale parameter

rfun <- function(n){ loc + scale * ifelse(runif(n) > 0.5, 1, -1) * rexp(n) }
body(rfun) <- substitute({ loc + scale * ifelse(runif(n) > 0.5, 1, -1) * 
rexp(n) },
                         list(loc = loc, scale = sca))

dfun <- function(x){ exp(-abs(x-loc)/scale)/(2*scale) }
body(dfun) <- substitute({ exp(-abs(x-loc)/scale)/(2*scale) }, list(loc 
= loc, scale = sca))

pfun <- function(x){ 0.5*(1 + sign(x-loc)*(1-exp(-abs(x-loc)/scale))) }
body(pfun) <- substitute({ 0.5*(1 + 
sign(x-loc)*(1-exp(-abs(x-loc)/scale))) },
                         list(loc = loc, scale = sca))
                        
qfun <- function(x){ loc - scale*sign(x-0.5)*log(1 - 2*abs(x-0.5)) }
body(qfun) <- substitute({ loc - scale*sign(x-0.5)*log(1 - 2*abs(x-0.5)) },
                         list(loc = loc, scale = sca))

D1 <- new("AbscontDistribution", r = rfun, d = dfun, p = pfun, q = qfun)
plot(D1)

D2 <- D1 + D1 # convolution based on FFT
plot(D2)

hth,
Matthias
#
Hi,

It is quite easy to integrate ((q+x)^n)*exp(-2x) over x from 0 to Inf, if
you are familiar with the following basic Laplace transform results:
(a) L_s[x^n] = \Gamma(n+1)/s^(n+1)
(b) L_s[f(x+q)] = exp(-qs) F(s)

Using (a) and (b) and substituting s=2, you get your desired integral, I:
I = exp(-2q) \Gamma(n+1) / 2^(n+1).

Hope this helps,
Ravi.

--------------------------------------------------------------------------
Ravi Varadhan, Ph.D.
Assistant Professor,  The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email:  rvaradhan at jhmi.edu
--------------------------------------------------------------------------
#
Mathematica says

Assuming[q âˆˆ Reals && q > 0, Integrate[(q + x)^n*Exp[-2*x], {x, 0,
Infinity}]]

2^(-1 - n)*Exp[2*q]*Gamma[1 + n, 2*q]

and 2-argument Gamma is the incomplete Gamma function
(integration starting at 2*q)
Duncan Murdoch wrote: