Nothing
###############################################################################
## get optimally robust IC for convex asymptotic risks
###############################################################################
setMethod("getInfRobIC", signature(L2deriv = "UnivariateDistribution",
risk = "asGRisk",
neighbor = "UncondNeighborhood"),
function(L2deriv, risk, neighbor, symm, Finfo, trafo, upper = NULL,
lower = NULL, maxiter, tol,
warn, noLow = FALSE, verbose = NULL, ...){
if(missing(verbose)|| is.null(verbose))
verbose <- getRobAStBaseOption("all.verbose")
if(missing(warn)|| is.null(warn)) warn <- FALSE
biastype <- biastype(risk)
normtype <- normtype(risk)
radius <- neighbor@radius
p <- nrow(trafo)
k <- ncol(trafo)
## non-standard norms
FI <- distr::solve(trafo%*%matrix(1/Finfo,1,1)%*%t(trafo))
if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") ){
QuadForm(normtype) <- PosSemDefSymmMatrix(FI)
normtype(risk) <- normtype
}
std <- if(is(normtype,"QFNorm"))
QuadForm(normtype) else diag(p)
FI0 <- sum(diag(std%*%as.matrix(FI)))
if(identical(all.equal(radius, 0), TRUE)){
if(warn) cat("'radius == 0' => (classical) optimal IC\n",
"in sense of Cramer-Rao bound is returned\n")
res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(),
neighbor = neighbor, Finfo = Finfo, trafo = trafo,
verbose = verbose)
res <- c(res, list(biastype = biastype, normtype = NormType()))
if(!is(risk, "asMSE")){
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, normtype = normtype,
clip = res$b, cent = res$a,
stand = res$A, trafo = trafo, FI = FI0)
res$risk <- c(Risk, res$risk)
}
Cov <- res$risk$asCov
res$risk$asBias <- list(value = b, biastype = biastype,
normtype = NormType(),
neighbortype = class(neighbor))
res$risk$asMSE <- list(value = Cov + radius^2*b^2,
r = radius,
at = neighbor)
return(res)
}
z <- 0
c0 <- 0
iter <- 0
if(is(symm, "SphericalSymmetry"))
S <- symm@SymmCenter == 0
else
S <- FALSE
### print ---
## assign("l2D",L2deriv,.GlobalEnv)
###
prec <- 1
problem <- FALSE
repeat{
iter <- iter + 1
z.old <- z
c0.old <- c0
## new
if(is(risk,"asMSE")){
L1n <- getL1normL2deriv(L2deriv = L2deriv, cent = z)
lower0 <- L1n/(1 + radius^2)
# if(is(neighbor,"TotalVarNeighborhood")) {
# lower0 <- (L1n-z)/(1 + radius^2)/2}
upper0 <- max(L1n/radius,
sqrt( as.numeric( Finfo + z^2 )/(( 1 + radius^2)^2 - 1) ))
if (is.null(lower))
lower <- .Machine$double.eps^0.75
else {if(iter>1) lower <- min(lower0,2*lower)}
if (is.null(upper))#|(iter == 1))
upper <- getUp(L2deriv)
else {if(iter>1) upper <- max(0.5*upper,3*upper0)}
## print(c(lower,upper))
#lower <- 0; upper <- 100
##
}else{
if(is.null(lower)) lower <- 10^-9
if(is.null(upper)) upper <- 40
}
c0 <- try(uniroot(getInfClip,
## new
lower = lower, upper = upper,
##
tol = tol, L2deriv = L2deriv, risk = risk,
neighbor = neighbor, biastype = biastype,
cent = z, symm = S,
trafo = trafo)$root, silent = TRUE)
if(!is.numeric(c0)){
if(warn) cat("The IC algorithm did not converge!\n",
"'radius >= maximum radius' for the given risk?\n",
"=> the minimum asymptotic bias (lower case) solution is returned\n")
res <- getInfRobIC(L2deriv = L2deriv, risk = asBias(biastype = biastype(risk),
normtype = normtype(risk)),
neighbor = neighbor, Finfo = Finfo,
symm = symm, trafo = trafo, upper = upper,
maxiter = maxiter, tol = tol, warn = warn,
verbose = verbose)
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, normtype = normtype,
clip = res$b, cent = res$a,
stand = res$A, trafo = trafo, FI = FI0)
res$risk <- c(Risk, res$risk)
return(res)
}
z <- getInfCent(L2deriv = L2deriv, neighbor = neighbor, biastype = biastype,
clip = c0, cent = z, symm = S, trafo = trafo, tol.z = tol)
## cat("c0:\t", c0, "c0.old:\t", c0.old, "z:\t", z, "z.old:\t", z.old, "\n")
if(S) break
prec.old <- prec
## print(c(c0,z))
prec <- max(abs(z - z.old), abs(c0-c0.old))
if(iter>1){
if(verbose && iter%%5==1){
cat("current precision in IC algo:\t", prec, "\n")
print(round(c(r=radius,c=c0,z=z),3))
}
}
if(prec < tol) break
if(abs(prec.old - prec) < 1e-10){
if(iter>1)
problem <- TRUE
if(warn) cat("algorithm did not converge!\n", "achieved precision:\t", prec, "\n")
break
}
if(iter > maxiter){
if(iter>1)
problem <- TRUE
if(warn) cat("maximum iterations reached!\n", "achieved precision:\t", prec, "\n")
break
}
}
info <- paste("optimally robust IC for", sQuote(class(risk)[1]))
A <- getInfStand(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, clip = c0, cent = z, trafo = trafo)
a <- as.vector(A)*z
b <- abs(as.vector(A))*c0
if(!is(risk, "asMSE")){
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, normtype = normtype,
clip = b, cent = a, stand = A,
trafo = trafo, FI = FI0)
}else{
Risk <- NULL
}
Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, clip = c0, cent = z, stand = A)
Risk <- c(Risk, list(asCov = Cov,
asBias = list(value = b, biastype = biastype,
normtype = normtype(risk),
neighbortype = class(neighbor)),
trAsCov = list(value = Cov, normtype = normtype(risk)),
asMSE = list(value = Cov + radius^2*b^2,
r = radius,
at = neighbor)))
if(is(neighbor,"ContNeighborhood")){
w <- new("HampelWeight")
clip(w) <- b
cent(w) <- as.numeric(z)
stand(w) <- A
}else if (is(neighbor,"TotalVarNeighborhood")){
w <- new("BdStWeight")
clip(w) <- c(0,b)+a
stand(w) <- A
}
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
normW = NormType())
## print(list(A = A, a = a, b = b))
return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info, w = w,
biastype = biastype, normtype = normtype(risk), problem = problem ))
})
###################################################################################
# multivariate solution G-Risk --- new 10-08-09
###################################################################################
setMethod("getInfRobIC", signature(L2deriv = "RealRandVariable",
risk = "asGRisk",
neighbor = "UncondNeighborhood"),
function(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
L2derivDistrSymm, Finfo, trafo, onesetLM = FALSE,
z.start, A.start, upper = NULL, lower = NULL,
OptOrIter = "iterate",
maxiter, tol, warn, verbose = NULL, withPICcheck = TRUE,
..., .withEvalAsVar = TRUE){
if(missing(verbose)|| is.null(verbose))
verbose <- getRobAStBaseOption("all.verbose")
mc <- match.call()
## some abbreviations / checks
radius <- neighbor@radius
biastype <- biastype(risk)
normtype <- normtype(risk)
p <- nrow(trafo)
k <- ncol(trafo)
if( is(neighbor,"TotalVarNeighborhood") && p>1)
stop("Not yet implemented")
## non-standard norms
FI <- distr::solve(trafo%*%distr::solve(Finfo)%*%t(trafo))
if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") ){
QuadForm(normtype) <- PosSemDefSymmMatrix(FI)
normtype(risk) <- normtype
}
std <- if(is(normtype,"QFNorm"))
QuadForm(normtype) else diag(p)
FI0 <- sum(diag(std%*%as.matrix(FI)))
## starting values
if(is.null(z.start)) z.start <- numeric(k)
if(is.null(A.start)) A.start <- trafo %*% distr::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, radius = radius,
risk = risk, neighbor = neighbor,
biastype = biastype, normtype = normtype,
Distr = Distr, Finfo = Finfo, trafo = trafo,
QuadForm = std, verbose = verbose, warn = warn))
## determine which entries must be computed
# 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(.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
if (is.na(pM)) OptOrIter <- 1
## initialize
if(is(neighbor,"ContNeighborhood")){
w <- new("HampelWeight")
}else{ if (is(neighbor,"TotalVarNeighborhood"))
w <- new("BdStWeight")
}
z <- z.start
A <- A.start
b <- 0
a <- a.start
iter <- 0
prec <- 1
iter.In <- 0
problem <- FALSE
## determining A,a,b with either optimization of iteration:
if(OptOrIter == 1){
if(is.null(lower)){
lowBerg <- .getLowerSol(L2deriv = L2deriv, risk = risk,
neighbor = neighbor, Distr = Distr,
DistrSymm = DistrSymm,
L2derivSymm = L2derivSymm,
L2derivDistrSymm = L2derivDistrSymm,
z.start = z.start, A.start = A.start,
trafo = trafo, maxiter = maxiter, tol = tol,
warn = FALSE, Finfo = Finfo, QuadForm = std,
verbose = verbose, ...)
lower <- lowBerg$b}
#if(is.null(upper))
upper <- 5*max(distr::solve(Finfo))
OptIterCall <- numeric(1)
Cov <- 0
Risk <- 1e10
normtype.old <- normtype
a <- as.numeric(A %*% z)
normtype.opt <- normtype
# asGRiskb <- function(b0){
# iter <<- iter + 1
erg <- getLagrangeMultByOptim(b = 1, L2deriv = L2deriv, risk = risk,
FI = Finfo, trafo = trafo, neighbor = neighbor,
biastype = biastype, normtype = normtype, Distr = Distr,
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),
verbose = verbose, ...)
w <- erg$w
A <- erg$A
a <- erg$a
z <- erg$z
b <- erg$b
OptIterCall <- erg$call
std <- if(is.null(erg$std)) std else erg$std
biastype <- erg$biastype
normtype.old <- erg$normtype.old
normtype <- erg$normtype
std.n <- if(is(normtype,"QFNorm"))
QuadForm(normtype) else diag(p)
FI0 <- sum(diag(std.n%*%as.matrix(FI)))
risk <- erg$risk
iter <- erg$iter
prec <- prec.In <- iter.In <- NULL
}else{
repeat{
iter <- iter + 1
a.old <- a
z.old <- z
b.old <- b
A.old <- A
##
## get interval for b-search:
LUB <- .getLowUpB(L2deriv = L2deriv, Finfo = Finfo, Distr = Distr,
normtype = normtype, z = z, A = A, radius = radius,
iter = iter)
if (!is.null(upper)|(iter == 1))
{lower <- .Machine$double.eps^0.75;
if(is.null(upper)) upper <- 10*LUB$upper
}else{ lower <- LUB$lower; upper <- LUB$upper}
##
# lower <- LUB$lower
# upper <- if(is.null(upper)) LUB$upper else min(upper,LUB$upper)
# print(c(lower,upper))
## solve for b
b <- try(uniroot(getInfClip,
lower = lower, upper = upper,
tol = tol, L2deriv = L2deriv, risk = risk,
biastype = biastype, Distr = Distr, neighbor = neighbor,
stand = A, cent = z, trafo = trafo, ...)$root, silent = TRUE)
## if no solution return lower case:
if(!is.numeric(b))
return(.getLowerSol(L2deriv = L2deriv, risk = risk,
neighbor = neighbor, Distr = Distr,
DistrSymm = DistrSymm,
L2derivSymm = L2derivSymm,
L2derivDistrSymm = L2derivDistrSymm,
z.start = z.start, A.start = A.start,
trafo = trafo, maxiter = maxiter, tol = tol,
warn = warn, Finfo = Finfo, QuadForm = std,
verbose = verbose, ...)
)
maxit2 <- if(OptOrIter==2) 0 else maxiter
erg <- getLagrangeMultByIter(b = b, L2deriv = L2deriv, risk = risk,
trafo = trafo, neighbor = neighbor, biastype = biastype,
normtype = normtype, Distr = Distr,
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 = warn&(OptOrIter!=2))
## read out solution
w <- erg$w
A <- erg$A
a <- erg$a
z <- erg$z
biastype <- erg$biastype
normtype.old <- erg$normtype.old
normtype <- erg$normtype
risk <- erg$risk
iter.In <- iter.In + erg$iter
prec.In <- erg$prec
OptIterCall <- erg$call
std <- if(is.null(erg$std)) std else erg$std
# print(list(z=z,A=A,b=b))
FI0 <- sum(diag(std%*%as.matrix(FI)))
## 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), max(abs(a-a.old))))
if(verbose && iter%%5==1){
cat("current precision in IC algo:\t", prec, "\n")
print(round(c(r=radius,b=b,A=A,a=a),3))
}
if(prec < tol) break
if(abs(prec.old - prec) < 1e-10){
problem <- TRUE
if(warn) cat("algorithm did not converge!\n", "achieved precision:\t", prec, "\n")
break
}
if(iter > maxiter){
problem <- TRUE
if(warn) cat("maximum iterations reached!\n", "achieved precision:\t", prec, "\n")
break
}
}
if (onesetLM){
if(is(neighbor,"ContNeighborhood"))
cent(w) <- as.numeric(z)
if(is(neighbor,"TotalVarNeighborhood"))
clip(w) <- c(0,b)+a
stand(w) <- A
weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
normW = normtype)
}
else normtype <- normtype.old
### issue some diagnostics if wanted
if(verbose && withPICcheck){
cat("Iterations needed: outer (b-loop):",
iter," inner (A,a-loop):", iter.In,"\n")
cat("Precision achieved: all in all (b+A,a-loop):",
prec," inner (A,a-loop):", prec.In,"\n")
}
if(verbose && withPICcheck) print(list(A=A,a=a,w=w))
}
Cov <- substitute(do.call(getInfV, args = c(list(L2deriv = L2deriv0,
neighbor = neighbor0, biastype = biastype0,
Distr = Distr0, V.comp = A.comp0, cent = a0,
stand = A0, w = w0), dots0)), list(L2deriv0 = L2deriv,
neighbor0 = neighbor, biastype0 = biastype,
Distr0 = Distr, A.comp0 = A.comp, a0 = a,
A0 = A, w0 = w, dots0=list(...)))
rifct <- function(std0, Cov0, rad0, b0){
sum(diag(std0%*%eval(Cov0))) + rad0^2 * b0^2}
asMSE.0 <- substitute(do.call(ri.fct, args=list(std0=std1, Cov0=Cov1,
rad0 = rad1, b0=b1)), list(ri.fct = rifct,
std1=std, Cov1=Cov, rad1=radius, b1=b))
if(!is(risk, "asMSE")){
Risk <- substitute(do.call(getAsRisk, args =list(risk = risk0,
L2deriv = L2deriv0, neighbor = neighbor0,
biastype = biastype0, normtype = normtype0,
clip = b0, cent = a0, stand = A0,
trafo = trafo0, FI = FI000)), list(risk0=risk,
L2deriv0=L2deriv, neighbor0 = neighbor,
biastype0 = biastype, normtype0 = normtype,
b0 = b, a0 = a, A0 = A,
trafo0 = trafo, FI000 = FI0))
}else{ Risk <- asMSE.0
}
### add some further informations for the pIC-slots info and risk
info <- paste("optimally robust IC for", sQuote(class(risk)[1]))
trAsCov.fct <- function(std0, Cov0) sum(diag(std0%*%eval(Cov0)))
trAsCov <- substitute(do.call(tr.fct, args=list(std0=std1, Cov0=Cov1)),
list(tr.fct = trAsCov.fct, std1=std, Cov1=Cov))
Risk <- c(Risk, list(asCov = Cov,
asBias = list(value = b, biastype = biastype,
normtype = normtype,
neighbortype = class(neighbor)),
trAsCov = list(value = trAsCov,
normtype = normtype),
asMSE = list(value = asMSE.0,
r = radius,
at = neighbor)))
if(.withEvalAsVar) Risk <- .evalListRec(Risk)
if(verbose && withPICcheck)
.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,
iter.In = iter.In, prec.In = prec.In, problem = problem ))
})
### helper function to return the upper case solution if r=0
.getUpperSol <- function(L2deriv, radius, risk, neighbor, biastype,
normtype, Distr, Finfo, trafo,
QuadForm, verbose, warn, ...){
if(warn) cat("'radius == 0' => (classical) optimal IC\n",
"in sense of Cramer-Rao bound is returned\n")
res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(), neighbor = neighbor,
Distr = Distr, Finfo = Finfo, trafo = trafo,
QuadForm = QuadForm, verbose = verbose)
b <- res$b
res <- c(res, list(biastype = biastype, normtype = normtype))
if(!is(risk, "asMSE")){
FI <- trafo%*%distr::solve(Finfo)%*%t(trafo)
FI <- sum(diag(QuadForm %*% FI))
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype, normtype = normtype,
cent = res$a, stand = res$A, trafo = trafo,
FI = FI, ...)
res$risk <- c(Risk, res$risk)
}
trAsCov <- sum(diag(QuadForm%*%res$risk$asCov));
res$risk$trAsCov <- list(value = trAsCov, normtype = normtype)
res$risk$asBias <- list(value = b, biastype = biastype,
normtype = normtype,
neighbortype = class(neighbor))
res$risk$asMSE <- list(value = trAsCov + radius^2*b^2,
r = radius,
at = neighbor)
return(res)
}
### helper function to return the lower case solution if b-search was not successful
.getLowerSol <- function(L2deriv, risk, neighbor, Distr, DistrSymm,
L2derivSymm, L2derivDistrSymm,
z.start, A.start, trafo,
maxiter, tol, warn, Finfo, QuadForm, verbose,...){
if(warn) cat("Could not determine optimal clipping bound!\n",
"'radius >= maximum radius' for the given risk?\n",
"=> the minimum asymptotic bias (lower case) solution is returned\n",
"If 'no' => Try again with modified starting values ",
"'z.start' and 'A.start'\n")
res <- getInfRobIC(L2deriv = L2deriv,
risk = asBias(biastype = biastype(risk),
normtype = normtype(risk)),
neighbor = neighbor, Distr = Distr, DistrSymm = DistrSymm,
L2derivSymm = L2derivSymm, L2derivDistrSymm = L2derivDistrSymm,
z.start = z.start, A.start = A.start, trafo = trafo,
maxiter = round(maxiter), tol = tol, warn = warn, Finfo = Finfo,
verbose = verbose,...)
normtype(risk) <- res$normtype
if(!is(risk, "asMSE")){
FI <- trafo%*%distr::solve(Finfo)%*%t(trafo)
FI <- sum(diag(QuadForm %*% FI))
Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
biastype = biastype(risk), normtype = normtype(risk),
clip = NULL, cent = res$a, stand = res$A, Distr = Distr,
FI = FI,w = res$w)
res$risk <- c(Risk, res$risk)
}
return(res)
}
### helper function to return upper & lower bounds for b for b-search
.getLowUpB <- function(L2deriv, Finfo, Distr, normtype, z, A, radius, iter){
L1n <- getL1normL2deriv(L2deriv = L2deriv, cent = z, stand = A,
Distr = Distr, normtype = normtype)
lower0 <- L1n/(1+radius^2)
if(is(neighbor,"TotalVarNeighborhood")) {
lower0 <- (L1n-A%*%z)/(1 + radius^2)/2}
QF <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(nrow(A))
upper0 <- max(L1n/radius,
sqrt( (sum( diag(QF%*%A%*%Finfo%*%t(A))) + t(A%*%z)%*%QF%*%(A%*%z)) /
((1 + radius^2)^2-1)))
if (iter == 1){
lower <- .Machine$double.eps^.6
upper <- 10*upper0
}else{
lower <- lower0
upper <- 2*upper0
}
##
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, ...){
dotsI <- .filterEargsWEargList(list(...))
if(is.null(dotsI$useApply)) dotsI$useApply <- FALSE
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] <- do.call(E,c(list(Distr,integrand0,ixx=i), dotsI))
cent1 <- numeric(pvalues)
for(i in 1:pvalues)
cent1[i] <- do.call(E,c(list(Distr,integrand1,ixx=i), dotsI))
consist <- 0*trafo
for(i in 1:pvalues){
for(j in 1:nrvalues){
if(A.comp[i,j])
consist[i,j] <- do.call(E,c(list(Distr,integrand2,ixx=i,jxx=j), dotsI))
}
}
consist <- consist-trafo
cat("centering (k-space):",cent0,"\n")
cat("centering (p-space):",cent1,"\n")
cat("Fisher consistency:\n")
print(consist)
}
################################################################################
## old routine
################################################################################
#
#setMethod("getInfRobIC", signature(L2deriv = "RealRandVariable",
# risk = "asGRisk",
# neighbor = "UncondNeighborhood"),
# function(L2deriv, risk, neighbor, Distr, DistrSymm, L2derivSymm,
# L2derivDistrSymm, Finfo, trafo, onesetLM = FALSE,
# z.start, A.start, upper = NULL, maxiter, tol, warn, verbose = FALSE){
# biastype <- biastype(risk)
# normtype <- normtype(risk)
# p <- nrow(trafo)
# if( is(neighbor,"TotalVarNeighborhood") && p>1)
# stop("Not yet implemented")
#
# FI <- distr::solve(trafo%*%distr::solve(Finfo)%*%t(trafo))
# if(is(normtype,"InfoNorm") || is(normtype,"SelfNorm") )
# {QuadForm(normtype) <- PosSemDefSymmMatrix(FI);
# normtype(risk) <- normtype}
# QF <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(nrow(trafo))
#
# if(is.null(z.start)) z.start <- numeric(ncol(trafo))
# if(is.null(A.start)) A.start <- trafo %*% distr::solve(Finfo)
#
# radius <- neighbor@radius
# if(identical(all.equal(radius, 0), TRUE)){
# if(warn) cat("'radius == 0' => (classical) optimal IC\n",
# "in sense of Cramer-Rao bound is returned\n")
# res <- getInfRobIC(L2deriv = L2deriv, risk = asCov(), neighbor = neighbor,
# Distr = Distr, Finfo = Finfo, trafo = trafo,
# QuadForm = QF, verbose = verbose)
# res <- c(res, list(biastype = biastype, normtype = normtype))
# if(!is(risk, "asMSE")){
# Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, cent = res$a,
# stand = res$A, trafo = trafo)
# res$risk <- c(Risk, res$risk)
# }
# trAsCov <- sum(diag(QF%*%res$risk$asCov));
# res$risk$trAsCov <- list(value = trAsCov, normtype = normtype)
# res$risk$asBias <- list(value = b, biastype = biastype,
# normtype = normtype,
# neighbortype = class(neighbor))
# res$risk$asMSE <- list(value = trAsCov + radius^2*b^2,
# r = radius,
# at = neighbor)
# return(res)
# }
#
# comp <- .getComp(L2deriv, DistrSymm, L2derivSymm, L2derivDistrSymm)
#
# z.comp <- comp$"z.comp"
# A.comp <- comp$"A.comp"
#
# if(is(neighbor,"ContNeighborhood")){
# w <- new("HampelWeight")
# }else{ if (is(neighbor,"TotalVarNeighborhood"))
# w <- new("BdStWeight")
# }
# z <- z.start
# A <- A.start
# b <- 0
# a <- 0
# iter <- 0
# prec <- 1
# repeat{
# iter <- iter + 1
# z.old <- z
# b.old <- b
# A.old <- A
# ##
# if(is(neighbor,"ContNeighborhood"))
# cent(w) <- z
# stand(w) <- A
#
# ## new
# L1n <- getL1normL2deriv(L2deriv = L2deriv, cent = z, stand = A,
# Distr = Distr, normtype = normtype)
# lower0 <- L1n/(1+radius^2)
# if(is(neighbor,"TotalVarNeighborhood")) {
# lower0 <- (L1n-A%*%z)/(1 + radius^2)/2}
#
# QF <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(nrow(A))
# upper0 <- max(L1n/radius,
# sqrt( (sum( diag(QF%*%A%*%Finfo%*%t(A))) + t(A%*%z)%*%QF%*%(A%*%z)) /
# ((1 + radius^2)^2-1)))
#
# if (!is.null(upper)|(iter == 1))
# {lower <- .Machine$double.eps^0.6;
# if(is.null(upper)) upper <- 10*upper0
# }else{ lower <- lower0; upper <- upper0}
#
# ##
#
# b <- try(uniroot(getInfClip,
# ## new
# lower = lower, upper = upper,
# ##
# tol = tol, L2deriv = L2deriv, risk = risk,
# biastype = biastype, Distr = Distr, neighbor = neighbor,
# stand = A, cent = z, trafo = trafo)$root, silent = TRUE)
# if(!is.numeric(b)){
# if(warn) cat("Could not determine optimal clipping bound!\n",
# "'radius >= maximum radius' for the given risk?\n",
# "=> the minimum asymptotic bias (lower case) solution is returned\n",
# "If 'no' => Try again with modified starting values ",
# "'z.start' and 'A.start'\n")
# res <- getInfRobIC(L2deriv = L2deriv,
# risk = asBias(biastype = biastype(risk),
# normtype = normtype(risk)),
# neighbor = neighbor, Distr = Distr, DistrSymm = DistrSymm,
# L2derivSymm = L2derivSymm, L2derivDistrSymm = L2derivDistrSymm,
# z.start = z.start, A.start = A.start, trafo = trafo,
# maxiter = maxiter, tol = tol, warn = warn, Finfo = Finfo,
# verbose = verbose)
# normtype(risk) <- res$normtype
# if(!is(risk, "asMSE")){
# Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, clip = NULL,
# cent = res$a, stand = res$A, trafo = trafo)
# res$risk <- c(Risk, res$risk)
# }
# return(res)
# }
# if(is(neighbor,"ContNeighborhood")){
# clip(w) <- b
# }else if(is(neighbor,"TotalVarNeighborhood")){
# clip(w) <- if(iter==1) c(-1,1)*b/2 else c(0,b)+a
#
# }
#
#
# weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
# normW = normtype)
#
# z <- getInfCent(L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, Distr = Distr, z.comp = z.comp,
# w = w)
#
# if(is(neighbor,"TotalVarNeighborhood")){
# a <- z
# z <- distr::solve(A,a)
# zc <- numeric(ncol(trafo))
# }else if(is(neighbor,"ContNeighborhood")) {
# zc <- z
# }
#
# A <- getInfStand(L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, Distr = Distr, A.comp = A.comp,
# cent = zc, trafo = trafo, w = w)
#
# 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 = as.vector(A %*% z),
# stand = A, w = w)}
#
# prec.old <- prec
# prec <- max(abs(b-b.old), max(abs(A-A.old)), max(abs(z-z.old)))
# if(verbose)
# cat("current precision in IC algo:\t", prec, "\n")
# if(prec < tol) break
# if(abs(prec.old - prec) < 1e-10){
# cat("algorithm did not converge!\n", "achieved precision:\t", prec, "\n")
# break
# }
# if(iter > maxiter){
# cat("maximum iterations reached!\n", "achieved precision:\t", prec, "\n")
# break
# }
# }
# if (onesetLM){
# if(is(neighbor,"ContNeighborhood"))
# cent(w) <- z
# if(is(neighbor,"TotalVarNeighborhood"))
# clip(w) <- c(0,b)+a
# stand(w) <- A
# weight(w) <- getweight(w, neighbor = neighbor, biastype = biastype,
# normW = normtype)
# }else normtype <- normtype.old
#
# if(is(neighbor,"ContNeighborhood"))
# a <- as.vector(A %*% z)
#
# info <- paste("optimally robust IC for", sQuote(class(risk)[1]))
# if(!is(risk, "asMSE")){
# Risk <- getAsRisk(risk = risk, L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, clip = b, cent = a, stand = A,
# trafo = trafo)
# }else{
# Risk <- NULL
# }
# Cov <- getInfV(L2deriv = L2deriv, neighbor = neighbor,
# biastype = biastype, Distr = Distr,
# V.comp = A.comp, cent = a,
# stand = A, w = w)
#
# QF <- if(is(normtype,"QFNorm")) QuadForm(normtype) else diag(nrow(A))
#
# trAsCov <- sum(diag(QF%*%Cov))
# Risk <- c(Risk, list(asCov = Cov,
# asBias = list(value = b, biastype = biastype,
# normtype = normtype,
# neighbortype = class(neighbor)),
# trAsCov = list(value = trAsCov,
# normtype = normtype),
# asMSE = list(value = trAsCov + radius^2*b^2,
# r = radius,
# at = neighbor)))
#
# return(list(A = A, a = a, b = b, d = NULL, risk = Risk, info = info, w = w,
# biastype = biastype, normtype = normtype))
# })
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.