[Distr-commits] r654 - branches/distr-2.3/pkg/distr/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon May 17 00:39:26 CEST 2010
Author: ruckdeschel
Date: 2010-05-17 00:39:26 +0200 (Mon, 17 May 2010)
New Revision: 654
Added:
branches/distr-2.3/pkg/distr/R/Convpow.R
Removed:
branches/distr-2.3/pkg/distr/R/Convpow0.R
Log:
a renaming schedule : enforcing ending .R in convpow.R II
Copied: branches/distr-2.3/pkg/distr/R/Convpow.R (from rev 652, branches/distr-2.3/pkg/distr/R/Convpow0.R)
===================================================================
--- branches/distr-2.3/pkg/distr/R/Convpow.R (rev 0)
+++ branches/distr-2.3/pkg/distr/R/Convpow.R 2010-05-16 22:39:26 UTC (rev 654)
@@ -0,0 +1,292 @@
+##########################################################
+## Function for n-fold convolution
+## -- absolute continuous distribution --
+##########################################################
+
+##implentation of Algorithm 3.4. of
+# Kohl, M., Ruckdeschel, P., Stabla, T. (2005):
+# General purpose convolution algorithm for distributions
+# in S4-Classes by means of FFT.
+# Technical report, Feb. 2005. Also available in
+# http://www.uni-bayreuth.de/departments/math/org/mathe7/
+# /RUCKDESCHEL/pubs/comp.pdf
+
+
+setMethod("convpow",
+ signature(D1 = "AbscontDistribution"),
+ function(D1, N){
+ if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if (N==1) return(D1)
+ ##STEP 1
+
+ lower <- getLow(D1); upper <- getUp(D1);
+
+ ##STEP 2
+
+ ## binary logarithm of the effective number of gridpoints
+ m <- max(getdistrOption("DefaultNrFFTGridPointsExponent") -
+ floor(log(N)/log(2)),5)
+ M <- 2^m
+ Nl <-2^ceiling(log(N)/log(2))
+
+ h <- (upper-lower)/M
+ dp1 <- .discretizeP(D1, lower, upper, h)
+
+ ##STEP 3
+
+ dpn0 <- c(dp1, numeric((Nl-1)*M))
+ ##STEP 4
+
+ ## computation of DFT
+ fftdpn <- fft(dpn0)
+
+ ##STEP 5
+
+ ## convolution theorem for DFTs
+ dpn <- c(0,(Re(fft(fftdpn^N, inverse = TRUE)) / (Nl*M))[1:(N*M-N+2)])
+
+ x <- seq(from = N*lower+N/2*h, to = N*upper-N/2*h, by = h)
+ x <- c(x[1]-h, x[1], x+h)
+
+ ## density (steps 5--7)
+
+ dfun <- .makeDNew(x, dpn, h)
+
+ ## cdf (steps 5--7)
+ pfun <- .makePNew(x, dpn, h, .notwithLArg(D1))
+
+ ## continuity correction by h/2
+
+ ## quantile function
+ yL <- if (q(D1)(0) == -Inf) -Inf else N*lower
+ yR <- if (q(D1)(1) == Inf) Inf else N*upper
+ px.l <- pfun(x + 0.5*h)
+ px.u <- pfun(x + 0.5*h, lower.tail = FALSE)
+ qfun <- .makeQNew(x + 0.5*h, px.l, px.u, .notwithLArg(D1), yL, yR)
+
+ rfun = function(n) colSums(matrix(r(D1)(n*N), ncol=n))
+
+ object <- new("AbscontDistribution", r = rfun, d = dfun, p = pfun,
+ q = qfun, .withArith = TRUE, .withSim = FALSE)
+
+ if(is(D1 at Symmetry,"SphericalSymmetry"))
+ object at Symmetry <- SphericalSymmetry(N*SymmCenter(D1 at Symmetry))
+
+ rm(m, dpn, dp1, dpn0, fftdpn)
+ rm(h, px.u, px.l, rfun, dfun, qfun, pfun, upper, lower)
+ return(object)
+})
+
+setMethod("convpow",
+ signature(D1 = "LatticeDistribution"),
+ function(D1, N, ep = getdistrOption("TruncQuantile")){
+ if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if (N==1) return(D1)
+
+ if(!is.numeric(ep)) stop("argument 'ep' must be a numeric.")
+ if(length(ep)!=1) stop("argument 'ep' must be a numeric of length 1.")
+ if((ep<0)||(ep>1)) stop("argument 'ep' must be in (0,1).")
+
+ w <- width(lattice(D1))
+
+ supp0 <- support(D1)
+ supp1 <- seq(by=abs(w),from=N*min(supp0),to=N*max(supp0))
+
+ d1 <- d(D1)(supp0); d1 <- c(d1,numeric((length(supp0)-1)*(N-1)))
+
+ ## computation of DFT
+ ftde1 <- fft(d1)
+
+ ## convolution theorem for DFTs
+ newd <- Re(fft(ftde1^N, inverse = TRUE)) / length(ftde1)
+ newd <- (newd >= .Machine$double.eps)*newd
+
+ rsum.u <- min( sum( cumsum(rev(newd))) > ep/2)+1, length(supp1))
+ rsum.l <- max( sum( cumsum(newd) < ep/2), 1)
+
+ newd <- newd[rsum.l:rsum.u]
+ newd <- newd/sum(newd)
+ supp1 <- supp1[rsum.l:rsum.u]
+
+ supp2 <- supp1[newd>ep]
+ newd2 <- newd[newd>ep]
+ newd2 <- newd2/sum(newd2)
+
+ Symmetry <- NoSymmetry()
+ if(is(D1 at Symmetry,"SphericalSymmetry"))
+ Symmetry <- SphericalSymmetry(N*SymmCenter(D1 at Symmetry))
+
+ if( length(supp1) >= 2 * length(supp2))
+ return(DiscreteDistribution(supp = supp2, prob = newd2,
+ .withArith = TRUE, Symmetry = Symmetry))
+ else
+ return(LatticeDistribution(supp = supp1, prob = newd,
+ .withArith = TRUE, Symmetry = Symmetry))
+})
+
+###############################################################################
+#
+# new from 2.0: convpov for AcDcLcDistribution
+#
+###############################################################################
+#
+setMethod("convpow",
+ signature(D1 = "AcDcLcDistribution"),
+ function(D1, N, ep = getdistrOption("TruncQuantile")){
+ if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if (N==1) return(D1)
+ e1 <- as(D1, "UnivarLebDecDistribution")
+ if(is(e1,"DiscreteDistribution")) return(convpow(e1,N))
+ if(is(e1,"AbscontDistribution")) return(convpow(e1,N))
+
+ if(!is.numeric(ep)) stop("argument 'ep' must be a numeric.")
+ if(length(ep)!=1) stop("argument 'ep' must be a numeric of length 1.")
+ if((ep<0)||(ep>1)) stop("argument 'ep' must be in (0,1).")
+
+ aw1 <- acWeight(e1)
+ dw1 <- 1-aw1
+ dD1 <- discretePart(e1)
+ aD1 <- acPart(e1)
+ dD1 <- discretePart(e1)
+ if(is(dD1,"LatticeDistribution"))
+ dD1 <- as(dD1,"LatticeDistribution")
+ # dDm <- max(d.discrete(e1)(support(e1)))*dw1
+
+ if(aw1<ep) return(convpow(dD1,N))
+ if(dw1<ep) return(convpow(aD1,N))
+
+ maxN <- ceiling(2*log(ep)/log(dw1))
+ Nm <- min(maxN,N)
+ Mm <- N%/%Nm
+ Rm <- N-Mm*Nm
+
+ sumM <- function(mm){
+ db <- dbinom(0:mm, size = mm, prob = aw1)
+ im <- (0:mm)[db>ep^2]
+ db <- db[db>ep^2]
+ db <- db/sum(db)
+ if(length(im)>1){
+ DList <- lapply(im,
+ function(x) {
+ S.a <- convpow(aD1, x)
+ S.d <- convpow(dD1, mm-x) #as(dD1,
+ # "DiscreteDistribution"), mm-x)
+ as(S.a+S.d,"UnivarLebDecDistribution")
+ })
+ erg <- do.call(flat.LCD, c(DList, alist(mixCoeff = db)))
+ }else{
+ DList <- as(convpow(aD1,im)+convpow(dD1,mm-im),"UnivarLebDecDistribution")
+ erg <- flat.LCD(DList, mixCoeff = 1)
+ }
+ return(erg)
+ }
+
+ erg <- sumM(Nm)
+ if(Mm>1) erg <- convpow(erg,Mm,ep=ep)
+ if(Rm>0) erg <- sumM(Rm)+ as(erg,"UnivarLebDecDistribution")
+ if(is(erg,"UnivarLebDecDistribution")) erg <- simplifyD(erg)
+
+ if(is(D1 at Symmetry,"SphericalSymmetry"))
+ erg at Symmetry <- SphericalSymmetry(N*SymmCenter(D1 at Symmetry))
+ return(erg)
+})
+#
+###############################################################################
+setMethod("convpow",
+ signature(D1 = "DiscreteDistribution"),
+ function(D1, N){
+ if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if (N==1) return(D1)
+ if (N==2) return(D1+D1)
+ DN1 <- convpow(D1,N%/%2)
+ DN1 <- DN1 + DN1
+ if (N%%2==1) DN1 <- DN1+D1
+ return(DN1)
+ })
+###############################################################################
+
+setMethod("convpow",
+ signature(D1 = "Norm"),
+ function(D1, N)
+ {if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if(N==1) D1 else Norm(mean = N*mean(D1), sd = sqrt(N)*sd(D1))}
+ )
+
+setMethod("convpow",
+ signature(D1 = "Pois"),
+ function(D1, N)
+ {if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if(N==1) D1 else Pois(lambda=N*lambda(D1))
+ }
+ )
+
+setMethod("convpow",
+ signature(D1 = "Binom"),
+ function(D1, N)
+ {if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if(N==1) D1 else Binom(size=N*size(D1),prob=prob(D1))}
+ )
+
+setMethod("convpow",
+ signature(D1 = "Nbinom"),
+ function(D1, N)
+ {if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if(N==1) D1 else Nbinom(size=N*size(D1),prob=prob(D1))}
+ )
+
+#setMethod("convpow",
+# signature(D1 = "Gammad"),
+# function(D1, N)
+# {if((N<1)||(abs(floor(N)-N)>.Machine$double.eps))
+# stop("N has to be a natural greater than or equal to 1")
+# if(N==1) D1 else Gammad(shape=N*shape(D1),scale=scale(D1))}
+# )
+
+setMethod("convpow",
+ signature(D1 = "Dirac"),
+ function(D1, N)
+ {if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ Dirac(shape=N*location(D1))}
+ )
+
+setMethod("convpow",
+ signature(D1 = "ExpOrGammaOrChisq"),
+ function(D1, N)
+ {if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if(N==1) return(D1)
+ else if(is(D1,"Gammad"))
+ {D1 <- as(D1,"Gammad")
+ return(Gammad(shape = N*shape(D1),
+ scale = scale(D1))) }
+ else convpow(as(D1, "AbscontDistribution"),N)}
+ )
+
+ setMethod("convpow",
+ signature(D1 = "Cauchy"),
+ function(D1, N)
+ {if( !.isNatural0(N))
+ stop("N has to be a natural (or 0)")
+ if (N==0) return(Dirac(0))
+ if(N==1) D1 else Cauchy(location = N*location(D1),
+ scale = N*scale(D1))}
+ )
Deleted: branches/distr-2.3/pkg/distr/R/Convpow0.R
===================================================================
--- branches/distr-2.3/pkg/distr/R/Convpow0.R 2010-05-16 22:38:48 UTC (rev 653)
+++ branches/distr-2.3/pkg/distr/R/Convpow0.R 2010-05-16 22:39:26 UTC (rev 654)
@@ -1,292 +0,0 @@
-##########################################################
-## Function for n-fold convolution
-## -- absolute continuous distribution --
-##########################################################
-
-##implentation of Algorithm 3.4. of
-# Kohl, M., Ruckdeschel, P., Stabla, T. (2005):
-# General purpose convolution algorithm for distributions
-# in S4-Classes by means of FFT.
-# Technical report, Feb. 2005. Also available in
-# http://www.uni-bayreuth.de/departments/math/org/mathe7/
-# /RUCKDESCHEL/pubs/comp.pdf
-
-
-setMethod("convpow",
- signature(D1 = "AbscontDistribution"),
- function(D1, N){
- if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if (N==1) return(D1)
- ##STEP 1
-
- lower <- getLow(D1); upper <- getUp(D1);
-
- ##STEP 2
-
- ## binary logarithm of the effective number of gridpoints
- m <- max(getdistrOption("DefaultNrFFTGridPointsExponent") -
- floor(log(N)/log(2)),5)
- M <- 2^m
- Nl <-2^ceiling(log(N)/log(2))
-
- h <- (upper-lower)/M
- dp1 <- .discretizeP(D1, lower, upper, h)
-
- ##STEP 3
-
- dpn0 <- c(dp1, numeric((Nl-1)*M))
- ##STEP 4
-
- ## computation of DFT
- fftdpn <- fft(dpn0)
-
- ##STEP 5
-
- ## convolution theorem for DFTs
- dpn <- c(0,(Re(fft(fftdpn^N, inverse = TRUE)) / (Nl*M))[1:(N*M-N+2)])
-
- x <- seq(from = N*lower+N/2*h, to = N*upper-N/2*h, by = h)
- x <- c(x[1]-h, x[1], x+h)
-
- ## density (steps 5--7)
-
- dfun <- .makeDNew(x, dpn, h)
-
- ## cdf (steps 5--7)
- pfun <- .makePNew(x, dpn, h, .notwithLArg(D1))
-
- ## continuity correction by h/2
-
- ## quantile function
- yL <- if (q(D1)(0) == -Inf) -Inf else N*lower
- yR <- if (q(D1)(1) == Inf) Inf else N*upper
- px.l <- pfun(x + 0.5*h)
- px.u <- pfun(x + 0.5*h, lower.tail = FALSE)
- qfun <- .makeQNew(x + 0.5*h, px.l, px.u, .notwithLArg(D1), yL, yR)
-
- rfun = function(n) colSums(matrix(r(D1)(n*N), ncol=n))
-
- object <- new("AbscontDistribution", r = rfun, d = dfun, p = pfun,
- q = qfun, .withArith = TRUE, .withSim = FALSE)
-
- if(is(D1 at Symmetry,"SphericalSymmetry"))
- object at Symmetry <- SphericalSymmetry(N*SymmCenter(D1 at Symmetry))
-
- rm(m, dpn, dp1, dpn0, fftdpn)
- rm(h, px.u, px.l, rfun, dfun, qfun, pfun, upper, lower)
- return(object)
-})
-
-setMethod("convpow",
- signature(D1 = "LatticeDistribution"),
- function(D1, N, ep = getdistrOption("TruncQuantile")){
- if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if (N==1) return(D1)
-
- if(!is.numeric(ep)) stop("argument 'ep' must be a numeric.")
- if(length(ep)!=1) stop("argument 'ep' must be a numeric of length 1.")
- if((ep<0)||(ep>1)) stop("argument 'ep' must be in (0,1).")
-
- w <- width(lattice(D1))
-
- supp0 <- support(D1)
- supp1 <- seq(by=abs(w),from=N*min(supp0),to=N*max(supp0))
-
- d1 <- d(D1)(supp0); d1 <- c(d1,numeric((length(supp0)-1)*(N-1)))
-
- ## computation of DFT
- ftde1 <- fft(d1)
-
- ## convolution theorem for DFTs
- newd <- Re(fft(ftde1^N, inverse = TRUE)) / length(ftde1)
- newd <- (newd >= .Machine$double.eps)*newd
-
- rsum.u <- min( sum( cumsum(rev(newd))) > ep/2)+1, length(supp1))
- rsum.l <- max( sum( cumsum(newd) < ep/2), 1)
-
- newd <- newd[rsum.l:rsum.u]
- newd <- newd/sum(newd)
- supp1 <- supp1[rsum.l:rsum.u]
-
- supp2 <- supp1[newd>ep]
- newd2 <- newd[newd>ep]
- newd2 <- newd2/sum(newd2)
-
- Symmetry <- NoSymmetry()
- if(is(D1 at Symmetry,"SphericalSymmetry"))
- Symmetry <- SphericalSymmetry(N*SymmCenter(D1 at Symmetry))
-
- if( length(supp1) >= 2 * length(supp2))
- return(DiscreteDistribution(supp = supp2, prob = newd2,
- .withArith = TRUE, Symmetry = Symmetry))
- else
- return(LatticeDistribution(supp = supp1, prob = newd,
- .withArith = TRUE, Symmetry = Symmetry))
-})
-
-###############################################################################
-#
-# new from 2.0: convpov for AcDcLcDistribution
-#
-###############################################################################
-#
-setMethod("convpow",
- signature(D1 = "AcDcLcDistribution"),
- function(D1, N, ep = getdistrOption("TruncQuantile")){
- if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if (N==1) return(D1)
- e1 <- as(D1, "UnivarLebDecDistribution")
- if(is(e1,"DiscreteDistribution")) return(convpow(e1,N))
- if(is(e1,"AbscontDistribution")) return(convpow(e1,N))
-
- if(!is.numeric(ep)) stop("argument 'ep' must be a numeric.")
- if(length(ep)!=1) stop("argument 'ep' must be a numeric of length 1.")
- if((ep<0)||(ep>1)) stop("argument 'ep' must be in (0,1).")
-
- aw1 <- acWeight(e1)
- dw1 <- 1-aw1
- dD1 <- discretePart(e1)
- aD1 <- acPart(e1)
- dD1 <- discretePart(e1)
- if(is(dD1,"LatticeDistribution"))
- dD1 <- as(dD1,"LatticeDistribution")
- # dDm <- max(d.discrete(e1)(support(e1)))*dw1
-
- if(aw1<ep) return(convpow(dD1,N))
- if(dw1<ep) return(convpow(aD1,N))
-
- maxN <- ceiling(2*log(ep)/log(dw1))
- Nm <- min(maxN,N)
- Mm <- N%/%Nm
- Rm <- N-Mm*Nm
-
- sumM <- function(mm){
- db <- dbinom(0:mm, size = mm, prob = aw1)
- im <- (0:mm)[db>ep^2]
- db <- db[db>ep^2]
- db <- db/sum(db)
- if(length(im)>1){
- DList <- lapply(im,
- function(x) {
- S.a <- convpow(aD1, x)
- S.d <- convpow(dD1, mm-x) #as(dD1,
- # "DiscreteDistribution"), mm-x)
- as(S.a+S.d,"UnivarLebDecDistribution")
- })
- erg <- do.call(flat.LCD, c(DList, alist(mixCoeff = db)))
- }else{
- DList <- as(convpow(aD1,im)+convpow(dD1,mm-im),"UnivarLebDecDistribution")
- erg <- flat.LCD(DList, mixCoeff = 1)
- }
- return(erg)
- }
-
- erg <- sumM(Nm)
- if(Mm>1) erg <- convpow(erg,Mm,ep=ep)
- if(Rm>0) erg <- sumM(Rm)+ as(erg,"UnivarLebDecDistribution")
- if(is(erg,"UnivarLebDecDistribution")) erg <- simplifyD(erg)
-
- if(is(D1 at Symmetry,"SphericalSymmetry"))
- erg at Symmetry <- SphericalSymmetry(N*SymmCenter(D1 at Symmetry))
- return(erg)
-})
-#
-###############################################################################
-setMethod("convpow",
- signature(D1 = "DiscreteDistribution"),
- function(D1, N){
- if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if (N==1) return(D1)
- if (N==2) return(D1+D1)
- DN1 <- convpow(D1,N%/%2)
- DN1 <- DN1 + DN1
- if (N%%2==1) DN1 <- DN1+D1
- return(DN1)
- })
-###############################################################################
-
-setMethod("convpow",
- signature(D1 = "Norm"),
- function(D1, N)
- {if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if(N==1) D1 else Norm(mean = N*mean(D1), sd = sqrt(N)*sd(D1))}
- )
-
-setMethod("convpow",
- signature(D1 = "Pois"),
- function(D1, N)
- {if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if(N==1) D1 else Pois(lambda=N*lambda(D1))
- }
- )
-
-setMethod("convpow",
- signature(D1 = "Binom"),
- function(D1, N)
- {if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if(N==1) D1 else Binom(size=N*size(D1),prob=prob(D1))}
- )
-
-setMethod("convpow",
- signature(D1 = "Nbinom"),
- function(D1, N)
- {if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if(N==1) D1 else Nbinom(size=N*size(D1),prob=prob(D1))}
- )
-
-#setMethod("convpow",
-# signature(D1 = "Gammad"),
-# function(D1, N)
-# {if((N<1)||(abs(floor(N)-N)>.Machine$double.eps))
-# stop("N has to be a natural greater than or equal to 1")
-# if(N==1) D1 else Gammad(shape=N*shape(D1),scale=scale(D1))}
-# )
-
-setMethod("convpow",
- signature(D1 = "Dirac"),
- function(D1, N)
- {if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- Dirac(shape=N*location(D1))}
- )
-
-setMethod("convpow",
- signature(D1 = "ExpOrGammaOrChisq"),
- function(D1, N)
- {if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if(N==1) return(D1)
- else if(is(D1,"Gammad"))
- {D1 <- as(D1,"Gammad")
- return(Gammad(shape = N*shape(D1),
- scale = scale(D1))) }
- else convpow(as(D1, "AbscontDistribution"),N)}
- )
-
- setMethod("convpow",
- signature(D1 = "Cauchy"),
- function(D1, N)
- {if( !.isNatural0(N))
- stop("N has to be a natural (or 0)")
- if (N==0) return(Dirac(0))
- if(N==1) D1 else Cauchy(location = N*location(D1),
- scale = N*scale(D1))}
- )
More information about the Distr-commits
mailing list