[Distr-commits] r183 - branches/distr-2.0/pkg/distrSim/R pkg/distr/R pkg/distr/chm pkg/distr/man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Wed Jul 9 15:56:14 CEST 2008

Author: ruckdeschel
Date: 2008-07-09 15:56:14 +0200 (Wed, 09 Jul 2008)
New Revision: 183

Modified:
branches/distr-2.0/pkg/distrSim/R/subsetting-methods.R
pkg/distr/R/99.R
pkg/distr/R/Convpow.r
pkg/distr/R/DiscreteDistribution.R
pkg/distr/R/LatticeDistribution.R
pkg/distr/R/internalUtils.R
pkg/distr/R/plot-methods.R
pkg/distr/R/plot-methods_LebDec.R
pkg/distr/chm/Distr.chm
pkg/distr/chm/Distr.toc
pkg/distr/chm/distroptions.html
pkg/distr/chm/internals.html
pkg/distr/man/distroptions.Rd
pkg/distr/man/internals.Rd
Log:
just in "trunc" (not in 2.0 "branch"):

enhanced version of convpow for DiscreteDistributions
(needs atmost 2 log_2(N) convolutions) --- without recursions!
(corresponding enhancements for plotting routines)

new options:
\item{\code{DistrCollapse}}{logical; in convolving discrete distributions, shall support points
with distance smaller than \code{DistrResolution} be collapsed; default value:  \code{TRUE}}
\item{\code{DistrConvPoints}}{numeric; maximal number of support points in convolving discrete distributions;
if \code{DistrCollapse = TRUE} each summand is discretized to this number of support points; default value:  \code{700}}
\item{\code{DistrMaxPlotPoints}}{numeric; maximal number of support points of discrete distributions to be plotted;
if support is larger, then for plotting it is thinned out; default value:  \code{200}}

Modified: branches/distr-2.0/pkg/distrSim/R/subsetting-methods.R
===================================================================
--- branches/distr-2.0/pkg/distrSim/R/subsetting-methods.R	2008-07-01 18:15:28 UTC (rev 182)
+++ branches/distr-2.0/pkg/distrSim/R/subsetting-methods.R	2008-07-09 13:56:14 UTC (rev 183)
@@ -58,6 +58,8 @@
## Is the following line correct?
## should it be: if(!is.list(i)) i <- lapply(kl0, function(y) y)
## or: if(!is.list(i)) i <- as.list(kl0)?
+
+
if(!is.list(i)) i <- lapply(kl0, function(y) i)

if(is(value,"atomic"))

Modified: pkg/distr/R/99.R
===================================================================
--- pkg/distr/R/99.R	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/R/99.R	2008-07-09 13:56:14 UTC (rev 183)
@@ -10,7 +10,9 @@
## new Items from 2.0:
withgaps = TRUE,
simplifyD = TRUE,
-                      DistrCollapse = TRUE
+                      DistrCollapse = TRUE,
+                      DistrConvPoints = 700,
+                      DistrMaxPlotPoints = 200
)

.OkTyp <-  c("DiscreteDistribution","AbscontDistribution",

Modified: pkg/distr/R/Convpow.r
===================================================================
--- pkg/distr/R/Convpow.r	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/R/Convpow.r	2008-07-09 13:56:14 UTC (rev 183)
@@ -112,6 +112,48 @@
return(LatticeDistribution(supp=supp1,prob=newd))
})

+setMethod("convpow",
+          signature(D1 = "DiscreteDistribution"),
+          function(D1, N){
+            if (N==0) return(Dirac(0))
+            if (N==1) return(D1)
+            if (N==2) return(D1+D1)
+            #
+            l2N <- floor(log(N,2))
+
+            binrep <- numeric(l2N+1)
+            Dlist <- vector("list", l2N+1)
+            N0 <- N
+            i <- 1
+            j <- 0
+
+            invbinrep <- numeric(l2N+1)
+            while(N0 >1) {
+               binrep[i] <- N0 %%2
+               if(binrep[i]) {
+                  j<-j+1
+                  invbinrep[j] <- i
+                  }
+               N0 <- N0%/%2
+               i <- i+1}
+
+            invbinrep <- invbinrep[1:j]
+
+            D0 <- Dlist[[1]] <- D1
+            for( i in 2 : (l2N+1))
+               { Dlist[[i]] <- D0 + D0
+                 D0 <- Dlist[[i]]
+               }
+
+            D0 <- Dlist[[invbinrep[1]]]
+            if( j > 1 ){
+                for (i in 2 : j)
+                     D0 <- D0 + Dlist[[invbinrep[j]]]
+               }
+            return(D0)
+})
+
+
###############################################################################
#
# new from 2.0: convpov for  AcDcLcDistribution

