[Robkalman-commits] r28 - in branches/robKalman_bs: . robKalman robKalman/R robKalman/chm robKalman/demo robKalman/inst robKalman/man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Jun 9 20:08:08 CEST 2009
Author: ruckdeschel
Date: 2009-06-09 20:08:08 +0200 (Tue, 09 Jun 2009)
New Revision: 28
Added:
branches/robKalman_bs/robKalman/
branches/robKalman_bs/robKalman/DESCRIPTION
branches/robKalman_bs/robKalman/NAMESPACE
branches/robKalman_bs/robKalman/R/
branches/robKalman_bs/robKalman/R/0AllClass.R
branches/robKalman_bs/robKalman/R/ACMfilt.R
branches/robKalman_bs/robKalman/R/ACMfilter.R
branches/robKalman_bs/robKalman/R/AllGeneric.R
branches/robKalman_bs/robKalman/R/AllGenerics.R
branches/robKalman_bs/robKalman/R/AllInitialize.R
branches/robKalman_bs/robKalman/R/AllPlot.R
branches/robKalman_bs/robKalman/R/AllShow.R
branches/robKalman_bs/robKalman/R/Psi.R
branches/robKalman_bs/robKalman/R/Util.R
branches/robKalman_bs/robKalman/R/arGM.R
branches/robKalman_bs/robKalman/R/arGMinternal.R
branches/robKalman_bs/robKalman/R/calibrateRLS.R
branches/robKalman_bs/robKalman/R/classKalman.R
branches/robKalman_bs/robKalman/R/mACMfilter.R
branches/robKalman_bs/robKalman/R/mACMinternal.R
branches/robKalman_bs/robKalman/R/rLSfilter.R
branches/robKalman_bs/robKalman/R/recFilter.R
branches/robKalman_bs/robKalman/R/simulateSScont.R
branches/robKalman_bs/robKalman/chm/
branches/robKalman_bs/robKalman/chm/00Index.html
branches/robKalman_bs/robKalman/chm/0robKalman-package.html
branches/robKalman_bs/robKalman/chm/ACMfilt.html
branches/robKalman_bs/robKalman/chm/Logo.jpg
branches/robKalman_bs/robKalman/chm/Rchm.css
branches/robKalman_bs/robKalman/chm/arGM.html
branches/robKalman_bs/robKalman/chm/calibrateRLS.html
branches/robKalman_bs/robKalman/chm/internalACM.html
branches/robKalman_bs/robKalman/chm/internalKalman.html
branches/robKalman_bs/robKalman/chm/internalarGM.html
branches/robKalman_bs/robKalman/chm/internalpsi.html
branches/robKalman_bs/robKalman/chm/internalrLS.html
branches/robKalman_bs/robKalman/chm/recFilter.html
branches/robKalman_bs/robKalman/chm/robKalman.chm
branches/robKalman_bs/robKalman/chm/robKalman.hhp
branches/robKalman_bs/robKalman/chm/robKalman.toc
branches/robKalman_bs/robKalman/chm/simulateSScont.html
branches/robKalman_bs/robKalman/chm/util.html
branches/robKalman_bs/robKalman/demo/
branches/robKalman_bs/robKalman/demo/00Index
branches/robKalman_bs/robKalman/demo/ACMdemo.R
branches/robKalman_bs/robKalman/demo/rLSdemo.R
branches/robKalman_bs/robKalman/inst/
branches/robKalman_bs/robKalman/inst/NEWS
branches/robKalman_bs/robKalman/man/
branches/robKalman_bs/robKalman/man/0robKalman-package.Rd
branches/robKalman_bs/robKalman/man/ACMfilt.Rd
branches/robKalman_bs/robKalman/man/arGM.Rd
branches/robKalman_bs/robKalman/man/calibrateRLS.Rd
branches/robKalman_bs/robKalman/man/internalACM.Rd
branches/robKalman_bs/robKalman/man/internalKalman.Rd
branches/robKalman_bs/robKalman/man/internalarGM.Rd
branches/robKalman_bs/robKalman/man/internalpsi.Rd
branches/robKalman_bs/robKalman/man/internalrLS.Rd
branches/robKalman_bs/robKalman/man/recFilter.Rd
branches/robKalman_bs/robKalman/man/simulateSScont.Rd
branches/robKalman_bs/robKalman/man/util.Rd
Log:
@ Bernhard --- irgendwie hat es nicht gewollt wie gesollt ... sorry!
habe das jetzt gefixt...
Added: branches/robKalman_bs/robKalman/DESCRIPTION
===================================================================
--- branches/robKalman_bs/robKalman/DESCRIPTION (rev 0)
+++ branches/robKalman_bs/robKalman/DESCRIPTION 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,14 @@
+Package: robKalman
+Version: 0.3
+Date: 2009-03-18
+Title: Robust Kalman Filtering
+Description: Routines for Robust Kalman Filtering --- the ACM- and rLS-filter
+Author: Peter Ruckdeschel, Bernhard Spangl, Irina Ursachi, Cezar Chirila
+Maintainer: Peter Ruckdeschel <Peter.Ruckdeschel at itwm.fraunhofer.de>
+Depends: R(>= 2.3.0), methods, graphics, startupmsg, dse1, dse2, MASS, limma, robustbase, numDeriv
+Imports: stats, MASS
+LazyLoad: yes
+License: LGPL-3
+URL: http://robkalman.r-forge.r-project.org/
+LastChangedDate: {$LastChangedDate: 2009-03-31 15:31:30 +0200 (Di, 31 Mrz 2009) $}
+LastChangedRevision: {$LastChangedRevision: 447 $}
Added: branches/robKalman_bs/robKalman/NAMESPACE
===================================================================
--- branches/robKalman_bs/robKalman/NAMESPACE (rev 0)
+++ branches/robKalman_bs/robKalman/NAMESPACE 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,8 @@
+import("methods")
+import("stats")
+import("startupmsg")
+
+export("ACMfilt", "ACMfilter", "arGM", "Euclidnorm",
+ "simulateState", "simulateObs", "rcvmvnorm", "Huberize",
+ "rLScalibrateB", "limitS", "rLSFilter", "KalmanFilter",
+ "recursiveFilter")
Added: branches/robKalman_bs/robKalman/R/0AllClass.R
===================================================================
--- branches/robKalman_bs/robKalman/R/0AllClass.R (rev 0)
+++ branches/robKalman_bs/robKalman/R/0AllClass.R 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,68 @@
+.onLoad <- function(lib, pkg){
+ require("methods", character = TRUE, quietly = TRUE)
+
+}
+
+
+.onAttach <- function(library, pkg)
+{
+buildStartupMessage(pkg="robKalman", library=library, packageHelp=TRUE #,
+ #MANUAL=""
+ )
+ invisible()
+}
+
+#.onUnload <- function(libpath)
+#{
+# library.dynam.unload("distrEx", libpath)
+#}
+#
+#
+## register zoo as "S4"-class
+setOldClass("zoo")
+#
+setClassUnion("Hyperparamtype",
+ c("NULL","matrix", "array", "zoo")
+ )
+
+#
+#
+#
+#
+## positive definite, symmetric matrices with finite entries
+#setClass("ACMcontrol",representation())
+#setClass("rLScontrol",representation())
+#
+#
+#
+#
+# class SSM --- State space model
+setClass("SSM",
+ representation = representation(
+ name = "character", ## name of the ssm
+ F = "Hyperparamtype", ## transition matrix/ces or NULL
+ Z = "Hyperparamtype", ## observation matrix/ces or NULL
+ Q = "Hyperparamtype", ## innovation covariance or NULL
+ V = "Hyperparamtype", ## observation error covariance or NULL
+ p = "numeric", ## state dimension
+ q = "numeric"), ## observation dimension
+ prototype = prototype(name = gettext("a state space"),
+ F = NULL,
+ Z = NULL,
+ Q = NULL,
+ V = NULL,
+ p = 1,
+ q = 1),
+ contains = "VIRTUAL")
+
+# class TimeInvariantSSM
+setClass("TimeInvariantSSM",
+ prototype = prototype(name = gettext("a time-invariant state space"),
+ F = 1,
+ Z = 1,
+ Q = 1,
+ V = 1,
+ p = 1,
+ q = 1),
+ contains = "SSM")
+
\ No newline at end of file
Added: branches/robKalman_bs/robKalman/R/ACMfilt.R
===================================================================
--- branches/robKalman_bs/robKalman/R/ACMfilt.R (rev 0)
+++ branches/robKalman_bs/robKalman/R/ACMfilt.R 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,78 @@
+ACMfilt <- function (x, gm, s0=0,
+ psi="Hampel", a=2.5, b=a, c=5.0,
+ flag="weights", lagsmo=TRUE)
+{
+###########################################
+##
+## R-function: ACMfilt - approximate conditional-mean filtering (wrapper)
+## author: Bernhard Spangl
+## version: 1.1 (2007-08-13 and 2006-08-31)
+## References:
+## [Mart79c] R.D. Martin, Approximate Conditional-mean Type Smoothers
+## and Interpolators (1979)
+## [Mart81b] R.D. Martin, Robust Methods for Time Series (1981)
+## [MarT82b] R.D. Martin & D.J. Thomson, Robust-resistent Spectrum
+## Estimation (1982)
+##
+###########################################
+
+## Paramters:
+## x ... univariate time series (vector)
+## gm ... list as produced by function 'arGM' which includes components
+## 'ar' containing the AR(p) coefficient estimates, 'sinnov' containing
+## innovation scale estiamtes from AR(p) fits of orders 1 through p;
+## 'Cx' containing an estimate of the p by p autocovariance matrix,
+## and 'mu', the estimated mean of 'x'.
+## s0 ... scale of nominal Gaussian component of the additive noise
+## psi ... influence function to be used (default: "Hampel",
+## only Hampel's psi function available at the moment)
+## a, b, c ... tuning constants for Hampel's psi-function
+## (defaul: a=b=2.5, c=5.0)
+## flag ... character, if "weights" (default), use psi(t)/t to calculate
+## the weights; if "deriv", use psi'(t)
+## lagsmo ... logical, if TRUE (default) lag p-1 smoothing is performed;
+## if FALSE filtering from the top of ^X_t is performed
+
+## Variable definitions:
+
+ N <- length(x)
+ phi <- gm$ar
+ p <- length(phi)
+ si <- gm$sinnov[p]
+ Cx <- gm$Cx
+ Phi <- cbind(rbind(phi[-p], diag(rep(1, (p-1)))), c(phi[p], rep(0, (p-1))))
+ Q <- matrix(0, p, p)
+## Q <- diag(rep(0, p))
+ Q[1, 1] <- si^2
+
+ m0 <- rep(0, p)
+ H <- matrix(c(1, rep(0, (p-1))), 1, p)
+ V <- matrix(s0^2)
+ psi <- .psi(psi)
+
+ ## Centering:
+ x <- x - gm$mu
+ ACMres <- ACMfilter(Y=matrix(x,1,N), a=m0, S=Cx, F=Phi, Q=Q, Z=H, V=V, s0=s0, psi=psi, apsi=a, bpsi=b, cpsi=c, flag=flag)
+
+ X.ck <- ACMres$Xf; X.ck <- X.ck[,2:(N+1)]
+ X <- ACMres$Xrf; X <- X[,2:(N+1)]
+ st <- as.numeric(unlist(ACMres$rob1L))
+
+
+ if (!lagsmo) {
+ x.ck <- X.ck[1, ]
+ x <- X[1, ]
+ } else {
+ x.ck <- c(X.ck[p, p:N], X.ck[(p-1):1, N])
+ x <- c(X[p, p:N], X[(p-1):1, N])
+ }
+
+## ARmodel <- .ARmodel(x, p)
+## y <- ARmodel$y
+## Z <- ARmodel$Z
+## r <- resid(lm.fit(Z, y))
+
+ return(list(filt.ck=x.ck +gm$mu, filt=x + gm$mu, st=st)) #,
+## r=c(rep(NA, p), r)))
+
+}
Added: branches/robKalman_bs/robKalman/R/ACMfilter.R
===================================================================
--- branches/robKalman_bs/robKalman/R/ACMfilter.R (rev 0)
+++ branches/robKalman_bs/robKalman/R/ACMfilter.R 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,92 @@
+.getcorrCovACM <- function (S1, K, Z, W=diag(nrow(Z)))
+{
+###########################################
+##
+## R-function: .corrCov - computes filtering error covarince matrix
+## (internal function)
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-22)
+##
+###########################################
+
+## Paramters:
+## S1 ... prediction error covariance matrix
+## K ... Kalman gain
+## W ... weight matrix
+## Z ... observation matrix
+
+ S1 - K %*% W %*% Z %*% S1
+}
+
+##steps for classical Kalman filter (cK)
+.ACMinitstep <- function(a, S, ...)
+ {dots <- list(...)
+ if(hasArg("s0"))
+ s0<-dots$"s0"
+ else
+ s0<-NULL
+ list( x0 = a, S0 = S, s0 = s0)}
+
+
+.ACMpredstep <- function (x0, S0, F, Q, rob0, s0, ...) ### S=P F= Phi
+{
+###########################################
+##
+## R-function: .ACMpredstep - prediction step (internal function)
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-22)
+##
+###########################################
+
+## Paramters:
+## x0 ... state vector (filter estimate)
+## F=Phi ... design matrix of state equation
+## S0 ... filtering error covariance matrix
+## Q ... covariance matrix of state innovation process
+## rob0 ... general robust parameter --- here: scale s0 of nominal Gaussain component of additive noise
+ S1 <- .getpredCov(S0, F, Q)
+ return(list(x1 = F %*% x0, S1 = S1, rob1 = sqrt(S1[1, 1] + s0), Ind=1))
+}
+
+.ACMcorrstep <- function (y, x1, S1, Z, V, rob1, dum=NULL, psi, apsi, bpsi, cpsi, flag, ...)
+{
+###########################################
+##
+## R-function: .ACMcorrstep - correction step (internal function)
+## author: Bernhard Spangl
+## version: 1.0 (2006-05-22)
+##
+###########################################
+
+## Paramters:
+## y ... univariate time series
+## x1 ... state vector (one-step-ahead predictor)
+## rob1 ... general robust parameter --- here st ... time-dependent scale parameter
+## S1 ... prediction error covariance matrix
+## Z ... observation matrix
+## dum ... dummy variable for compatibility with ... argument of calling function
+## V ... covariance matrix of observation noise
+## psi ... influence function to be used
+## a, b, c ... tuning constants for Hampel's psi-function
+## (defaul: a=b=2.5, c=5.0)
+## flag ... character, if "weights" (default), use psi(t)/t to calculate
+## the weights; if "deriv", use psi'(t)
+
+ st <- rob1
+
+ K <- .getKG(S1, Z, V)
+
+ rst <- (y - x1[1])/st
+
+ ps <- psi(rst, apsi, bpsi, cpsi)
+ dx <- K * st * ps
+ x0 <- x1 + dx
+
+ ind <- (abs(rst-ps)>10^-8)
+
+ w <- psi(rst, apsi, bpsi, cpsi, flag)
+
+ S0 <- .getcorrCovACM(S1, K, Z, W = w*diag(rep(1, nrow(Z))))
+
+ return(list(x0 = x0, K = K, S0 = S0, Ind=ind, rob0=rob1))
+}
Added: branches/robKalman_bs/robKalman/R/AllGeneric.R
===================================================================
--- branches/robKalman_bs/robKalman/R/AllGeneric.R (rev 0)
+++ branches/robKalman_bs/robKalman/R/AllGeneric.R 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,10 @@
+#if(!isGeneric("type")){
+# setGeneric("type", function(object) standardGeneric("type"))
+#}
+#if(!isGeneric("center")){
+# setGeneric("center", function(object) standardGeneric("center"))
+#}
+#if(!isGeneric("center<-")){
+# setGeneric("center<-", function(object, value) standardGeneric("center<-"))
+#}
+#
Added: branches/robKalman_bs/robKalman/R/AllGenerics.R
===================================================================
--- branches/robKalman_bs/robKalman/R/AllGenerics.R (rev 0)
+++ branches/robKalman_bs/robKalman/R/AllGenerics.R 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,61 @@
+############################################################################
+# Access methods
+############################################################################
+
+if(!isGeneric("name"))
+ setGeneric("name", function(object) standardGeneric("name"))
+
+if(!isGeneric("getp"))
+ setGeneric("getp", function(object) standardGeneric("getp"))
+if(!isGeneric("getq"))
+ setGeneric("getq", function(object) standardGeneric("getq"))
+
+if(!isGeneric("getF"))
+ setGeneric("getF", function(object,t) standardGeneric("getF"))
+if(!isGeneric("getZ"))
+ setGeneric("getZ", function(object,t) standardGeneric("getZ"))
+if(!isGeneric("getQ"))
+ setGeneric("getQ", function(object,t) standardGeneric("getQ"))
+if(!isGeneric("getV"))
+ setGeneric("getV", function(object,t) standardGeneric("getV"))
+
+
+############################################################################
+# Replacement methods
+############################################################################
+
+if(!isGeneric("name<-"))
+ setGeneric("name<-",
+ function(object, value) standardGeneric("name<-"))
+
+############################################################################
+# generics to "usual" methods
+############################################################################
+
+
+# general methods
+
+if(!isGeneric("isOldVersion"))
+ setGeneric("isOldVersion", function(object) standardGeneric("isOldVersion"))
+
+if(!isGeneric("conv2NewVersion"))
+ setGeneric("conv2NewVersion",
+ function(object) standardGeneric("conv2NewVersion"))
+### setting gaps
+
+if(!isGeneric("setgaps"))
+ setGeneric("setgaps", function(object, ...) standardGeneric("setgaps"))
+
+#### generics for log, log10, lgamma, gamma
+
+
+if(!isGeneric("log"))
+ setGeneric("log") #, function(x, base) standardGeneric("log"))
+if(!isGeneric("log10"))
+ setGeneric("log10")
+if(!isGeneric("lgamma"))
+ setGeneric("lgamma")
+if(!isGeneric("gamma"))
+ setGeneric("gamma")
+
+
Added: branches/robKalman_bs/robKalman/R/AllInitialize.R
===================================================================
--- branches/robKalman_bs/robKalman/R/AllInitialize.R (rev 0)
+++ branches/robKalman_bs/robKalman/R/AllInitialize.R 2009-06-09 18:08:08 UTC (rev 28)
@@ -0,0 +1,1012 @@
+#### as to whether to use Generating functions or to use initialize methods:
+#### http://tolstoy.newcastle.edu.au/R/e2/devel/07/01/1976.html
+
+################################################################################
+## SPACES
+################################################################################
+
+setMethod("initialize", "Reals",
+ function(.Object) {
+ .Object at dimension <- 1
+ .Object at name <- gettext("Real Space")
+ .Object
+ })
+
+
+setMethod("initialize", "Naturals",
+ function(.Object) {
+ .Object at dimension <- 1
+ .Object at name <- gettext("Grid of Naturals")
+ .Object
+ })
+
+
+################################################################################
+## PARAMETERS
+################################################################################
+
+setMethod("initialize", "GeomParameter",
+ function(.Object, prob = .5) {
+ .Deprecated(new = "new(\"NbinomParameter\"(size = 1, prob, name)",
+ package = "distr",
+ msg = gettext(
+"Class 'GeomParameter' is no longer needed and will be replaced by \nclass 'NbinomParameter' soon."
+ ))
+ .Object at prob <- prob
+ .Object at name <- gettext("Parameter of a Geometric distribution")
+ .Object
+ })
+################################################################################
+## DISTRIBUTIONS
+################################################################################
+
+## Class: UnivariateDistribution
+###produces difficulties in coercing...:
+#
+#setMethod("initialize", "UnivariateDistribution",
+# function(.Object, r = NULL, d = NULL, p = NULL, q = NULL,
+# param = NULL, img = new("Reals"),
+# .withSim = FALSE, .withArith = FALSE) {
+# if(is.null(r)) {
+# stop("You have at least to give the slot r.")
+# return(invisible())}
+# ### Attention: no checking!!!
+# .Object at img <- img
+# .Object at param <- param
+# .Object at d <- d
+# .Object at p <- p
+# .Object at q <- q
+# .Object at r <- r
+# .Object at .withSim <- .withSim
+# .Object at .withArith <- .withArith
+# .Object })
+
+## class AbscontDistribution
+setMethod("initialize", "AbscontDistribution",
+ function(.Object, r = NULL, d = NULL, p = NULL, q = NULL,
+ gaps = NULL, param = NULL, img = new("Reals"),
+ .withSim = FALSE, .withArith = FALSE) {
+
+ ## don't use this if the call is new("AbscontDistribution")
+ LL <- length(sys.calls())
+ if(sys.calls()[[LL-3]] == "new(\"AbscontDistribution\")")
+ {return(.Object)}
+
+ if(is.null(r))
+ warning("you have to specify slot r at least")
+
+ ## TOBEDONE Errorkanal
+
+ dpq.approx <- 0
+
+ dfun <- d
+ pfun <- p
+ qfun <- q
+
+ if(is.null(d)) {
+ .withSim <- TRUE
+ dpq <- RtoDPQ(r)
+ dpq.approx <- 1
+ dfun <- dpq$dfun}
+
+ if(is.null(p)) {
+ .withSim <- TRUE
+ if(dpq.approx == 0) {dpq <- RtoDPQ(r)}
+ dpq.approx <- 1
+ pfun <- dpq$pfun}
+
+ if(is.null(q)) {
+ .withSim <- TRUE
+ if(dpq.approx == 0) {dpq <- RtoDPQ(r)}
+ qfun <- dpq$qfun}
+
+ .Object at img <- img
+ .Object at param <- param
+ .Object at d <- dfun
+ .Object at p <- pfun
+ .Object at q <- qfun
+ .Object at r <- r
+ .Object at gaps <- gaps
+ .Object at .withSim <- .withSim
+ .Object at .withArith <- .withArith
+ .Object })
+
+## class AffLinAbscontDistribution
+setMethod("initialize", "AffLinAbscontDistribution",
+ function(.Object, r = NULL, d = NULL, p = NULL, q = NULL, gaps = NULL,
+ a = 1, b = 0, X0 = Norm(), param = NULL, img = new("Reals"),
+ .withSim = FALSE, .withArith = FALSE) {
+ X <- new("AbscontDistribution", r = r, d = d, p = p, q = q, gaps = gaps,
+ param = param, img = img, .withSim = .withSim,
+ .withArith = .withArith)
+ .Object at gaps <- X at gaps
+ .Object at img <- X at img
+ .Object at param <- X at param
+ .Object at a <- a
+ .Object at b <- b
+ .Object at X0 <- X0
+ .Object at d <- X at d
+ .Object at p <- X at p
+ .Object at q <- X at q
+ .Object at r <- X at r
+ .Object at .withSim <- .withSim
+ .Object at .withArith <- .withArith
+ .Object
+})
+
+## Class: DiscreteDistribution
+setMethod("initialize", "DiscreteDistribution",
+ function(.Object, r = NULL, d = NULL, p = NULL, q = NULL,
+ support = NULL, param = NULL, img = new("Reals"),
+ .withSim = FALSE, .withArith = FALSE) {
+
+ ## don't use this if the call is new("DiscreteDistribution")
+ LL <- length(sys.calls())
+ if(sys.calls()[[LL-3]] == "new(\"DiscreteDistribution\")")
+ {return(.Object)}
+
+ if(is.null(r))
+ warning("you have to specify slot r at least")
+
+ if(is.null(support))
+ .Object at support <- as.numeric(names(table(r(10^6))))
+
+ else .Object at support <- support
+
+ len = length(support)
+
+ if(len > 1){
+ if(min(support[2:len] - support[1:(len - 1)]) <
+ getdistrOption("DistrResolution"))
+ stop("grid too narrow --> change DistrResolution")
+ }
+
+ dpq.approx <- 0
+
+ dfun <- d
+ pfun <- p
+ qfun <- q
+
+ if(is.null(d)) {
+ .withSim <- TRUE
+ dpq <- RtoDPQ.d(r)
+ dpq.approx <- 1
+ dfun <- dpq$dfun
+ }
+
+ if(is.null(p)) {
+ .withSim <- TRUE
+ if(dpq.approx==0) dpq <- RtoDPQ.d(r)
+ dpq.approx <- 1
+ pfun <- dpq$pfun
+ }
+
+ if(is.null(q)) {
+ .withSim <- TRUE
+ if(dpq.approx==0) dpq <- RtoDPQ.d(r)
+ qfun <- dpq$qfun
+ }
+
+ .Object at img <- img
+ .Object at param <- param
+ .Object at d <- dfun
+ .Object at p <- pfun
+ .Object at q <- qfun
+ .Object at r <- r
+ .Object at .withSim <- .withSim
+ .Object at .withArith <- .withArith
+ .Object
+ })
+
+## Class: AffLinDiscreteDistribution
+setMethod("initialize", "AffLinDiscreteDistribution",
+ function(.Object, r = NULL, d = NULL, p = NULL, q = NULL,
+ support = NULL, a = 1, b = 0, X0 = Binom(), param = NULL,
+ img = new("Reals"), .withSim = FALSE, .withArith = FALSE) {
+ ## don't use this if the call is new("DiscreteDistribution")
+ LL <- length(sys.calls())
+ if(sys.calls()[[LL-3]] == "new(\"AffLinDiscreteDistribution\")")
+ X <- new("DiscreteDistribution")
+ else X <- new("DiscreteDistribution", r = r, d = d, p = p, q = q, support = support,
+ param = param, img = img, .withSim = .withSim,
+ .withArith = .withArith)
+ .Object at support <- X at support
+ .Object at img <- X at img
+ .Object at param <- X at param
+ .Object at a <- a
+ .Object at b <- b
+ .Object at X0 <- X0
+ .Object at d <- X at d
+ .Object at p <- X at p
+ .Object at q <- X at q
+ .Object at r <- X at r
+ .Object at .withSim <- .withSim
+ .Object at .withArith <- .withArith
+ .Object
+})
+
+## Class: LatticeDistribution
+setMethod("initialize", "LatticeDistribution",
+ function(.Object, r = NULL, d = NULL, p = NULL, q = NULL,
+ support = NULL, lattice = NULL, param = NULL,
+ img = new("Reals"), .withSim = FALSE, .withArith = FALSE) {
+
+
+ LL <- length(sys.calls())
+ if(sys.calls()[[LL-3]] == "new(\"LatticeDistribution\")")
+ D <- new("DiscreteDistribution")
+ else
+ D <- new("DiscreteDistribution", r = r, d = d, p = p,
+ q = q, support = support, param = param, img = img,
+ .withSim = .withSim, .withArith = .withArith)
+
+
+ OS <- D at support
+
+ if(is.null(lattice))
+ { if(! .is.vector.lattice(OS))
+ stop("Support as given/generated is not a lattice.")
+ .Object at lattice <- .make.lattice.es.vector(OS)
+ }else{
+ .Object at lattice <- lattice
+ }
+
+
+ .Object at support <- OS
+ .Object at img <- D at img
+ .Object at param <- D at param
+ .Object at d <- D at d
+ .Object at p <- D at p
+ .Object at q <- D at q
+ .Object at r <- D at r
+ .Object at .withSim <- .withSim
+ .Object at .withArith <- .withArith
+ .Object
+ })
+
+## Class: AffLinLatticeDistribution
+setMethod("initialize", "AffLinLatticeDistribution",
+ function(.Object, r = NULL, d = NULL, p = NULL, q = NULL,
+ support = NULL, lattice = NULL, a = 1, b = 0, X0 = Binom(),
+ param = NULL, img = new("Reals"), .withSim = FALSE,
+ .withArith = FALSE) {
+
+ LL <- length(sys.calls())
+ if(sys.calls()[[LL-3]] == "new(\"AffLinLatticeDistribution\")")
+ X <- new("LatticeDistribution")
+ else X <- new("LatticeDistribution", r = r, d = d, p = p, q = q,
+ support = support, lattice = lattice, param = param,
+ img = img, .withSim = .withSim,
+ .withArith = .withArith)
+
+ .Object at support <- X at support
+ .Object at lattice <- X at lattice
+ .Object at img <- X at img
+ .Object at param <- X at param
+ .Object at a <- a
+ .Object at b <- b
+ .Object at X0 <- X0
+ .Object at d <- X at d
+ .Object at p <- X at p
+ .Object at q <- X at q
+ .Object at r <- X at r
+ .Object at .withSim <- .withSim
+ .Object at .withArith <- .withArith
+ .Object
+})
+
+######### particular discrete distributions
+
+### Class: Dirac distribution
+setMethod("initialize", "Dirac",
+ function(.Object, location = 0, .withArith = FALSE) {
+ .Object at img <- new("Reals")
+ .Object at param <- new("DiracParameter", location = location)
+ .Object at r <- function(n){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute({ rep(locationSub, n)},
+ list(locationSub = location)
+ )
+ body(.Object at d) <- substitute(
+ { y <- rep(locationSub, length(x))
+ d0 <- mapply(function(x,y)
+ as.numeric(isTRUE(all.equal(x,y))),
+ x = x, y = y)
+ if (log) d0 <- log(d0)
+ return(d0)
+ }, list(locationSub = location)
+ )
+ body(.Object at p) <- substitute(
+ {p0 <-as.numeric(q + .Machine$double.eps^.5 >=
+ locationSub)
+ if (!lower.tail) p0 <- 1-p0
+ if (log.p) p0 <- log(p0)
+ return(p0)
+ },
+ list(locationSub = location)
+ )
+ body(.Object at q) <- substitute(
+ { if (log.p) p <- exp(p)
+ if(any((p < 0)|(p > 1)))
+ warning("q Method of class Dirac produced NaN's.")
+ q0 <- ifelse((p < 0)|(p > 1), NaN, locationSub)
+ return(q0)
+ },
+ list(locationSub = location)
+ )
+ .Object at support <- location
+ .Object at lattice <- new("Lattice", pivot = location, width = 1,
+ Length = 1)
+ .Object at .withArith <- .withArith
+ .Object
+ })
+
+## Class: binomial distribution
+setMethod("initialize", "Binom",
+ function(.Object, size = 1, prob = 0.5, .withArith = FALSE) {
+ .Object at img <- new("Naturals")
+ .Object at param <- new("BinomParameter", size = size, prob = prob)
+ .Object at support <- 0:size
+ .Object at r <- function(n){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute(
+ { rbinom(n, size = sizeSub, prob = probSub) },
+ list(sizeSub = size, probSub = prob)
+ )
+ body(.Object at d) <- substitute(
+ { dbinom(x, size = sizeSub, prob = probSub,
+ log = log) },
+ list(sizeSub = size, probSub = prob)
+ )
+ body(.Object at p) <- substitute(
+ { pbinom(q, size = sizeSub, prob = probSub,
+ lower.tail = lower.tail, log.p = log.p) },
+ list(sizeSub = size, probSub = prob)
+ )
+ body(.Object at q) <- substitute(
+ { qbinom(p, size = sizeSub, prob = probSub,
+ lower.tail = lower.tail, log.p = log.p) },
+ list(sizeSub = size, probSub = prob)
+ )
+ .Object at support = 0:size
+ .Object at lattice = new("Lattice", pivot = 0, width = 1,
+ Length = size+1)
+ .Object at .withArith <- .withArith
+ .Object
+ })
+
+## Class: hypergeometric distribution
+setMethod("initialize", "Hyper",
+ function(.Object, m = 1, n = 1, k = 1, .withArith = FALSE) {
+ .Object at img <- new("Naturals")
+ .Object at param <- new("HyperParameter", m = m, n = n, k = k)
+ .Object at support <- 0:k
+ .Object at r <- function(nn){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute(
+ { rhyper(nn, m = mSub, n = nSub, k = kSub) },
+ list(mSub = m, nSub = n, kSub = k)
+ )
+ body(.Object at d) <- substitute(
+ { dhyper(x, m = mSub, n = nSub, k = kSub,
+ log = log) },
+ list(mSub = m, nSub = n, kSub = k)
+ )
+ body(.Object at p) <- substitute(
+ { phyper(q, m = mSub, n = nSub, k = kSub,
+ lower.tail = lower.tail, log.p = log.p)
+ },
+ list(mSub = m, nSub = n, kSub = k)
+ )
+ body(.Object at q) <- substitute(
+ { qhyper(p, m = mSub, n = nSub, k = kSub,
+ lower.tail = lower.tail, log.p = log.p)
+ },
+ list(mSub = m, nSub = n, kSub = k)
+ )
+ .Object at support <- seq(from = 0, to = min(k,m), by = 1)
+ .Object at lattice <- new("Lattice", pivot = 0, width = 1,
+ Length = min(k,m)+1 )
+ .Object at .withArith <- .withArith
+ .Object
+ })
+
+## Class: Poisson distribution
+setMethod("initialize", "Pois",
+ function(.Object, lambda = 1, .withArith=FALSE) {
+ .Object at img <- new("Naturals")
+ .Object at param <- new("PoisParameter", lambda = lambda)
+ .Object at r <- function(n){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute({ rpois(n, lambda = lambdaSub) },
+ list(lambdaSub = lambda)
+ )
+ body(.Object at d) <- substitute({ dpois(x, lambda = lambdaSub,
+ log = log) },
+ list(lambdaSub = lambda)
+ )
+ body(.Object at p) <- substitute({ ppois(q, lambda = lambdaSub,
+ lower.tail = lower.tail,
+ log.p = log.p) },
+ list(lambdaSub = lambda)
+ )
+ body(.Object at q) <- substitute({ qpois(p, lambda = lambdaSub,
+ lower.tail = lower.tail,
+ log.p = log.p) },
+ list(lambdaSub = lambda)
+ )
+ .Object at support <- seq(from = 0, by = 1, to =
+ qpois(getdistrOption("TruncQuantile"),
+ lambda = lambda, lower.tail = FALSE)
+ + 2
+ )
+ .Object at lattice <- new("Lattice", pivot = 0, width = 1,
+ Length = Inf)
+ .Object at .withArith <- .withArith
+ .Object
+ })
+
+## Class: negative binomial distribution
+setMethod("initialize", "Nbinom",
+ function(.Object, size = 1, prob = 0.5, .withArith = FALSE) {
+ .Object at img <- new("Naturals")
+ .Object at param <- new("NbinomParameter", size = size, prob = prob)
+ .Object at r <- function(n){}
+ .Object at d <- function(x, log = FALSE){}
+ .Object at p <- function(q, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
+ body(.Object at r) <- substitute(
+ { rnbinom(n, size = sizeSub, prob = probSub) },
+ list(sizeSub = size, probSub = prob)
+ )
+ body(.Object at d) <- substitute(
+ { dnbinom(x, size = sizeSub, prob = probSub,
+ log = log) },
+ list(sizeSub = size, probSub = prob)
+ )
+ body(.Object at p) <- substitute(
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robkalman -r 28
More information about the Robkalman-commits
mailing list