[R] convolution of the double exponential distribution
Bickel, David
DAVID.BICKEL at pioneer.com
Wed Dec 28 00:12:43 CET 2005
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): 159183.
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 HiBred International (DuPont)
Bioinformatics and Exploratory Research
7200 NW 62nd Ave.; PO Box 184
Johnston, IA 501310184
5153344739 Tel
5153344473 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; rhelp 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
 >
 >______________________________________________
 >Rhelp at stat.math.ethz.ch mailing list
 >https://stat.ethz.ch/mailman/listinfo/rhelp
 >PLEASE do read the posting guide!
 http://www.Rproject.org/postingguide.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(xloc)/scale)/(2*scale) }
 body(dfun) < substitute({ exp(abs(xloc)/scale)/(2*scale)
 }, list(loc
 = loc, scale = sca))

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

 qfun < function(x){ loc  scale*sign(x0.5)*log(1  2*abs(x0.5)) }
 body(qfun) < substitute({ loc  scale*sign(x0.5)*log(1 
 2*abs(x0.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}}
More information about the Rhelp
mailing list