[Distr-commits] r694 - in pkg: SweaveListingUtils SweaveListingUtils/inst SweaveListingUtils/man distr distr/R distr/inst distr/inst/doc distr/man distrDoc distrDoc/inst distrDoc/man distrEllipse distrEllipse/inst distrEllipse/man distrEx distrEx/R distrEx/chm distrEx/inst distrEx/man distrEx/src distrMod distrMod/R distrMod/inst distrMod/man distrSim distrSim/R distrSim/inst distrSim/man distrTEst distrTEst/R distrTEst/inst distrTEst/man distrTeach distrTeach/R distrTeach/inst distrTeach/man startupmsg utils
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Dec 2 14:41:24 CET 2010
Author: ruckdeschel
Date: 2010-12-02 14:41:24 +0100 (Thu, 02 Dec 2010)
New Revision: 694
Added:
pkg/distr/R/Convpow.R
pkg/distr/R/igamma.R
pkg/distr/man/igamma.Rd
pkg/distrEx/R/GEV.R
pkg/distrEx/R/kMAD.R
pkg/distrEx/man/GEV-class.Rd
pkg/distrEx/man/GEV.Rd
pkg/distrEx/man/GEVParameter-class.Rd
pkg/distrEx/man/kMAD.Rd
pkg/distrEx/src/kMad.c
pkg/distrMod/man/NBinomFamily.Rd
pkg/utils/getRevNr.R
Removed:
pkg/distr/R/Convpow.r
Modified:
pkg/SweaveListingUtils/DESCRIPTION
pkg/SweaveListingUtils/inst/NEWS
pkg/SweaveListingUtils/man/0SweaveListingUtils-package.Rd
pkg/distr/DESCRIPTION
pkg/distr/NAMESPACE
pkg/distr/R/AllGenerics.R
pkg/distr/R/ContDistribution.R
pkg/distr/R/DiscreteDistribution.R
pkg/distr/R/NegbinomDistribution.R
pkg/distr/R/UnivarLebDecDistribution.R
pkg/distr/R/UtilitiesDistributions.R
pkg/distr/R/internalUtils.R
pkg/distr/R/internals-qqplot.R
pkg/distr/R/plot-methods.R
pkg/distr/R/plot-methods_LebDec.R
pkg/distr/R/qqbounds.R
pkg/distr/R/qqplot.R
pkg/distr/inst/NEWS
pkg/distr/inst/doc/newDistributions.Rnw
pkg/distr/man/0distr-package.Rd
pkg/distr/man/AbscontDistribution-class.Rd
pkg/distr/man/DiscreteDistribution-class.Rd
pkg/distr/man/Math-methods.Rd
pkg/distr/man/RtoDPQ.LC.Rd
pkg/distr/man/RtoDPQ.Rd
pkg/distr/man/UnivarLebDecDistribution.Rd
pkg/distr/man/internals-qqplot.Rd
pkg/distr/man/qqplot.Rd
pkg/distrDoc/DESCRIPTION
pkg/distrDoc/inst/NEWS
pkg/distrDoc/man/0distrDoc-package.Rd
pkg/distrEllipse/DESCRIPTION
pkg/distrEllipse/inst/NEWS
pkg/distrEllipse/man/0distrEllipse-package.Rd
pkg/distrEx/DESCRIPTION
pkg/distrEx/NAMESPACE
pkg/distrEx/R/AllClass.R
pkg/distrEx/R/AllGeneric.R
pkg/distrEx/R/AllInitialize.R
pkg/distrEx/R/CvMDist.R
pkg/distrEx/R/Expectation.R
pkg/distrEx/R/Functionals.R
pkg/distrEx/R/Kurtosis.R
pkg/distrEx/R/Skewness.R
pkg/distrEx/chm/E.html
pkg/distrEx/inst/NEWS
pkg/distrEx/man/0distrEx-package.Rd
pkg/distrEx/man/E.Rd
pkg/distrEx/man/Var.Rd
pkg/distrEx/src/distrEx.dll
pkg/distrMod/DESCRIPTION
pkg/distrMod/NAMESPACE
pkg/distrMod/R/AllClass.R
pkg/distrMod/R/AllPlot.R
pkg/distrMod/R/AllReturnClasses.R
pkg/distrMod/R/SimpleL2ParamFamilies.R
pkg/distrMod/R/qqplot.R
pkg/distrMod/inst/NEWS
pkg/distrMod/man/0distrMod-package.Rd
pkg/distrMod/man/qqplot.Rd
pkg/distrSim/DESCRIPTION
pkg/distrSim/R/plot-methods.R
pkg/distrSim/inst/NEWS
pkg/distrSim/man/0distrSim-package.Rd
pkg/distrTEst/DESCRIPTION
pkg/distrTEst/R/plot-methods.R
pkg/distrTEst/inst/NEWS
pkg/distrTEst/man/0distrTEst-package.Rd
pkg/distrTeach/DESCRIPTION
pkg/distrTeach/R/illustCLT.R
pkg/distrTeach/R/illustLLN.R
pkg/distrTeach/inst/NEWS
pkg/distrTeach/man/0distrTeach-package.Rd
pkg/startupmsg/DESCRIPTION
pkg/utils/DESCRIPTIONutils.R
pkg/utils/NEWS
pkg/utils/README-R-utils
pkg/utils/setNewVersion.r
pkg/utils/showAllFiles.R
pkg/utils/showsvnlog.R
Log:
merged versions 2.3 into trunk
Modified: pkg/SweaveListingUtils/DESCRIPTION
===================================================================
--- pkg/SweaveListingUtils/DESCRIPTION 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/SweaveListingUtils/DESCRIPTION 2010-12-02 13:41:24 UTC (rev 694)
@@ -1,17 +1,21 @@
Package: SweaveListingUtils
Title: Utilities for Sweave together with TeX listings package
Encoding: latin1
-Version: 0.4.5
+Version: 0.5
Depends: R(>= 2.10.0), startupmsg
Suggests: distr
Imports: stats
LazyLoad: yes
Author: Peter Ruckdeschel
-Description: provides utilities for defining R / Rd as Tex-package-listings "language" and including R / Rd source file (sniplets) copied
- from R-forge in its most recent version (or another url) thereby avoiding inconsistencies between vignette and documented
- source code
+Description: provides utilities for defining R / Rd as
+ Tex-package-listings "language" and including R / Rd source
+ file (sniplets) copied from R-forge in its most recent version
+ (or another url) thereby avoiding inconsistencies between
+ vignette and documented source code
Maintainer: Peter Ruckdeschel <Peter.Ruckdeschel at itwm.fraunhofer.de>
License: LGPL-3
-Date: 2010-05-15
-LastChangedDate: {$LastChangedDate$}
+Date: 1
+LastChangedDate: {$LastChangedDate: 2010-05-17 01:42:46 +0200 (Mo, 17
+ Mai 2010) $}
LastChangedRevision: {$LastChangedRevision$}
+SVNRevision: 693
Modified: pkg/SweaveListingUtils/inst/NEWS
===================================================================
--- pkg/SweaveListingUtils/inst/NEWS 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/SweaveListingUtils/inst/NEWS 2010-12-02 13:41:24 UTC (rev 694)
@@ -3,6 +3,22 @@
######################################################################
##############
+v 0.5
+##############
+
+under the hood:
++ followed change in naming convention on r-forge,ie
+ base.url = "http://r-forge.r-project.org/scm/viewvc.php/*checkout*/pkg/"
+ instead of
+ base.url = "http://r-forge.r-project.org/plugins/scmsvn/viewcvs.php/*checkout*/pkg/"
++ quotes are now \texttt in listings by default (as required in jss)
+ --- fixes a quote - typesetting issue as required by JSS
++ DESCRIPTION files and package-help files gain a tag
+ SVNRevision to be filled by get[All]RevNr.R from utils in distr
+
+
+
+##############
v 0.4
##############
Modified: pkg/SweaveListingUtils/man/0SweaveListingUtils-package.Rd
===================================================================
--- pkg/SweaveListingUtils/man/0SweaveListingUtils-package.Rd 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/SweaveListingUtils/man/0SweaveListingUtils-package.Rd 2010-12-02 13:41:24 UTC (rev 694)
@@ -15,11 +15,12 @@
\details{
\tabular{ll}{
Package: \tab SweaveListingUtils \cr
-Version: \tab 0.4.4 \cr
-Date: \tab 2010-04-19 \cr
+Version: \tab 0.5 \cr
+Date: \tab 1 \cr
Depends: \tab R(>= 2.10.0), startupmsg \cr
LazyLoad: \tab yes \cr
License: \tab LGPL-3 \cr
+SVNRevision: \tab 693 \cr
}
TeX-package \file{listings}, confer \url{http://www.ctan.org/tex-archive/macros/latex/contrib/listings/},
Modified: pkg/distr/DESCRIPTION
===================================================================
--- pkg/distr/DESCRIPTION 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/distr/DESCRIPTION 2010-12-02 13:41:24 UTC (rev 694)
@@ -1,16 +1,20 @@
Package: distr
-Version: 2.2.3
-Date: 2010-06-15
+Version: 2.3
+Date: 2010-12-02
Title: Object oriented implementation of distributions
-Description: Object oriented implementation of distributions
-Author: Florian Camphausen, Matthias Kohl, Peter Ruckdeschel, Thomas Stabla
+Description: S4 Classes and Methods for distributions
+Author: Florian Camphausen, Matthias Kohl, Peter Ruckdeschel, Thomas
+ Stabla
Maintainer: Peter Ruckdeschel <Peter.Ruckdeschel at itwm.fraunhofer.de>
-Depends: R(>= 2.2.0), methods, graphics, startupmsg, sfsmisc
+Depends: R(>= 2.2.0), methods, graphics, startupmsg, sfsmisc,
+ SweaveListingUtils
Suggests: distrEx, SweaveListingUtils
Imports: stats
LazyLoad: yes
Encoding: latin1
License: LGPL-3
URL: http://distr.r-forge.r-project.org/
-LastChangedDate: {$LastChangedDate$}
+LastChangedDate: {$LastChangedDate: 2010-06-15 20:51:07 +0200 (Di, 15
+ Jun 2010) $}
LastChangedRevision: {$LastChangedRevision$}
+SVNRevision: 693
Modified: pkg/distr/NAMESPACE
===================================================================
--- pkg/distr/NAMESPACE 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/distr/NAMESPACE 2010-12-02 13:41:24 UTC (rev 694)
@@ -64,7 +64,7 @@
"shape1", "shape1<-", "shape2", "shape2<-",
"size", "size<-", "support", "initialize",
"print", "plot", "+", "-", "/", "*", "coerce",
- "Math", "log", "log10", "gamma", "lgamma",
+ "Math", "log", "log10", "gamma", "lgamma", "digamma",
"dim", "show", "convpow", "pivot", "sign",
"lattice", "width", "Length", "pivot<-",
"width<-", "Length<-", "liesInSupport",
@@ -91,5 +91,5 @@
export("PosDefSymmMatrix","PosSemDefSymmMatrix")
export("NoSymmetry", "EllipticalSymmetry", "SphericalSymmetry",
"DistrSymmList")
-export("qqbounds")
+export("qqbounds","igamma")
exportMethods("qqplot")
Modified: pkg/distr/R/AllGenerics.R
===================================================================
--- pkg/distr/R/AllGenerics.R 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/distr/R/AllGenerics.R 2010-12-02 13:41:24 UTC (rev 694)
@@ -320,7 +320,7 @@
if(!isGeneric("setgaps"))
setGeneric("setgaps", function(object, ...) standardGeneric("setgaps"))
-#### generics for log, log10, lgamma, gamma
+#### generics for log, log10, lgamma, gamma, digamma
if(getRversion()<'2.9.0'){
@@ -330,6 +330,8 @@
setGeneric("log10")
if(!isGeneric("lgamma"))
setGeneric("lgamma")
+if(!isGeneric("digamma"))
+ setGeneric("digamma")
if(!isGeneric("gamma"))
setGeneric("gamma")
if(!isGeneric("sign"))
Modified: pkg/distr/R/ContDistribution.R
===================================================================
--- pkg/distr/R/ContDistribution.R 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/distr/R/ContDistribution.R 2010-12-02 13:41:24 UTC (rev 694)
@@ -418,10 +418,19 @@
## Group Math for absolutly continuous distributions
setMethod("Math", "AbscontDistribution",
- function(x){ rnew <- function(n, ...){}
+ function(x){
+
+ rnew <- function(n, ...){}
body(rnew) <- substitute({ f(g(n, ...)) },
list(f = as.name(.Generic), g = x at r))
- object <- AbscontDistribution( r = rnew,
+
+ n <- 10^getdistrOption("RtoDPQ.e")+1
+ u <- seq(0,1,length=n+1); u <- (u[1:n]+u[2:(n+1)])/2
+ y <- callGeneric(q(x)(u))
+ DPQnew <- RtoDPQ(r=rnew, y=y)
+
+ object <- AbscontDistribution(d = DPQnew$d, p = DPQnew$p,
+ r = rnew, q = DPQnew$q,
.withSim = TRUE, .withArith = TRUE)
object
})
@@ -587,23 +596,82 @@
+setMethod("digamma", "AbscontDistribution",
+ function(x){
+ rnew <- function(n, ...){}
+ body(rnew) <- substitute({ digamma(g(n, ...)) }, list(g = x at r))
+ px0 <- p(x)(0)
+ if(px0>0) stop("argument of 'digamma' must be concentrated on positive values")
+ xx <- x
+
+ pnew <- function(q, lower.tail = TRUE, log.p = FALSE){
+ iq <- igamma(q)
+ px <- p(xx)(iq, lower.tail = lower.tail, log.p = log.p)
+ return(px)
+ }
+ dnew <- function(x, log = FALSE){
+ ix <- igamma(x)
+ dx <- d(xx)(ix, log = log)
+ nx <- trigamma(ix)
+ if(log) dx <- dx - log(nx)
+ else dx <- dx/nx
+ return(dx)
+ }
+
+ .x <- sort(c(qexp(unique(pmin(seq(0,1,length=5e4)+1e-10,1-1e-10))),
+ -abs(rnorm(1e4)),
+ qcauchy(seq(0.999,1-1e-10,length=5e3),lower.tail=FALSE)))
+ i <- 0; x0 <- 1
+ while(pnew(x0,lower.tail = FALSE)> getdistrOption("TruncQuantile") && i < 20)
+ x0 <- x0 * 2
+ up1 <- x0
+ i <- 0; x0 <- -1
+ while(pnew(x0)> getdistrOption("TruncQuantile") && i < 20)
+ x0 <- x0 * 2
+ low1 <- x0
+
+
+
+ qnew <- .P2Q(p = pnew, xx =.x,
+ ql = low1, qu=up1,
+ ngrid = getdistrOption("DefaultNrGridPoints"),
+ qL = -Inf, qU = Inf)
+
+
+ object <- AbscontDistribution( r = rnew, d = dnew, p = pnew, q=qnew,
+ .withSim = TRUE, .withArith = TRUE, .logExact = FALSE)
+ object
+ })
+
setMethod("lgamma", "AbscontDistribution",
function(x){
- rnew = function(n, ...){}
+ rnew <- function(n, ...){}
body(rnew) <- substitute({ lgamma(g(n, ...)) }, list(g = x at r))
- object <- AbscontDistribution( r = rnew,
- .withSim = TRUE, .withArith = TRUE)
+
+ n <- 10^getdistrOption("RtoDPQ.e")+1
+ u <- seq(0,1,length=n+1); u <- (u[1:n]+u[2:(n+1)])/2
+ y <- lgamma(q(x)(u))
+ DPQnew <- RtoDPQ(r=rnew, y=y)
+
+ object <- AbscontDistribution( r = rnew, d = DPQnew$d, p = DPQnew$p,
+ q=DPQnew$q, .withSim = TRUE, .withArith = TRUE)
object
})
setMethod("gamma", "AbscontDistribution",
function(x){
- rnew = function(n, ...){}
+ rnew <- function(n, ...){}
body(rnew) <- substitute({ gamma(g(n, ...)) }, list(g = x at r))
- object <- AbscontDistribution( r = rnew,
- .withSim = TRUE, .withArith = TRUE)
+ n <- 10^getdistrOption("RtoDPQ.e")+1
+ u <- seq(0,1,length=n+1); u <- (u[1:n]+u[2:(n+1)])/2
+ y <- gamma(q(x)(u))
+ DPQnew <- RtoDPQ(r=rnew, y=y)
+
+ object <- AbscontDistribution( r = rnew, d = DPQnew$d, p = DPQnew$p,
+ q=DPQnew$q, .withSim = TRUE, .withArith = TRUE)
object
})
+
setMethod("sqrt", "AbscontDistribution",
function(x) x^0.5)
Copied: pkg/distr/R/Convpow.R (from rev 693, branches/distr-2.3/pkg/distr/R/Convpow.R)
===================================================================
--- pkg/distr/R/Convpow.R (rev 0)
+++ pkg/distr/R/Convpow.R 2010-12-02 13:41:24 UTC (rev 694)
@@ -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: pkg/distr/R/Convpow.r
===================================================================
--- pkg/distr/R/Convpow.r 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/distr/R/Convpow.r 2010-12-02 13:41:24 UTC (rev 694)
@@ -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))}
- )
Modified: pkg/distr/R/DiscreteDistribution.R
===================================================================
--- pkg/distr/R/DiscreteDistribution.R 2010-12-02 13:14:26 UTC (rev 693)
+++ pkg/distr/R/DiscreteDistribution.R 2010-12-02 13:41:24 UTC (rev 694)
@@ -519,6 +519,19 @@
})
+setMethod("digamma", "DiscreteDistribution",
+ function(x){
+ px0 <- p(x)(0)
+ if(px0>0) stop("argument of 'digamma' must be concentrated on positive values")
+ rnew <- function(n, ...){}
+ body(rnew) <- substitute({ digamma(g(n, ...)) }, list(g = x at r))
+
+ object <- DiscreteDistribution(
+ supp=digamma(support(x)),
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/distr -r 694
More information about the Distr-commits
mailing list