[Robast-commits] r135 - in branches/robast-0.6/pkg: ROptEst/R ROptEst/inst/scripts RobAStBase/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jul 31 20:33:23 CEST 2008
Author: stamats
Date: 2008-07-31 20:33:22 +0200 (Thu, 31 Jul 2008)
New Revision: 135
Modified:
branches/robast-0.6/pkg/ROptEst/R/optIC.R
branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R
branches/robast-0.6/pkg/ROptEst/inst/scripts/BinomialModel.R
branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R
branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R
Log:
k-step estimation suboptimal - but seems to work ...
Modified: branches/robast-0.6/pkg/ROptEst/R/optIC.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/R/optIC.R 2008-07-31 16:54:29 UTC (rev 134)
+++ branches/robast-0.6/pkg/ROptEst/R/optIC.R 2008-07-31 18:33:22 UTC (rev 135)
@@ -14,8 +14,13 @@
Finfo = model at center@FisherInfo, trafo = model at center@param at trafo,
upper = upper, maxiter = maxiter, tol = tol, warn = warn,
noLow = noLow)
- options(ow)
+ options(ow)
res$info <- c("optIC", res$info)
+ modIC <- function(L2Fam, IC){
+ infMod <- InfRobModel(L2Fam, model at neighbor)
+ optIC(infMod, risk)
+ }
+ res <- c(res, modifyIC = modIC)
return(generateIC(model at neighbor, model at center, res))
}else{
if(is(model at center@distribution, "UnivariateDistribution")){
@@ -46,8 +51,13 @@
L2derivDistrSymm = L2derivDistrSymm, Finfo = model at center@FisherInfo,
trafo = model at center@param at trafo, z.start = z.start, A.start = A.start,
upper = upper, maxiter = maxiter, tol = tol, warn = warn)
- options(ow)
+ options(ow)
res$info <- c("optIC", res$info)
+ modIC <- function(L2Fam, IC){
+ infMod <- InfRobModel(L2Fam, model at neighbor)
+ optIC(infMod, risk)
+ }
+ res <- c(res, modifyIC = modIC)
return(generateIC(model at neighbor, model at center, res))
}else{
stop("not yet implemented")
@@ -71,11 +81,16 @@
symm = model at center@L2derivDistrSymm[[1]],
Finfo = model at center@FisherInfo, trafo = model at center@param at trafo,
upper = upper, maxiter = maxiter, tol = tol, warn = warn)
- options(ow)
+ options(ow)
if(is(model at neighbor, "ContNeighborhood"))
res$info <- c("optIC", "optIC", res$info, "Optimal IC for 'InfRobModel' with 'ContNeighborhood'!!!")
else
- res$info <- c("optIC", res$info)
+ res$info <- c("optIC", res$info)
+ modIC <- function(L2Fam, IC){
+ infMod <- InfRobModel(L2Fam, model at neighbor)
+ optIC(infMod, risk)
+ }
+ res <- c(res, modifyIC = modIC)
return(generateIC(TotalVarNeighborhood(radius = model at neighbor@radius), model at center, res))
}else{
stop("restricted to 1-dimensional parameteric models")
@@ -98,11 +113,16 @@
neighbor = model at neighbor, risk = risk,
sampleSize = sampleSize, upper = upper, maxiter = maxiter,
tol = tol, warn = warn, Algo = Algo, cont = cont)
- options(ow)
+ options(ow)
if(is(model at neighbor, "ContNeighborhood"))
res$info <- c("optIC", "optIC", res$info, "Optimal IC for 'FixRobModel' with 'ContNeighborhood'!!!")
else
- res$info <- c("optIC", res$info)
+ res$info <- c("optIC", res$info)
+ modIC <- function(L2Fam, IC){
+ infMod <- InfRobModel(L2Fam, model at neighbor)
+ optIC(infMod, risk)
+ }
+ res <- c(res, modifyIC = modIC)
return(generateIC(TotalVarNeighborhood(radius = model at neighbor@radius), model at center, res))
}else{
stop("restricted to 1-dimensional parametric models")
Modified: branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R 2008-07-31 16:54:29 UTC (rev 134)
+++ branches/robast-0.6/pkg/ROptEst/R/radiusMinimaxIC.R 2008-07-31 18:33:22 UTC (rev 135)
@@ -15,10 +15,10 @@
stop("'upRad < loRad' is not fulfilled")
biastype <- biastype(risk)
L2derivDim <- numberOfMaps(L2Fam at L2deriv)
-
+
if(is(normtype(risk),"SelfNorm")||is(normtype(risk),"InfoNorm"))
upRad <- min(upRad,10)
-
+
if(L2derivDim == 1){
ow <- options("warn")
options(warn = -1)
@@ -71,7 +71,7 @@
risk = risk, symm = L2Fam at L2derivSymm[[1]],
Finfo = L2Fam at FisherInfo, upper = upper.b,
trafo = L2Fam at param@trafo, maxiter = maxiter, tol = tol, warn = warn)
- options(ow)
+ options(ow)
res$info <- c("radiusMinimaxIC", paste("radius minimax IC for radius interval [",
round(loRad, 3), ", ", round(upRad, 3), "]", sep=""))
res$info <- rbind(res$info, c("radiusMinimaxIC",
@@ -79,6 +79,11 @@
res$info <- rbind(res$info, c("radiusMinimaxIC",
paste("maximum ", sQuote(class(risk)[1]), "-inefficiency: ",
round(ineff, 3), sep="")))
+ modIC <- function(L2Fam, IC){
+ infMod <- InfRobModel(L2Fam, neighbor)
+ optIC(infMod, risk)
+ }
+ res <- c(res, modifyIC = modIC)
return(generateIC(neighbor, L2Fam, res))
}else{
if(is(L2Fam at distribution, "UnivariateDistribution")){
@@ -109,13 +114,12 @@
p <- nrow(trafo)
FI0 <- trafo%*%solve(Finfo)%*%t(trafo)
FI <- solve(FI0)
-
+
if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
{QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
normtype(risk) <- normtype}
std <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(p)
-
ow <- options("warn")
options(warn = -1)
upper.b <- upper
@@ -125,7 +129,7 @@
if(identical(all.equal(loRad, 0), TRUE)){
loRad <- 0
loRisk <- sum(diag(std%*%FI0))
- loNorm <- normtype
+ loNorm <- normtype
}else{
neighbor at radius <- loRad
resLo <- getInfRobIC(L2deriv = L2deriv, neighbor = neighbor, risk = risk,
@@ -140,7 +144,7 @@
neighbor = neighbor, biastype = biastype,
clip = resLo$b, cent = resLo$a,
stand = resLo$A, trafo = trafo)[[1]]
- loNorm <- resLo$normtype
+ loNorm <- resLo$normtype
}
if(upRad == Inf){
@@ -150,7 +154,7 @@
Distr = L2Fam at distribution,
DistrSymm = L2Fam at distrSymm,
L2derivSymm = L2derivSymm,
- L2derivDistrSymm= L2derivDistrSymm,
+ L2derivDistrSymm= L2derivDistrSymm,
trafo = trafo, z.start = z.start,
A.start = A.start,
maxiter = maxiter, tol = tol,
@@ -170,7 +174,7 @@
normtype(riskUp) <- resUp$normtype
upRisk <- getAsRisk(risk = riskUp, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, clip = resUp$b, cent = resUp$a, stand = resUp$A, trafo = trafo)[[1]]
- upNorm <- resUp$normtype
+ upNorm <- resUp$normtype
}
leastFavR <- uniroot(getIneffDiff, lower = lower, upper = upper,
tol = .Machine$double.eps^0.25, L2Fam = L2Fam, neighbor = neighbor,
@@ -185,7 +189,7 @@
Finfo = L2Fam at FisherInfo, trafo = trafo, z.start = z.start,
A.start = A.start, upper = upper.b, maxiter = maxiter,
tol = tol, warn = warn)
- options(ow)
+ options(ow)
res$info <- c("radiusMinimaxIC", paste("radius minimax IC for radius interval [",
round(loRad, 3), ", ", round(upRad, 3), "]", sep=""))
res$info <- rbind(res$info, c("radiusMinimaxIC",
@@ -193,6 +197,11 @@
res$info <- rbind(res$info, c("radiusMinimaxIC",
paste("maximum ", sQuote(class(risk)[1]), "-inefficiency: ",
round(ineff, 3), sep="")))
+ modIC <- function(L2Fam, IC){
+ infMod <- InfRobModel(L2Fam, neighbor)
+ optIC(infMod, risk)
+ }
+ res <- c(res, modifyIC = modIC)
return(generateIC(neighbor, L2Fam, res))
}else{
stop("not yet implemented")
Modified: branches/robast-0.6/pkg/ROptEst/inst/scripts/BinomialModel.R
===================================================================
--- branches/robast-0.6/pkg/ROptEst/inst/scripts/BinomialModel.R 2008-07-31 16:54:29 UTC (rev 134)
+++ branches/robast-0.6/pkg/ROptEst/inst/scripts/BinomialModel.R 2008-07-31 18:33:22 UTC (rev 135)
@@ -111,22 +111,38 @@
## 2. Kolmogorov(-Smirnov) minimum distance estimator
(est0 <- MDEstimator(x=x, BinomFamily(size=25), interval = c(0, 1)))
-## 3. one-step estimation: radius known
-RobB3 <- InfRobModel(center=BinomFamily(size=25, prob=est0$estimate),
+## 3.1. one-step estimation: radius known
+RobB3 <- InfRobModel(center=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(radius=0.5))
IC9 <- optIC(model=RobB3, risk=asMSE())
-(est1c <- oneStepEstimator(x, IC=IC9, start=est0$estimate))
+(est1c <- oneStepEstimator(x, IC=IC9, start=est0))
-RobB4 <- InfRobModel(center=BinomFamily(size=25, prob=est0$estimate),
+RobB4 <- InfRobModel(center=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(radius=0.25))
IC10 <- optIC(model=RobB4, risk=asMSE())
-(est1v <- oneStepEstimator(x, IC=IC10, start=est0$estimate))
+(est1v <- oneStepEstimator(x, IC=IC10, start=est0))
-## 4. one-step estimation: radius interval
-IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=est0$estimate),
+## 3.2. k-step estimation: radius known
+IC9 <- optIC(model=RobB3, risk=asMSE())
+(est2c <- kStepEstimator(x, IC=IC9, start=est0, steps = 3L))
+
+IC10 <- optIC(model=RobB4, risk=asMSE())
+(est2v <- kStepEstimator(x, IC=IC10, start=est0, steps = 3L))
+
+## 4.1. one-step estimation: radius interval
+IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
-(est2c <- oneStepEstimator(x, IC=IC11, start=est0$estimate))
+(est3c <- oneStepEstimator(x, IC=IC11, start=est0))
-IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=est0$estimate),
+IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
-(est2v <- oneStepEstimator(x, IC=IC12, start=est0$estimate))
+(est3v <- oneStepEstimator(x, IC=IC12, start=est0))
+
+## 4.2. k-step estimation: radius interval
+IC11 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
+ neighbor=ContNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
+(est4c <- kStepEstimator(x, IC=IC11, start=est0, steps = 3L))
+
+IC12 <- radiusMinimaxIC(L2Fam=BinomFamily(size=25, prob=estimate(est0)),
+ neighbor=TotalVarNeighborhood(), risk=asMSE(), loRad=0, upRad=Inf)
+(est4v <- kStepEstimator(x, IC=IC12, start=est0, steps = 3L))
Modified: branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R 2008-07-31 16:54:29 UTC (rev 134)
+++ branches/robast-0.6/pkg/RobAStBase/R/kStepEstimator.R 2008-07-31 18:33:22 UTC (rev 135)
@@ -43,12 +43,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -86,12 +92,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -151,12 +163,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -194,12 +212,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -257,12 +281,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -300,12 +330,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -365,12 +401,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -408,12 +450,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
Modified: branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R
===================================================================
--- branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R 2008-07-31 16:54:29 UTC (rev 134)
+++ branches/robast-0.6/pkg/RobAStBase/R/oneStepEstimator.R 2008-07-31 18:33:22 UTC (rev 135)
@@ -37,12 +37,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -100,12 +106,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -162,12 +174,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
@@ -226,12 +244,18 @@
"computation of IC, asvar and asbias via useLast = FALSE"))
}
if("asCov" %in% names(Risks(IC)))
- asVar <- Risks(IC)$asCov
+ if(length(Risks(IC)$asCov) == 1)
+ asVar <- Risks(IC)$asCov
+ else
+ asVar <- Risks(IC)$asCov$value
else
asVar <- getRiskIC(IC, risk = asCov())$asCov$value
if("asBias" %in% names(Risks(IC))){
- asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ if(length(Risks(IC)$asBias) == 1)
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias
+ else
+ asBias <- neighborRadius(IC)*Risks(IC)$asBias$value
}else{
if(is(IC, "HampIC")){
r <- neighborRadius(IC)
More information about the Robast-commits
mailing list