[Distr-commits] r143 - in pkg: distr/R distr/man distrMod/R

[noreply at r-forge.r-project.org]

Author: ruckdeschel
Date: 2008-05-08 19:31:38 +0200 (Thu, 08 May 2008)
New Revision: 143

Modified:
   pkg/distr/R/internalUtils.R
   pkg/distr/man/internals.Rd
   pkg/distrMod/R/MLEstimator.R
Log:
minor changes in distr:
+new S3 class MLEstimator
+new help functions for prime functions by means of simpson rule (to be used by RtoDPQ -> .makePNew)
+documentation to these functions

Modified: pkg/distr/R/internalUtils.R
===================================================================
--- pkg/distr/R/internalUtils.R	2008-05-08 15:55:55 UTC (rev 142)
+++ pkg/distr/R/internalUtils.R	2008-05-08 17:31:38 UTC (rev 143)
@@ -617,32 +617,87 @@
             return(dfun)
 }
 
+.primefun <- function(f,x, nm = NULL){
+ 
+ h <- diff(x)
+ l <- length(x)
 
+ xm <- (x[-l]+x[-1])/2
 
+ fxm <- f(xm)
+ fx <- f(x)
+ 
+ 
+ fxs  <- 2 * cumsum(fx) - fx - fx[1]
+ fxsm <- 4 * cumsum(fxm)
+
+ fxx <- c(0, (fxs[-1]+fxsm)* h / 6 )
+
+ if (is.null(nm)) nm <- fxx[l]
+
+ fx1 <- approxfun(x, fxx, yright = nm, yleft = 0)
+
+ ffx <- function(u){
+      ffy <- fx1(u) 
+      ffy[u > max(x)] <- nm 
+      ffy[u < min(x)] <- 0
+      return(ffy)
+     }
+
+ return(ffx)
+}
+
+.csimpsum <- function(fx){
+ l <- length(fx)
+ if (l%%2 == 0) {fx <- c(0,fx); l <- l+1}
+ f.even <- fx[seq(l) %% 2 == 0]
+ f.odd  <- fx[seq(l) %% 2 == 1]
+ fs    <- 2 * cumsum(f.odd) - f.odd - f.odd[1]
+ fsm   <- 4 * cumsum(f.even)
+ ff=c(0,(fs[-1]+fsm)/6 )
+ ff
+}
+
 .makePNew <- function(x, dx, h = NULL, notwithLLarg = FALSE,
                       Cont = TRUE, myPf = NULL, pxl = NULL, pxu = NULL){
 
-  p.l <- if(!is.null(pxl)) pxl else cumsum(dx)
+  xm <- min(x); xM <- max(x)
+  xs <- pmin(xM, pmax(xm,x))
+
   if (is.null (h)) h <- 0
 
   if (Cont){
          mfun <- if (is.null (myPf)) approxfun else myPf
-  }else  mfun <- .makePd
+         l <- length(x)
+         xs.u <- xs.l <- xs
+         if(l%%2==0) {
+               xs.l <- c(2*xs[1]-xs[2],xs)
+               xs.u <- c(xs, 2*xs[l]-xs[l-1])
+               l <- l+1
+               }
+         cfun <- .csimpsum
+         xs.l <- xs.l[seq(l)%%2==1]
+         xs.u <- xs.u[seq(l)%%2==1]
+  }else    {
+         mfun <- .makePd
+         cfun <- cumsum
+         xs.u <- xs.l <- xs
+  }       
 
+  p.l <- if(!is.null(pxl)) pxl else cfun(dx)
+
   ## continuity correction by h/2
   nm <- max(p.l)
