[Distr-commits] r653 - 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:38:48 CEST 2010
Author: ruckdeschel
Date: 2010-05-17 00:38:48 +0200 (Mon, 17 May 2010)
New Revision: 653
Removed:
branches/distr-2.3/pkg/distr/R/Convpow.r
Log:
a renaming schedule : enforcing ending .R in convpow.R
Deleted: branches/distr-2.3/pkg/distr/R/Convpow.r
===================================================================
--- branches/distr-2.3/pkg/distr/R/Convpow.r 2010-05-16 22:37:17 UTC (rev 652)
+++ branches/distr-2.3/pkg/distr/R/Convpow.r 2010-05-16 22:38:48 UTC (rev 653)
@@ -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 <- (abs(newd) >= .Machine$double.eps)*newd
-
- rsum.u <- min( sum( rev(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