Re: [R] convolution of the double exponential distribution

From: Ravi Varadhan <RVARADHAN_at_jhmi.edu>
Date: Sat 24 Dec 2005 - 02:54:45 EST


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@jhmi.edu

> -----Original Message-----
> From: r-help-bounces@stat.math.ethz.ch [mailto:r-help-
> bounces@stat.math.ethz.ch] On Behalf Of Matthias Kohl
> Sent: Friday, December 23, 2005 10:09 AM
> To: Bickel, David
> Cc: r-help@stat.math.ethz.ch; Duncan Murdoch
> Subject: Re: [R] convolution of the double exponential distribution
>
> Duncan Murdoch schrieb:
>
> >On 12/22/2005 7:56 PM, Bickel, David wrote:
> >
> >
> >>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"
> >>
> >>
> >
> >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
> >
> >______________________________________________
> >R-help@stat.math.ethz.ch mailing list
> >
https://stat.ethz.ch/mailman/listinfo/r-help
> >PLEASE do read the posting guide! http://www.R-project.org/posting-
> guide.html
> >
> >
>
> 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
>
> --
> StaMatS - Statistik + Mathematik Service
> Dipl.Math.(Univ.) Matthias Kohl
> www.stamats.de
>
> ______________________________________________
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-
> guide.html


R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Sat Dec 24 03:04:16 2005

This archive was generated by hypermail 2.1.8 : Sat 24 Dec 2005 - 04:38:07 EST