[Robast-commits] r344 - in branches/robast-0.7/pkg/ROptEst: R chm inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Aug 20 21:15:57 CEST 2009
Author: ruckdeschel
Date: 2009-08-20 21:15:56 +0200 (Thu, 20 Aug 2009)
New Revision: 344
Modified:
branches/robast-0.7/pkg/ROptEst/R/getInfLM.R
branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R
branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R
branches/robast-0.7/pkg/ROptEst/chm/ROptEst.chm
branches/robast-0.7/pkg/ROptEst/chm/getinfLM.html
branches/robast-0.7/pkg/ROptEst/chm/internals.html
branches/robast-0.7/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
branches/robast-0.7/pkg/ROptEst/man/getinfLM.Rd
branches/robast-0.7/pkg/ROptEst/man/internals.Rd
Log:
fixed several bugs in ROptEst
-getLagrangeMultBy[...]() gain arg a.start to transmit the actual centering / clipping modification in transformed models
+in MSE-method (as iter <= 1 always), clipping bounds in weights in Total Variation case (p=1, k>1) were never updated
-> new helper functino .isVirginW
+in verbose mode now an automatic check for [p]IC is done
+std (standardization in non-standard norms) was not updated correctly in getLagrangeMultBy[...]
- in getinfRobIC_asGRisk/asHampel
+Symmetry slots did not work correctly for non-trivial trafos (at least would have to be reworked); now: all entries are computed
+new helper function .checkPIC for checking within getInfRobIC
+in ..asHampel: bmax had to be converted to matrix in case p=1, k>1
- new example for trafo-based model in script NormalLocationScaleModel.R
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfLM.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfLM.R 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfLM.R 2009-08-20 19:15:56 UTC (rev 344)
@@ -4,7 +4,8 @@
getLagrangeMultByIter <- function(b, L2deriv, risk, trafo,
neighbor, biastype, normtype, Distr,
- z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
+ a.start, z.start, A.start, w.start, std,
+ z.comp, A.comp, maxiter, tol,
verbose, warnit = TRUE){
LMcall <- match.call()
@@ -13,12 +14,13 @@
A <- A.start
w <- w.start
iter <- 0
- a <- 0
+ a <- a.start
## iteration-loop
repeat{
## increment
iter <- iter + 1
+ a.old <- a
z.old <- z
A.old <- A
@@ -28,17 +30,18 @@
cent(w) <- as.numeric(z)
stand(w) <- A
}else if(is(neighbor,"TotalVarNeighborhood")){
- clip(w) <- if(iter==1) c(-b,b)/2 else c(0,b)+a
+ clip(w) <- if(.isVirginW(w)) c(-b,b)/2 else c(0,b)+a
stand(w) <- A
}
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
normW = normtype)
- #print(w)
+ # print(w)
## update centering
z <- getInfCent(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, Distr = Distr, z.comp = z.comp,
w = w, tol.z = .Machine$double.eps^.5)
+ # print(c("z"=z))
if(is(neighbor,"TotalVarNeighborhood")){
a <- z
z <- as.numeric(solve(A,a))
@@ -52,17 +55,23 @@
biastype = biastype, Distr = Distr, A.comp = A.comp,
cent = zc, trafo = trafo, w = w)
+ # print(c("A"=A))
## in case of self-standardization: update norm
normtype.old <- normtype
if(is(normtype,"SelfNorm")){
normtype(risk) <- normtype <- updateNorm(normtype = normtype,
L2 = L2deriv, neighbor = neighbor, biastype = biastype,
- Distr = Distr, V.comp = A.comp, cent = z, stand = A, w = w)
+ Distr = Distr, V.comp = A.comp, cent = zc, stand = A, w = w)
}
## precision and iteration counting
- prec <- max(max(abs(A-A.old)), max(abs(z-z.old)))
- if(verbose)
+ prec <- max(max(abs(A-A.old)), max(abs(a-a.old)),max(abs(z-z.old)))
+# if(verbose)
+# .checkPIC(L2deriv = L2deriv, neighbor = neighbor,
+# Distr = Distr, trafo = trafo, z = zc, A = A, w = w,
+# z.comp = z.comp, A.comp = A.comp)
+
+ if(verbose && iter < maxiter)
cat("current precision in IC algo:\t", prec, "\n")
if(prec < tol) break
if(iter > maxiter){
@@ -75,12 +84,11 @@
## determine LM a
if(is(neighbor,"ContNeighborhood"))
- a <- as.vector(A %*% z)
+ a <- as.vector(A %*% zc)
- std <- if(is(normtype,"QFNorm"))
- QuadForm(normtype) else NULL
+ if(is(normtype,"QFNorm")) std <- QuadForm(normtype)
- return(list(A = A, a = a, z = z, w = w,
+ return(list(A = A, a = a, z = zc, w = w,
biastype = biastype, normtype = normtype,
normtype.old = normtype.old,
risk = risk, std = std,
@@ -90,7 +98,7 @@
getLagrangeMultByOptim <- function(b, L2deriv, risk, FI, trafo,
neighbor, biastype, normtype, Distr,
- z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
+ a.start, z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
verbose, ...){
LMcall <- match.call()
@@ -154,7 +162,7 @@
cent(w0) <- as.numeric(z0)
stand(w0) <- A0
}else if(is(neighbor,"TotalVarNeighborhood")){
- clip(w0) <- if(iter1==1) c(-b,b)/2 else c(0,b)+a0
+ clip(w0) <- if(.isVirginW(w0)) c(-b,b)/2 else c(0,b)+a0
stand(w0) <- A0
}
weight(w0) <- getweight(w0, neighbor = neighbor, biastype = biastype,
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asGRisk.R 2009-08-20 19:15:56 UTC (rev 344)
@@ -187,8 +187,8 @@
## starting values
if(is.null(z.start)) z.start <- numeric(k)
if(is.null(A.start)) A.start <- trafo %*% solve(Finfo)
+ a.start <- as.numeric(A.start %*% z.start)
-
## sort out upper solution if radius = 0
if(identical(all.equal(radius, 0), TRUE))
return(.getUpperSol(L2deriv = L2deriv, b = b, radius = radius,
@@ -198,10 +198,16 @@
QF = std, verbose = verbose, warn = warn))
## determine which entries must be computed
- comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
- z.comp <- comp$"z.comp"
- A.comp <- comp$"A.comp"
+ # by default everything
+ z.comp <- rep(TRUE,k)
+ A.comp <- matrix(rep(TRUE,k*k),nrow=k)
+ # otherwise if trafo == unitMatrix may use symmetry info
+ if(distrMod:::.isUnitMatrix(trafo)){
+ comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
+ z.comp <- comp$"z.comp"
+ A.comp <- comp$"A.comp"
+ }
## selection of the algorithm
pM <- pmatch(tolower(OptOrIter),c("optimize","iterate", "doubleiterate"))
OptOrIter <- pM
@@ -216,7 +222,7 @@
z <- z.start
A <- A.start
b <- 0
- a <- 0
+ a <- a.start
iter <- 0
prec <- 1
iter.In <- 0
@@ -243,7 +249,7 @@
Cov <- 0
Risk <- 1e10
normtype.old <- normtype
- a <- as.numeric(A%*% z)
+ a <- as.numeric(A %*% z)
normtype.opt <- normtype
asGRiskb <- function(b0){
@@ -251,7 +257,7 @@
erg <- getLagrangeMultByOptim(b = b0, L2deriv = L2deriv, risk = risk,
FI = Finfo, trafo = trafo, neighbor = neighbor,
biastype = biastype, normtype = normtype, Distr = Distr,
- z.start = z, A.start = A, w.start = w, std = std,
+ a.start = a, z.start = z, A.start = A, w.start = w, std = std,
z.comp = z.comp, A.comp = A.comp,
maxiter = round(maxiter/50*iter^5), tol = tol^(iter^5/40),
onesetLM = onesetLM, verbose = verbose, ...)
@@ -310,6 +316,7 @@
}else{
repeat{
iter <- iter + 1
+ a.old <- a
z.old <- z
b.old <- b
A.old <- A
@@ -353,7 +360,7 @@
erg <- getLagrangeMultByIter(b = b, L2deriv = L2deriv, risk = risk,
trafo = trafo, neighbor = neighbor, biastype = biastype,
normtype = normtype, Distr = Distr,
- z.start = z, A.start = A, w.start = w,
+ a.start = a, z.start = z, A.start = A, w.start = w,
std = std, z.comp = z.comp,
A.comp = A.comp, maxiter = maxit2, tol = tol,
verbose = verbose, warnit = (OptOrIter!=2))
@@ -376,7 +383,8 @@
## check precision and number of iterations in outer b-loop
prec.old <- prec
- prec <- max(abs(b-b.old), max(abs(A-A.old)), max(abs(z-z.old)))
+ prec <- max(abs(b-b.old), max(abs(A-A.old)),
+ max(abs(z-z.old), max(abs(a-a.old))))
if(verbose)
cat("current precision in IC algo:\t", prec, "\n")
if(prec < tol) break
@@ -415,7 +423,7 @@
biastype = biastype, Distr = Distr,
V.comp = A.comp, cent = a,
stand = A, w = w)
- if(verbose) print(list(Cov=Cov,A=A,c=a,w=w))
+ if(verbose) print(list(Cov=Cov,A=A,a=a,w=w))
if(!is(risk, "asMSE")){
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, clip = b, cent = a, stand = A,
@@ -440,6 +448,11 @@
r = radius,
at = neighbor)))
+ if(verbose)
+ .checkPIC(L2deriv = L2deriv, neighbor = neighbor,
+ Distr = Distr, trafo = trafo, z = z, A = A, w = w,
+ z.comp = z.comp, A.comp = A.comp, ...)
+
return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info, w = w,
biastype = biastype, normtype = normtype,
call = mc, iter = iter, prec = prec, OIcall = OptIterCall,
@@ -530,6 +543,68 @@
return(list(lower=lower, upper=upper))
}
+
+### helper function to check whether (TotalVariation) weight w has already been modified
+.isVirginW <- function(w){
+ w0 <- new("BdStWeight")
+ identical(body(weight(w0)),body(weight(w)))
+}
+
+
+.checkPIC <- function(L2deriv, neighbor, Distr, trafo, z, A, w, z.comp, A.comp, ...){
+ cat("some check:\n-----------\n")
+ nrvalues <- ncol(trafo)
+ pvalues <- nrow(trafo)
+ if(is(neighbor,"ContNeighborhood"))
+ zx <- as.numeric(z)
+ else zx <- numeric(nrvalues)
+
+ L2v.f <- function(x)
+ evalRandVar(L2deriv, as.matrix(x)) [,,1]
+
+ w.f <- function(x) weight(w)(L2v.f(x))
+
+ integrand0 <- function(x,...,ixx){
+ L2v <- as.matrix(L2v.f(x)) - zx
+ wv <- w.f(x)
+ as.numeric(L2v[ixx,]*wv)
+ }
+
+ integrand1 <- function(x,...,ixx){
+ L2v <- as.matrix(L2v.f(x)) - zx
+ AL2v <- A %*% L2v
+ wv <- w.f(x)
+ as.numeric(AL2v[ixx,]*wv)
+ }
+
+ integrand2 <- function(x,...,ixx,jxx){
+ L2v <- as.matrix(L2v.f(x)) - zx
+ AL2v <- integrand1(x,...,ixx = ixx)
+ as.numeric(AL2v*L2v[jxx,])
+ }
+
+ cent0 <- numeric(nrvalues)
+ for(i in 1:nrvalues)
+ if(z.comp[i]) cent0[i] <- E(Distr,integrand0,...,ixx=i)
+
+ cent1 <- numeric(pvalues)
+ for(i in 1:pvalues)
+ cent1[i] <- E(Distr,integrand1,...,ixx=i)
+
+ consist <- 0*trafo
+ for(i in 1:pvalues){
+ for(j in 1:nrvalues){
+ if(A.comp[i,j])
+ consist[i,j] <- E(Distr,integrand2,...,ixx=i,jxx=j)
+ }
+ }
+ consist <- consist-trafo
+ cat("centering (k-space):",cent0,"\n")
+ cat("centering (p-space):",cent1,"\n")
+ cat("Fisher consistency:\n")
+ print(consist)
+}
+
################################################################################
## old routine
################################################################################
Modified: branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/R/getInfRobIC_asHampel.R 2009-08-20 19:15:56 UTC (rev 344)
@@ -167,7 +167,8 @@
## starting values
if(is.null(z.start)) z.start <- numeric(k)
- if(is.null(A.start)) A.start <- trafo
+ if(is.null(A.start)) A.start <- trafo%*%solve(Finfo)
+ a.start <- as.numeric(A.start %*% z.start)
## initialize
if(is(neighbor,"ContNeighborhood")){
@@ -192,10 +193,17 @@
}
## determine which entries must be computed
- comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
- z.comp <- comp$"z.comp"
- A.comp <- comp$"A.comp"
+ # by default everything
+ z.comp <- rep(TRUE,k)
+ A.comp <- matrix(rep(TRUE,k*k),nrow=k)
+ # otherwise if trafo == unitMatrix may use symmetry info
+ if(distrMod:::.isUnitMatrix(trafo)){
+ comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
+ z.comp <- comp$"z.comp"
+ A.comp <- comp$"A.comp"
+ }
+
## selection of the algorithm
pM <- pmatch(tolower(OptOrIter),c("optimize","iterate"))
OptOrIter <- pM
@@ -206,14 +214,16 @@
erg <- getLagrangeMultByOptim(b = b, L2deriv = L2deriv, risk = risk,
FI = Finfo, trafo = trafo, neighbor = neighbor,
biastype = biastype, normtype = normtype, Distr = Distr,
- z.start = z.start, A.start = A.start, w.start = w, std = std,
+ a.start = a.start, z.start = z.start, A.start = A.start,
+ w.start = w, std = std,
z.comp = z.comp, A.comp = A.comp, maxiter = maxiter,
tol = tol, verbose = verbose, ...)
else{
erg <- getLagrangeMultByIter(b = b, L2deriv = L2deriv, risk = risk,
trafo = trafo, neighbor = neighbor, biastype = biastype,
normtype = normtype, Distr = Distr,
- z.start = z.start, A.start = A.start, w.start = w,
+ a.start = a.start, z.start = z.start, A.start = A.start,
+ w.start = w,
std = std, z.comp = z.comp,
A.comp = A.comp, maxiter = maxiter, tol = tol,
verbose = verbose)
@@ -254,10 +264,10 @@
### determine Covariance of pIC
- Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
+ Cov <- as.matrix(getInfV(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, Distr = Distr,
V.comp = A.comp, cent = a,
- stand = A, w = w)
+ stand = A, w = w))
#getAsRisk(risk = asCov(), L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, Distr = Distr, clip = b, cent = a,
@@ -265,7 +275,7 @@
### add some further informations for the pIC-slots info and risk
info <- paste("optimally robust IC for 'asHampel' with bound =", round(b,3))
- trAsCov <- sum(diag(std%*%Cov)); r <- neighbor at radius
+ trAsCov <- sum(diag(std %*% Cov)); r <- neighbor at radius
Risk <- list(trAsCov = list(value = trAsCov,
normtype = normtype),
asCov = Cov,
@@ -276,6 +286,11 @@
r = r,
at = neighbor))
+ if(verbose)
+ .checkPIC(L2deriv = L2deriv, neighbor = neighbor,
+ Distr = Distr, trafo = trafo, z = z, A = A, w = w,
+ z.comp = z.comp, A.comp = A.comp, ...)
+
return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info,
w = w, biastype = biastype, normtype = normtype,
call = mc, iter = iter, prec = prec, OIcall = OptIterCall))
@@ -297,7 +312,7 @@
upper.x <- getUp(Distr)
x <- seq(from = lower.x, to = upper.x, length = nrvalpts)
bmax <- sapply(x,function(x) evalRandVar(ClassIC,x))
- bmax <- sqrt(max(colSums(bmax^2)))
+ bmax <- sqrt(max(colSums(as.matrix(bmax^2))))
cat("numerical approximation of maximal bound:\t", bmax, "\n")
if(b >= bmax){
Modified: branches/robast-0.7/pkg/ROptEst/chm/ROptEst.chm
===================================================================
(Binary files differ)
Modified: branches/robast-0.7/pkg/ROptEst/chm/getinfLM.html
===================================================================
--- branches/robast-0.7/pkg/ROptEst/chm/getinfLM.html 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/chm/getinfLM.html 2009-08-20 19:15:56 UTC (rev 344)
@@ -30,11 +30,11 @@
<pre>
getLagrangeMultByIter(b, L2deriv, risk, trafo,
neighbor, biastype, normtype, Distr,
- z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
+ a.start, z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
verbose, warnit = TRUE)
getLagrangeMultByOptim(b, L2deriv, risk, FI, trafo,
neighbor, biastype, normtype, Distr,
- z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
+ a.start, z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
verbose, ...)
</pre>
@@ -72,9 +72,12 @@
<tr valign="top"><td><code>Distr</code></td>
<td>
object of class <code>"Distribution"</code>. </td></tr>
+<tr valign="top"><td><code>a.start</code></td>
+<td>
+ initial value for the centering constant (in <code>p</code>-space). </td></tr>
<tr valign="top"><td><code>z.start</code></td>
<td>
- initial value for the centering constant. </td></tr>
+ initial value for the centering constant (in <code>k</code>-space). </td></tr>
<tr valign="top"><td><code>A.start</code></td>
<td>
initial value for the standardizing matrix. </td></tr>
@@ -191,7 +194,7 @@
<h3>See Also</h3>
-<p><code><a href="../../RobAStBase/html/InfRobModel-class.html">InfRobModel-class</a></code></p>
+<p><code><a onclick="findlink('RobAStBase', 'InfRobModel-class.html')" style="text-decoration: underline; color: blue; cursor: hand">InfRobModel-class</a></code></p>
<script Language="JScript">
function findlink(pkg, fn) {
Modified: branches/robast-0.7/pkg/ROptEst/chm/internals.html
===================================================================
--- branches/robast-0.7/pkg/ROptEst/chm/internals.html 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/chm/internals.html 2009-08-20 19:15:56 UTC (rev 344)
@@ -6,6 +6,11 @@
<table width="100%"><tr><td>internals_for_ROptEst(ROptEst)</td><td align="right">R Documentation</td></tr></table>
<object type="application/x-oleobject" classid="clsid:1e2a7bd0-dab9-11d0-b93a-00c04fc99f9e">
<param name="keyword" value="R: internals_for_ROptEst">
+<param name="keyword" value="R: .checkUpLow">
+<param name="keyword" value="R: .getUpperSol">
+<param name="keyword" value="R: .getLowerSol">
+<param name="keyword" value="R: .getLowUpB">
+<param name="keyword" value="R: .checkPIC">
<param name="keyword" value=" Internal / Helper functions of package ROptEst">
</object>
@@ -45,6 +50,11 @@
### helper function to return upper & lower bounds for b for b-search
.getLowUpB(L2deriv, Finfo, Distr, normtype, z, A, radius, iter)
+### helper function to check whether (TotalVariation) weight w has already been modified
+.isVirginW(w)
+
+### helper function to check whether (intermediate) results give a pIC
+.checkPIC(L2deriv, neighbor, Distr, trafo, z, A, w, z.comp, A.comp, ...)
</pre>
@@ -123,11 +133,25 @@
<tr valign="top"><td><code>A</code></td>
<td>
standardizing matrix.</td></tr>
+<tr valign="top"><td><code>w</code></td>
+<td>
+a weight of class <code>"BdStWeight"</code></td></tr>
+<tr valign="top"><td><code>z.comp</code></td>
+<td>
+logical vector: indicator which components of <code>z</code> need
+to be computed</td></tr>
+<tr valign="top"><td><code>A.comp</code></td>
+<td>
+logical matrix: indicator which components of <code>A</code> need
+to be computed</td></tr>
<tr valign="top"><td><code>iter</code></td>
<td>
the number of iterations computed so far; used for specifying
a different value of the clipping component of the weight in
total variation case in the very first iteration.</td></tr>
+<tr valign="top"><td><code>...</code></td>
+<td>
+further arguments to be passed on <code>E()</code>.</td></tr>
</table>
@@ -138,11 +162,14 @@
<i>(b_min,b_max)</i>;
<code>.getUpperSol</code> determines the upper case/classical solution and computes
corresponding risks
-<code>.nonNumericb</code> determines the lower case (minimax bias) solution and computes
+<code>.getLowerSol</code> determines the lower case (minimax bias) solution and computes
corresponding risks
<code>.getLowUpB</code> determines a search interval for <code>b</code> to given radius
<code>r</code>, i.e., lower and upper bounds for
<i>(b_min,b_max)</i>
+<code>.isVirginW</code> checks whether the (total variation) weight <code>w</code> in
+the argument has already been modified since creation (<code>TRUE</code> if not)
+<code>.checkPIC</code> checks whether (intermediate) results give a pIC
</p>
@@ -168,6 +195,12 @@
<tr valign="top"><td><code>.checkUpLow</code></td>
<td>
a list with items <code>lower</code> and <code>upper</code> (both numeric).</td></tr>
+<tr valign="top"><td><code>.isVirginW</code></td>
+<td>
+<code>TRUE</code> or <code>FALSE</code></td></tr>
+<tr valign="top"><td><code>.checkPIC</code></td>
+<td>
+nothing is returned; precision values are issued.</td></tr>
</table>
</p>
Modified: branches/robast-0.7/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R
===================================================================
--- branches/robast-0.7/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/inst/scripts/NormalLocationScaleModel.R 2009-08-20 19:15:56 UTC (rev 344)
@@ -92,8 +92,54 @@
## (may take quite some time!)
N0.r.rho1 <- leastFavorableRadius(L2Fam=N0, neighbor=ContNeighborhood(),
risk=asMSE(), rho=0.5)
+###############################################################################
+# a non-trivial trafo:
+###############################################################################
+tfct <- function(x){
+ nms0 <- c("mean","sd")
+ nms <- "comb"
+ fval0 <- x[1]+2*x[2]
+ names(fval0) <- nms
+ mat0 <- matrix(c(1,2), nrow = 1, dimnames = list(nms,nms0))
+ return(list(fval = fval0, mat = mat0))
+}
+
+## corresponding ideal and robust models
+N1.traf <- NormLocationScaleFamily(mean = 0, sd = 1, trafo= tfct)
+N1R.traf <- InfRobModel(center = N1.traf, neighbor = ContNeighborhood(radius = 0.5))
+N2R.traf <- InfRobModel(center = N1.traf, neighbor = TotalVarNeighborhood(radius = 0.5))
+
+### classical solution
+IC.traf.class <- optIC(model=N1.traf,risk=asCov())
+plot(IC.traf.class)
+checkIC(IC.traf.class)
+
+### Hampel solution *=c
+IC.traf.CV.H <- optIC(model = N1R.traf, risk = asHampel(bound = 3),verbose=TRUE)
+plot(IC.traf.CV.H)
+checkIC(IC.traf.CV.H)
+
+### MSE solution *=c
+IC.traf.CV.MSE <- optIC(model = N1R.traf, risk = asMSE(),verbose=TRUE)
+plot(IC.traf.CV.MSE)
+checkIC(IC.traf.CV.MSE)
+
+### Hampel solution *=v
+IC.traf.TV.H <- optIC(model = N2R.traf, risk = asHampel(bound = 6),
+ verbose=TRUE, checkBounds=FALSE)
+plot(IC.traf.TV.H)
+checkIC(IC.traf.TV.H)
+
+### MSE solution *=v
+IC.traf.TV.MSE <- optIC(model = N2R.traf, risk = asMSE(),verbose=TRUE)
+plot(IC.traf.TV.MSE)
+checkIC(IC.traf.TV.MSE)
+
+
+###############################################################################
## one-step estimation
+###############################################################################
## 1. generate a contaminated sample
ind <- rbinom(100, size=1, prob=0.05)
x <- rnorm(100, mean=0, sd=(1-ind) + ind*9)
Modified: branches/robast-0.7/pkg/ROptEst/man/getinfLM.Rd
===================================================================
--- branches/robast-0.7/pkg/ROptEst/man/getinfLM.Rd 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/man/getinfLM.Rd 2009-08-20 19:15:56 UTC (rev 344)
@@ -13,11 +13,11 @@
\usage{
getLagrangeMultByIter(b, L2deriv, risk, trafo,
neighbor, biastype, normtype, Distr,
- z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
+ a.start, z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
verbose, warnit = TRUE)
getLagrangeMultByOptim(b, L2deriv, risk, FI, trafo,
neighbor, biastype, normtype, Distr,
- z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
+ a.start, z.start, A.start, w.start, std, z.comp, A.comp, maxiter, tol,
verbose, ...)
}
@@ -33,7 +33,8 @@
\item{biastype}{object of class \code{"BiasType"} --- the bias type with we work.}
\item{normtype}{object of class \code{"NormType"} --- the norm type with we work.}
\item{Distr}{ object of class \code{"Distribution"}. }
- \item{z.start}{ initial value for the centering constant. }
+ \item{a.start}{ initial value for the centering constant (in \code{p}-space). }
+ \item{z.start}{ initial value for the centering constant (in \code{k}-space). }
\item{A.start}{ initial value for the standardizing matrix. }
\item{w.start}{ initial value for the weight function. }
\item{std}{ matrix of (or which may coerced to) class
Modified: branches/robast-0.7/pkg/ROptEst/man/internals.Rd
===================================================================
--- branches/robast-0.7/pkg/ROptEst/man/internals.Rd 2009-08-20 15:48:23 UTC (rev 343)
+++ branches/robast-0.7/pkg/ROptEst/man/internals.Rd 2009-08-20 19:15:56 UTC (rev 344)
@@ -1,5 +1,10 @@
\name{internals_for_ROptEst}
\alias{internals_for_ROptEst}
+\alias{.checkUpLow}
+\alias{.getUpperSol}
+\alias{.getLowerSol}
+\alias{.getLowUpB}
+\alias{.checkPIC}
\title{Internal / Helper functions of package ROptEst}
@@ -30,6 +35,11 @@
### helper function to return upper & lower bounds for b for b-search
.getLowUpB(L2deriv, Finfo, Distr, normtype, z, A, radius, iter)
+### helper function to check whether (TotalVariation) weight w has already been modified
+.isVirginW(w)
+
+### helper function to check whether (intermediate) results give a pIC
+.checkPIC(L2deriv, neighbor, Distr, trafo, z, A, w, z.comp, A.comp, ...)
}
\arguments{
@@ -59,9 +69,15 @@
\item{radius}{radius of the neighborhood.}
\item{z}{centering constant (in \code{k}-space)}
\item{A}{standardizing matrix.}
+ \item{w}{a weight of class \code{"BdStWeight"}}
+ \item{z.comp}{logical vector: indicator which components of \code{z} need
+ to be computed}
+ \item{A.comp}{logical matrix: indicator which components of \code{A} need
+ to be computed}
\item{iter}{the number of iterations computed so far; used for specifying
a different value of the clipping component of the weight in
total variation case in the very first iteration.}
+ \item{\dots}{further arguments to be passed on \code{E()}.}
}
\details{
@@ -69,11 +85,14 @@
\eqn{(b_{\rm\scriptstyle min},b_{\rm\scriptstyle min})}{(b_min,b_max)};
\code{.getUpperSol} determines the upper case/classical solution and computes
corresponding risks
-\code{.nonNumericb} determines the lower case (minimax bias) solution and computes
+\code{.getLowerSol} determines the lower case (minimax bias) solution and computes
corresponding risks
\code{.getLowUpB} determines a search interval for \code{b} to given radius
\code{r}, i.e., lower and upper bounds for
\eqn{(b_{\rm\scriptstyle min},b_{\rm\scriptstyle min})}{(b_min,b_max)}
+\code{.isVirginW} checks whether the (total variation) weight \code{w} in
+the argument has already been modified since creation (\code{TRUE} if not)
+\code{.checkPIC} checks whether (intermediate) results give a pIC
}
@@ -88,6 +107,8 @@
\item{.getUpperSol}{a return list for \code{getInfRobIC}}
\item{.getLowerSol}{a return list for \code{getInfRobIC}}
\item{.checkUpLow}{a list with items \code{lower} and \code{upper} (both numeric).}
+\item{.isVirginW}{\code{TRUE} or \code{FALSE}}
+\item{.checkPIC}{nothing is returned; precision values are issued.}
}
More information about the Robast-commits
mailing list