[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> — the probabilities for the grid-values</td></tr>
+<tr valign="top"><td><code>.discretizeD</code></td>
+<td>
+<code>LatticeDistribution</code> — 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