Modified: pkg/distr/R/DiscreteDistribution.R
===================================================================
--- pkg/distr/R/DiscreteDistribution.R	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/R/DiscreteDistribution.R	2008-07-09 13:56:14 UTC (rev 183)
@@ -194,6 +194,18 @@

setMethod("+", c("DiscreteDistribution","DiscreteDistribution"),
function(e1,e2){
+            if(getdistrOption("DistrCollapse"))
+               {
+                bothdisced <- TRUE
+                if(length(support(e1)) > getdistrOption("DistrConvPoints"))
+                     e1 <- .discretizeD(e1, getdistrOption("DistrConvPoints"))
+                else bothdisced <- FALSE
+                if(length(support(e2)) > getdistrOption("DistrConvPoints"))
+                     e2 <- .discretizeD(e2, getdistrOption("DistrConvPoints"))
+                else bothdisced <- FALSE
+                if(bothdisced) return(e1+e2)
+               }
+
convolutedsupport <- rep(support(e1), each = length(support(e2))) +
support(e2)

@@ -203,7 +215,7 @@
rm(gridvalues1,gridvalues2)

tmptable <- data.frame(x = convolutedsupport, dx = convolutedvalues)
-            rm(convolutedsupport,convolutedvalues)
+            rm(convolutedsupport, convolutedvalues)
tmp <- tapply(tmptable$dx, tmptable$x, sum)
rm(tmptable)

Modified: pkg/distr/R/LatticeDistribution.R
===================================================================
--- pkg/distr/R/LatticeDistribution.R	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/R/LatticeDistribution.R	2008-07-09 13:56:14 UTC (rev 183)
@@ -177,6 +177,18 @@

setMethod("+", c("LatticeDistribution", "LatticeDistribution"),
function(e1,e2){
+
+            ### Step 0 Discretizing if too many grid points:
+
+            if(getdistrOption("DistrCollapse"))
+               {
+                if(length(support(e1)) > getdistrOption("DistrConvPoints"))
+                     e1 <- .discretizeD(e1, getdistrOption("DistrConvPoints"))
+                if(length(support(e2)) > getdistrOption("DistrConvPoints"))
+                     e2 <- .discretizeD(e2, getdistrOption("DistrConvPoints"))
+               }
+
+
### Step 1

e1 <- as(e1, "LatticeDistribution")

Modified: pkg/distr/R/internalUtils.R
===================================================================
--- pkg/distr/R/internalUtils.R	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/R/internalUtils.R	2008-07-09 13:56:14 UTC (rev 183)
@@ -189,8 +189,20 @@
}
}

+.discretizeD <- function(D, n)
+ {upper <- rev(support(D))[1]
+  lower <- support(D)[1]
+  h <- (upper-lower)/n
+  dp <- .discretizeP(D, lower, upper, h)
+  s0 <- seq(lower, upper, by = h)
+  s <- (s0[1:n]+s0[2:(n+1)])/2
+  D0 <- DiscreteDistribution(supp = s, prob = dp,
+             .withArith = D at .withArith, .withSim = D at .withSim)
+  LatticeDistribution(DiscreteDistribution=D0)
+  }

+
#------------------------------------------------------------------------------
# .makeD, .makeP, .makeQ
#------------------------------------------------------------------------------

Modified: pkg/distr/R/plot-methods.R
===================================================================
--- pkg/distr/R/plot-methods.R	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/R/plot-methods.R	2008-07-09 13:56:14 UTC (rev 183)
@@ -399,6 +399,19 @@

dx <- d(x)(supp)

+       ### possibly thinning out the support
+       thin <- FALSE
+       ngrid <- length(supp)
+       ngrid0 <- getdistrOption("DistrMaxPlotPoints")
+       if (ngrid > ngrid0)
+            {nn <- seq(1,ngrid, by = ngrid %/% ngrid0)
+             thin <- TRUE
+            }else nn <- seq(ngrid)
+
+       ngrid0 <- length(nn)
+       suppn <- supp[nn]
+
+
if(hasArg(ylim))
{if(length(ylim) != 2)
stop("Wrong length of Argument ylim")
@@ -421,24 +434,25 @@
}

owarn <- getOption("warn"); options(warn = -1)
-       do.call(plot, c(list(x = supp, dx, type = "h", pch = pch.a,
+       do.call(plot, c(list(x = suppn, dx[nn], type = "h", pch = pch.a,
ylim = ylim1, xlim=xlim, ylab = "d(x)", xlab = "x",
log = logpd), dots.without.pch))
options(warn = owarn)

+       if (thin) text(suppn[1],0.3, adj=0, "[grid thinned\n out]",
+                               cex = 0.7)

title(main = inner.d, line = lineT, cex.main = cex.inner,
col.main = col.inner)

if(do.points)
-          do.call(points, c(list(x = supp, y = dx, pch = pch.a,
+          do.call(points, c(list(x = suppn, y = dx[nn], pch = pch.a,
cex = cex.points, col = col.points), dots.for.points))

owarn <- getOption("warn"); options(warn = -1)

-       ngrid <- length(supp)

-       supp1 <- if(ngrid>1) supp else c(-max(1,abs(supp))*.08,0)+supp
+       supp1 <- if(ngrid0>1) suppn else c(-max(1,abs(suppn))*.08,0)+suppn
psupp1 <- c(0,p(x)(supp1))

do.call(plot, c(list(x = stepfun(x = supp1, y = psupp1),
@@ -448,15 +462,15 @@
col.hor = col.hor, col.vert = col.vert,
log = logpd), dots.without.pch))
if(do.points)
-          {if(ngrid>1){
-              do.call(points, c(list(x = supp, y = psupp1[1:ngrid], pch = pch.u,
+          {if(ngrid0>1){
+              do.call(points, c(list(x = supp1, y = psupp1[1:ngrid0], pch = pch.u,
cex = cex.points, col = col.points), dots.for.points))
-              do.call(points, c(list(x = supp, y = psupp1[2:(ngrid+1)], pch = pch.a,
+              do.call(points, c(list(x = supp1, y = psupp1[2:(ngrid0+1)], pch = pch.a,
cex = cex.points, col = col.points), dots.for.points))
}else{
-              do.call(points, c(list(x = supp, y = 0, pch = pch.u,
+              do.call(points, c(list(x = suppn, y = 0, pch = pch.u,
cex = cex.points, col = col.points), dots.for.points))
-              do.call(points, c(list(x = supp, y = 1, pch = pch.a,
+              do.call(points, c(list(x = suppn, y = 1, pch = pch.a,
cex = cex.points, col = col.points), dots.for.points))
}
}
@@ -467,14 +481,14 @@
col.main = col.inner)

if(do.points)
-          do.call(points, c(list(x = supp,
-                  y = c(0,p(x)(supp[-length(supp)])), pch = pch.u,
+          do.call(points, c(list(x = suppn,
+                  y = c(0,p(x)(suppn[-length(suppn)])), pch = pch.u,
cex = cex.points, col = col.points), dots.for.points))

owarn <- getOption("warn"); options(warn = -1)
-       do.call(plot, c(list(x = stepfun(c(0,p(x)(supp)),
-                            c(NA,supp,NA), right = TRUE),
+       do.call(plot, c(list(x = stepfun(c(0,p(x)(suppn)),
+                            c(NA,suppn,NA), right = TRUE),
main = "", xlim = ylim2, ylim = c(min(supp),max(supp)),
ylab = "q(p)", xlab = "p",
verticals = verticals, do.points = do.points,
@@ -491,22 +505,23 @@
dots.without.pch0 <- dots.without.pch
dots.without.pch0 $col <- NULL - do.call(lines, c(list(x = c(0,p(x)(supp[1])), y = rep(supp[1],2), + do.call(lines, c(list(x = c(0,p(x)(suppn[1])), y = rep(suppn[1],2), col = col.vert), dots.without.pch0)) if(do.points) - {do.call(points, c(list(x = p(x)(supp[-length(supp)]), - y = supp[-1], pch = pch.u, cex = cex.points, + {do.call(points, c(list(x = p(x)(suppn[-length(suppn)]), + y = suppn[-1], pch = pch.u, cex = cex.points, col = col.points), dots.for.points)) - do.call(points, c(list(x = 0, y = supp[1], pch = pch.u, + do.call(points, c(list(x = 0, y = suppn[1], pch = pch.u, cex = cex.points, col = col.points), dots.for.points))} if(verticals && ngrid>1) {dots.without.pch0 <- dots.without.pch dots.without.pch0$col <- NULL

-           do.call(lines, c(list(x = rep(p(x)(supp[1]),2), y = c(supp[1],supp[2]),
-                  col = col.vert), dots.without.pch0))
+           do.call(lines, c(list(x = rep(p(x)(suppn[1]),2),
+                                 y = c(suppn[1],suppn[2]),
+                                 col = col.vert), dots.without.pch0))
}

Modified: pkg/distr/R/plot-methods_LebDec.R
===================================================================
--- pkg/distr/R/plot-methods_LebDec.R	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/R/plot-methods_LebDec.R	2008-07-09 13:56:14 UTC (rev 183)
@@ -186,10 +186,21 @@
supp <- supp[(supp >= xlim[1]) & (supp <= xlim[2])]
}else{grid <- seq(from = lower - 0.1 * dist, to = upper + 0.1 * dist,
length = ngrid)
-           supp <- support(x)
}
+
+     ### possibly thinning out the support
+     thin <- FALSE
+     ngrid <- length(supp)
+     ngrid0 <- getdistrOption("DistrMaxPlotPoints")
+     if (ngrid > ngrid0)
+            {nn <- seq(1,ngrid, by = ngrid %/% ngrid0)
+             thin <- TRUE
+            }else nn <- seq(ngrid)

-     grid <- unique(sort( c(supp, supp-del , grid )))
+     ngrid0 <- length(nn)
+     suppn <- supp[nn]
+
+     grid <- unique(sort( c(suppn, suppn-del , grid )))
pxg <- p(x)(grid)

@@ -209,11 +220,11 @@
}

if(!verticals){
-         grid <- unique(sort( c(supp-del/2, grid )))
-         grid[.isIn(grid,cbind(supp-del/2,supp-del/2))] <- NA
+         grid <- unique(sort( c(suppn-del/2, grid )))
+         grid[.isIn(grid,cbind(suppn-del/2,suppn-del/2))] <- NA
pxg <- p(x)(grid)
}else{
-         xv <- as.vector(t(cbind(supp-del,supp,NA)))
+         xv <- as.vector(t(cbind(suppn-del,suppn,NA)))
pxv <- p(x)(xv)
}

@@ -223,12 +234,12 @@
dots.without.pch))
options(warn = owarn)

-     pxg.d <- p(x)(supp)
-     pxg.d0 <- p(x)(supp-del)
+     pxg.d <- p(x)(suppn)
+     pxg.d0 <- p(x)(suppn-del)
if(do.points){
-        do.call(points, c(list(x = supp, y = pxg.d, pch = pch.a,
+        do.call(points, c(list(x = suppn, y = pxg.d, pch = pch.a,
cex = cex.points, col = col.points), dots.for.points))
-        do.call(points, c(list(x = supp-del, y = pxg.d0, pch = pch.u,
+        do.call(points, c(list(x = suppn-del, y = pxg.d0, pch = pch.u,
cex = cex.points, col = col.points), dots.for.points))
}
if(verticals){
@@ -241,7 +252,7 @@

### quantiles

-     ### fix finite support bounds
+     ### fix finite nort bounds
ixg  <-  grid>=max(q(x)(0),lower) & grid <= min(q(x)(1),upper)
pxg  <-   pxg[ixg]
grid <-  grid[ixg]

Modified: pkg/distr/chm/Distr.chm
===================================================================
(Binary files differ)

Modified: pkg/distr/chm/Distr.toc
===================================================================
--- pkg/distr/chm/Distr.toc	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/chm/Distr.toc	2008-07-09 13:56:14 UTC (rev 183)
@@ -278,6 +278,10 @@
<param name="Local" value="internals.html">
</OBJECT>
<LI> <OBJECT type="text/sitemap">
+<param name="Name" value=".discretizeD">
+<param name="Local" value="internals.html">
+</OBJECT>
+<LI> <OBJECT type="text/sitemap">
<param name="Name" value=".discretizeP">
<param name="Local" value="internals.html">
</OBJECT>

Modified: pkg/distr/chm/distroptions.html
===================================================================
--- pkg/distr/chm/distroptions.html	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/chm/distroptions.html	2008-07-09 13:56:14 UTC (rev 183)
@@ -75,6 +75,10 @@
<dt><code>DistrResolution</code></dt><dd>minimal spacing between two mass points in a discrete distribution, default value: <code>1e-6</code></dd>
<dt><code>DistrCollapse</code></dt><dd>logical; in convolving discrete distributions, shall support points
with distance smaller than <code>DistrResolution</code> be collapsed; default value:  <code>TRUE</code></dd>
+<dt><code>DistrConvPoints</code></dt><dd>numeric; maximal number of support points in convolving discrete distributions;
+if <code>DistrCollapse = TRUE</code> each summand is discretized to this number of support points; default value:  <code>2000</code></dd>
+<dt><code>DistrMaxPlotPoints</code></dt><dd>numeric; maximal number of support points of discrete distributions to be plotted;
+if support is larger, then for plotting it is thinned out; default value:  <code>200</code></dd>
<dt><code>TruncQuantile</code></dt><dd>argument for <code>q</code>-slot at which to truncate; also, for discrete distributions,
support is restricted to [<code>q(TruncQuantile)</code>,<code>q(1-TruncQuantile)</code>], default value: <code>1e-5</code></dd>
<dt><code>DefaultNrFFTGridPointsExponent</code></dt><dd>by default, for e = <code>DefaultNrFFTGridPointsExponent</code>,

Modified: pkg/distr/chm/internals.html
===================================================================
--- pkg/distr/chm/internals.html	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/chm/internals.html	2008-07-09 13:56:14 UTC (rev 183)
@@ -43,6 +43,7 @@
<param name="keyword" value="R:   .plusm">
<param name="keyword" value="R:   .getObjName">
<param name="keyword" value="R:   .discretizeP">
+<param name="keyword" value="R:   .discretizeD">
<param name="keyword" value="R:   .expm.d">
<param name="keyword" value="R:   .expm.c">
<param name="keyword" value="R:   .logm.d">
@@ -95,6 +96,7 @@
.notwithLArg(D)
.getObjName(i = 1)
.discretizeP(D, lower, upper, h)
+.discretizeD(D, n)
.fm(x,f)
.fM(x,f)
.fM2(x,f)
@@ -208,6 +210,9 @@
<tr valign="top"><td><code>D</code></td>
<td>
a distribution object</td></tr>
+<tr valign="top"><td><code>n</code></td>
+<td>
+an integer: number of grid points to discretize to</td></tr>
<tr valign="top"><td><code>i</code></td>
<td>
an integer</td></tr>
@@ -355,8 +360,12 @@
</p>
<p>
<code>.getObjName</code> returns the name of the object in the <code>i</code>th operand.
+</p>
+<p>
<code>.discretizeP</code> discretizes <code>D</code> to a grid of probabilities from
<code>lower</code> to <code>upper</code> with width <code>h</code>.
+<code>.discretizeD</code> discretizes <code>D</code> to a LatticeDistribution with
+n gridpoints.
</p>
<p>
<code>.fm</code>, <code>.fM</code> return the smallest / biggest value in (0,1) such that
@@ -494,6 +503,9 @@
<tr valign="top"><td><code>.discretizeP</code></td>
<td>
<code>numeric</code> &mdash; the probabilities for the grid-values</td></tr>
+<tr valign="top"><td><code>.discretizeD</code></td>
+<td>
+<code>LatticeDistribution</code> &mdash; the discretized distribution</td></tr>
<tr valign="top"><td><code>.makeDd,.makePd, .makeQd</code></td>
<td>
a function with args

Modified: pkg/distr/man/distroptions.Rd
===================================================================
--- pkg/distr/man/distroptions.Rd	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/man/distroptions.Rd	2008-07-09 13:56:14 UTC (rev 183)
@@ -48,6 +48,10 @@
\item{\code{DistrResolution}}{minimal spacing between two mass points in a discrete distribution, default value: \code{1e-6}}
\item{\code{DistrCollapse}}{logical; in convolving discrete distributions, shall support points
with distance smaller than \code{DistrResolution} be collapsed; default value:  \code{TRUE}}
+\item{\code{DistrConvPoints}}{numeric; maximal number of support points in convolving discrete distributions;
+if \code{DistrCollapse = TRUE} each summand is discretized to this number of support points; default value:  \code{700}}
+\item{\code{DistrMaxPlotPoints}}{numeric; maximal number of support points of discrete distributions to be plotted;
+if support is larger, then for plotting it is thinned out; default value:  \code{200}}
\item{\code{TruncQuantile}}{argument for \code{q}-slot at which to truncate; also, for discrete distributions,
support is restricted to [\code{q(TruncQuantile)},\code{q(1-TruncQuantile)}], default value: \code{1e-5}}
\item{\code{DefaultNrFFTGridPointsExponent}}{by default, for e = \code{DefaultNrFFTGridPointsExponent},

Modified: pkg/distr/man/internals.Rd
===================================================================
--- pkg/distr/man/internals.Rd	2008-07-01 18:15:28 UTC (rev 182)
+++ pkg/distr/man/internals.Rd	2008-07-09 13:56:14 UTC (rev 183)
@@ -37,6 +37,7 @@
\alias{.plusm}
\alias{.getObjName}
\alias{.discretizeP}
+\alias{.discretizeD}
\alias{.expm.d}
\alias{.expm.c}
\alias{.logm.d}
@@ -79,6 +80,7 @@
.notwithLArg(D)
.getObjName(i = 1)
.discretizeP(D, lower, upper, h)
+.discretizeD(D, n)
.fm(x,f)
.fM(x,f)
.fM2(x,f)
@@ -142,6 +144,7 @@
\item{Cont}{logical: \code{TRUE} if \code{object} is continuous}
\item{DClass}{character: name of distribution class}
\item{D}{a distribution object}
+  \item{n}{an integer: number of grid points to discretize to}
\item{i}{an integer}
\item{yleft, yright}{extrapolation value beyond left/right endpoint of grid}
\item{h}{numeric: grid width}
@@ -223,8 +226,11 @@
or if its slots \code{p,q} do not have \code{lower.tail} arguments.

\code{.getObjName} returns the name of the object in the \code{i}th operand.
+
\code{.discretizeP} discretizes \code{D} to a grid of probabilities from
\code{lower} to \code{upper} with width \code{h}.
+\code{.discretizeD} discretizes \code{D} to a LatticeDistribution with
+     n gridpoints.

\code{.fm}, \code{.fM} return the smallest / biggest value in (0,1) such that
\code{f}(x) is finite; \code{.fM2} is a variant of \code{.fM} using a
@@ -326,6 +332,7 @@
\code{AbscontDistribution} according to argument \code{DClass}}
\item{.getObjName}{\code{character}}
\item{.discretizeP}{\code{numeric} --- the probabilities for the grid-values}
+\item{.discretizeD}{\code{LatticeDistribution} --- the discretized distribution}
\item{.makeDd,.makePd, .makeQd}{a function with args
\code{x, y, yleft, yright}}
\item{.makeD,.makeDNew}{a function with args \code{x, log = FALSE}}



More information about the Distr-commits mailing list