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 >[EMAIL PROTECTED], [EMAIL PROTECTED] > > >| -----Original Message----- >| From: Matthias Kohl [mailto:[EMAIL PROTECTED] >| 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: [EMAIL PROTECTED] 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