Re: [R] convolution of the double exponential distribution

From: Matthias Kohl <Matthias.Kohl_at_stamats.de>
Date: Wed 28 Dec 2005 - 20:01:35 EST

Hi David,

you can even increase the precision of these computations by changing the variables
"DefaultNrFFTGridPointsExponent" (default: 12 -> 2^12 grid points) and "TruncQuantile"
(default: 1e-5) in our package "distr"; see distroptions()

You can for instance use

distroptions("DefaultNrFFTGridPointsExponent", 14) and
distroptions("TruncQuantile", 1e-8)

We checked our algorithm in case of Binomial, Poisson, Normal and Exponential distributions and found a very high precision; confer Section 5 of

http://www.stamats.de/comp.pdf

Moreover, for n-fold convolutions you can also use the function "convpow" which can be found under
http://www.uni-bayreuth.de/departments/math/org/mathe7/DISTR/nFoldConvolution.R

hth,
Matthias

>Ravi, Duncan, and Matthias,
>
>Thank you very much for the helpful replies. For convolutions with 2 or
>3 copies, I found that the CDFs from the distr package closely match the
>analytic results from this paper:
>K. Singh, M. Xie, and W. E. Strawderman, "Combining Information From
>Independent Sources Through Confidence Distributions," Annals of
>Statistics 33, no. 1 (2005): 159-183.
>
>That gives me confidence that the package will also work well for higher
>copy numbers. At least for me, using the package is much more convenient
>than programming all the needed integrals into R.
>
>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@pioneer.com, bickel@prueba.info
>
>
>| -----Original Message-----
>| From: Matthias Kohl [mailto:Matthias.Kohl@stamats.de]
>| Sent: Friday, December 23, 2005 9:09 AM
>| To: Bickel, David
>| Cc: Duncan Murdoch; r-help@stat.math.ethz.ch
>| 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
>|
>
>This communication is for use by the intended recipient and ...{{dropped}}
>
>______________________________________________
>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
>
>

-- 
StaMatS - Statistik + Mathematik Service
Dipl.Math.(Univ.) Matthias Kohl
Gottlieb-Keim-Straße 60
95448 Bayreuth
Phone: +49 921 50736 457
E-Mail: matthias.kohl@stamats.de
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
Received on Wed Dec 28 20:13:59 2005

This archive was generated by hypermail 2.1.8 : Wed 28 Dec 2005 - 23:58:02 EST