-  xm <- min(x); xM <- max(x)
-  xs <- pmin(xM,pmax(xm,x))
-  p1.l <- mfun(x = xs, y = p.l, yleft = 0, yright = nm)
+  p1.l <- mfun(x = xs.l, y = p.l, yleft = 0, yright = nm)
   nm <- p1.l(max(x))
 
   if(notwithLLarg){
       ifElsePS <- substitute(if (lower.tail) p1.l(q) else 1 - p1.l(q))
   }else{
-      p.u <- if(!is.null(pxu)) pxu else rev(cumsum(rev(dx)))
+      p.u <- if(!is.null(pxu)) pxu else rev(cfun(rev(dx)))
       ## continuity correction by h/2
       if (!Cont) p.u <- c(p.u[-1],0)
-      p1.u <- mfun(x = xs, y = p.u, yright = 0, yleft = nm)
+      p1.u <- mfun(x = xs.u, y = p.u, yright = 0, yleft = nm)
       rm(p.u)
       ifElsePS <- substitute(if (lower.tail) p1.l(q) else p1.u(q))
   }

Modified: pkg/distr/man/internals.Rd
===================================================================
--- pkg/distr/man/internals.Rd	2008-05-08 15:55:55 UTC (rev 142)
+++ pkg/distr/man/internals.Rd	2008-05-08 17:31:38 UTC (rev 143)
@@ -45,6 +45,8 @@
 \alias{.P2Q}
 \alias{.D2P}
 \alias{.Q2P}
+\alias{.csimpsum}
+\alias{.primefun}
 \alias{coerce,numeric,Integer-method}
 
 \title{Internal functions of package distr}
@@ -104,6 +106,8 @@
                 qL = -Inf, qU = Inf)
 .D2P (d, xx, ql, qu,  ngrid = getdistrOption("DefaultNrGridPoints"))
 .Q2P (q, ngrid = getdistrOption("DefaultNrGridPoints"))
+.csimpsum(fx)
+.primefun(f,x, nm = NULL)
 }
 
 
@@ -188,6 +192,9 @@
  of a regular grid of \code{ngrid} gridpoints to be used in place of \code{xx}.}  
  \item{qL,qU}{argmin / argmax of p()-method}  
  \item{ngrid}{number of gridpoints}
+ \item{fx}{a vector of function evaluations multiplied by the gridwidth}
+ \item{f}{a vector of function evaluations}
+ \item{nm}{an optional right asymptotic value}
 }
 
 \details{
@@ -292,6 +299,14 @@
 from function slots  \code{d} resp. \code{q} 
 by means of function \code{.makePNew} and explicite numeric inversion,  
 respectively.  
+
+\code{.csimpsum} is used internally in \code{.makePNew} to produce
+a prime function out of function evaluations by means of vectorized
+Simpson quadrature method, returning already the function values
+of the prime function on a grid; it is to mimick the behaviour
+of \code{cumsum}. \code{.primefun} is similar but more flexible and
+produces the prime function as a function.
+
 }
 
 
@@ -341,6 +356,8 @@
             lower.tail = TRUE, log.p = FALSE)}} 
 \code{.D2P, .Q2P}{a cdf \code{p} as function \code{function(q, 
             lower.tail = TRUE, log.p = FALSE)}}
+\code{.csimpsum}{a vector of evaluations of the prime function at the grid points}
+\code{.primefun}{the prime function as a function}
 }
 
 \section{Internal classes}{For the ease of method dispatch, there is an internal

Modified: pkg/distrMod/R/MLEstimator.R
===================================================================
--- pkg/distrMod/R/MLEstimator.R	2008-05-08 15:55:55 UTC (rev 142)
+++ pkg/distrMod/R/MLEstimator.R	2008-05-08 17:31:38 UTC (rev 143)
@@ -31,7 +31,7 @@
     res <- MCEstimator(x = x, ParamFamily = ParamFamily, criterion = negLoglikelihood,
                 interval = interval, par = par, ...)
     names(res$criterion) <- "negative log-likelihood"
-    class(res) <- c("MCEstimator", "Estimate")
+    class(res) <- c("MLEstimator", "MCEstimator", "Estimate")
 
     return(res)
 }