[R] convolution of the double exponential distribution

[Bickel, David]

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 at pioneer.com, bickel at prueba.info


| -----Original Message-----
| From: Matthias Kohl [mailto:Matthias.Kohl at stamats.de] 
| Sent: Friday, December 23, 2005 9:09 AM
| To: Bickel, David
| Cc: Duncan Murdoch; r-help at 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 at 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}}