Re: [R] convolution of the double exponential distribution

From: Bickel, David <DAVID.BICKEL_at_pioneer.com>
Date: Wed 28 Dec 2005 - 10:12:43 EST


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 Received on Wed Dec 28 10:18:04 2005

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:41:42 EST