[Robast-commits] r42 - pkg/RobLox/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 19 06:59:47 CET 2008
Author: stamats
Date: 2008-02-19 06:59:47 +0100 (Tue, 19 Feb 2008)
New Revision: 42
Modified:
pkg/RobLox/R/rlOptIC.R
pkg/RobLox/R/rlsOptIC_AL.R
pkg/RobLox/R/rlsOptIC_An1.R
pkg/RobLox/R/rlsOptIC_An2.R
pkg/RobLox/R/rlsOptIC_AnMad.R
pkg/RobLox/R/rlsOptIC_BM.R
pkg/RobLox/R/rlsOptIC_Ha3.R
pkg/RobLox/R/rlsOptIC_Ha4.R
pkg/RobLox/R/rlsOptIC_HaMad.R
pkg/RobLox/R/rlsOptIC_Hu1.R
pkg/RobLox/R/rlsOptIC_Hu2.R
pkg/RobLox/R/rlsOptIC_Hu2a.R
pkg/RobLox/R/rlsOptIC_Hu3.R
pkg/RobLox/R/rlsOptIC_HuMad.R
pkg/RobLox/R/rlsOptIC_M.R
pkg/RobLox/R/rlsOptIC_MM2.R
pkg/RobLox/R/rlsOptIC_Tu1.R
pkg/RobLox/R/rlsOptIC_Tu2.R
pkg/RobLox/R/rlsOptIC_TuMad.R
pkg/RobLox/R/roblox.R
pkg/RobLox/R/rsOptIC.R
Log:
just formatting
Modified: pkg/RobLox/R/rlOptIC.R
===================================================================
--- pkg/RobLox/R/rlOptIC.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlOptIC.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,14 +1,14 @@
###############################################################################
-# optimally robust IC for normal location
+## optimally robust IC for normal location
###############################################################################
rlOptIC <- function(r, mean = 0, sd = 1, bUp = 1000, computeIC = TRUE){
c0 <- uniroot(function(c0, r){return(r^2*c0 - 2*(dnorm(c0) - c0*pnorm(-c0)))},
lower = 1e-4, upper = bUp, tol = .Machine$double.eps^0.5, r = r)$root
-
+
A1 <- 1/(2*pnorm(c0) - 1)
b <- sd*A1*c0
A <- sd^2*A1
-
+
if(computeIC){
return(generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormLocationFamily(mean = mean, sd = sd),
Modified: pkg/RobLox/R/rlsOptIC_AL.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_AL.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_AL.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,15 +1,15 @@
###############################################################################
-# weight function
+## weight function
###############################################################################
.ALrlsGetw <- function(x, b, a1, a2, a3){
hvkt <- sqrt(a3^2*x^4 + (a1^2 - 2*a2*a3^2)*x^2 + a2^2*a3^2)
ind1 <- (hvkt < b)
-
+
return(ind1 + (1-ind1)*b/hvkt)
}
###############################################################################
-# computation of r
+## computation of r
###############################################################################
.ALrlsGetr <- function(b, r, a1, a2, a3){
integrandr <- function(x, b, a1, a2, a3){
@@ -25,7 +25,7 @@
###############################################################################
-# computation of a1, a2 and a3
+## computation of a1, a2 and a3
###############################################################################
.ALrlsGeta1a2a3 <- function(b, a1, a2, a3, inteps=1e-12){
integrand1 <- function(x, b, a1, a2, a3){
@@ -43,7 +43,7 @@
rel.tol = .Machine$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)$value
a2 <- Int1/Int2
-
+
integrand3 <- function(x, b, a1, a2, a3){
(x^2 - a2)^2*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
@@ -57,7 +57,7 @@
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.AL <- function(r, mean = 0, sd = 1, A.loc.start = 1, a.sc.start = 0,
A.sc.start = 0.5, bUp = 1000, delta = 1e-6, itmax = 100,
@@ -68,7 +68,7 @@
a3 = a3)$root
iter <- 0
- repeat{
+ repeat{
iter <- iter + 1
if(iter > itmax){
stop("Algorithm did not converge!\n",
@@ -76,7 +76,7 @@
"'A.loc.start', 'a.sc.start' and 'A.sc.start'\n")
}
a1.old <- a1; a2.old <- a2; a3.old <- a3; b.old <- b
-
+
a1a2a3 <- .ALrlsGeta1a2a3(b = b, a1 = a1, a2 = a2, a3 = a3)
a1 <- a1a2a3$a1; a2 <- a1a2a3$a2; a3 <- a1a2a3$a3
@@ -103,7 +103,7 @@
rel.tol = .Machine$double.eps^0.5, b = b, a1 = a1,
a2 = a2, a3 = a3)$value
ch2 <- a3*Int2
-
+
integrand3 <- function(x, b, a1, a2, a3){
(x^2 - a2)*.ALrlsGetw(x, b, a1, a2, a3)*dnorm(x)
}
@@ -113,19 +113,19 @@
ch3 <- a3*Int3
ch4 <- .ALrlsGetr(b = b, r = r, a1 = a1, a2 = a2, a3 = a3)
-
+
cat("Fisher consistency of eta.loc:\t", ch1-1, "\n")
cat("centering of eta.sc:\t", ch3, "\n")
cat("Fisher consistency of eta.sc:\t", ch2-1, "\n")
cat("MSE equation:\t", ch4, "\n")
}
-
+
A <- sd^2*diag(c(a1, a3))
a <- sd*c(0, a3*(a2-1))
b <- sd*b
mse <- sd^2*(a1 + a3)
-
+
if(computeIC){
return(generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormLocationScaleFamily(mean = mean, sd = sd),
Modified: pkg/RobLox/R/rlsOptIC_An1.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_An1.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_An1.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,10 +1,10 @@
###############################################################################
-# psi function
+## psi function
###############################################################################
lsAn1.psi <- function(x, a){ sin(x/a)*(abs(x) < a*pi) }
###############################################################################
-# chi function
+## chi function
###############################################################################
lsAn1.chi <- function(x, a){
A.loc <- 1/(2*integrate(f = function(x, a0){cos(x/a0)*dnorm(x)/a0}, lower = 0,
@@ -14,7 +14,7 @@
}
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.An1rlsGetbias <- function(x, a){
A.loc <- 1/(2*integrate(f = function(x, a0){cos(x/a0)*dnorm(x)/a0}, lower = 0,
@@ -24,13 +24,13 @@
Int2 <- integrate(f = function(x, a0){x^2*cos(x/a0)/a0*dnorm(x)}, lower = 0,
upper = a*pi, rel.tol = .Machine$double.eps^0.5, a0 = a)$value
A.sc <- 1/(2*A.loc*(Int1 + Int2))
-
+
return(sqrt(A.loc^2*(sin(x/a)*(abs(x) < a*pi))^2
+ A.sc^2*(A.loc*x*sin(x/a)*(abs(x) < a*pi) - 1)^2))
}
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.An1rlsGetvar <- function(a){
A.loc <- 1/(2*integrate(f = function(x, a0){cos(x/a0)*dnorm(x)/a0}, lower = 0,
@@ -50,11 +50,11 @@
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.An1rlsGetmse <- function(a, r){
Var <- .An1rlsGetvar(a = a)
-
+
x <- seq(from=0, to=a*pi, by = 0.01)
bias <- sapply(x, .An1rlsGetbias, a = a)
index <- which.max(bias)
@@ -66,7 +66,7 @@
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.An1 <- function(r, aUp = 2.5, delta = 1e-6){
res <- optimize(f = .An1rlsGetmse, lower = 1e-4, upper = aUp,
@@ -94,7 +94,7 @@
fct2 <- function(x){ A.sc*(A.loc*x*sin(x/a)*(abs(x) < a*pi) - 1) }
body(fct2) <- substitute({ A.sc*(A.loc*x*sin(x/a)*(abs(x) < a*pi) - 1) },
list(a = a, A.loc = A.loc, A.sc = A.sc))
-
+
return(IC(name = "IC of An1 type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$objective, asBias = b, asCov = res$objective - r^2*b^2),
Modified: pkg/RobLox/R/rlsOptIC_An2.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_An2.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_An2.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.An2rlsGetbias <- function(x, a, k){
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
@@ -12,7 +12,7 @@
}
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.An2rlsGetvar <- function(a, k){
h1 <- 2*integrate(f=function(x, a0){ sin(x/a0)^2*dnorm(x) }, lower = 0,
@@ -29,11 +29,11 @@
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.An2rlsGetmse <- function(ak, r, MAX){
a <- ak[1]; k <- ak[2]
-
+
# constraints
if(a < 0 || k < 0) return(MAX)
@@ -50,14 +50,14 @@
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.An2 <- function(r, a.start = 1.5, k.start = 1.5, delta = 1e-6, MAX = 100){
res <- optim(c(a.start, k.start), .An2rlsGetmse, method = "Nelder-Mead",
control = list(reltol=delta), r = r, MAX = MAX)
a <- res$par[1]; k <- res$par[2]
-
+
A.loc <- 1/(2*integrate(f = function(x, a0){ cos(x/a0)*dnorm(x)/a0 },
lower = 0, upper = a*pi, rel.tol = .Machine$double.eps^0.5,
a0 = a)$value)
@@ -70,7 +70,7 @@
fct1 <- function(x){ A.loc*sin(x/a)*(abs(x) < a*pi) }
body(fct1) <- substitute({ A.loc*sin(x/a)*(abs(x) < a*pi) },
- list(a = a, A.loc = A.loc))
+ list(a = a, A.loc = A.loc))
A.sc <- 1/(2*(2*pnorm(k) - 1) - 4*k*dnorm(k))
fct2 <- function(x){ A.sc*(pmin(x^2, k^2) - beta.k) }
body(fct2) <- substitute({ A.sc*(pmin(x^2, k^2) - beta.k) },
Modified: pkg/RobLox/R/rlsOptIC_AnMad.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_AnMad.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_AnMad.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.AnMadrlsGetvar <- function(a){
h1 <- 2*integrate(f = function(x, a0){ sin(x/a0)^2*dnorm(x) }, lower = 0,
@@ -13,7 +13,7 @@
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.AnMadrlsGetmse <- function(a, r){
A.loc <- 1/(2*integrate(f = function(x, a0){ cos(x/a0)*dnorm(x)/a0 }, lower = 0,
@@ -25,7 +25,7 @@
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.AnMad <- function(r, aUp = 2.5, delta = 1e-6){
res <- optimize(f = .AnMadrlsGetmse, lower = 1e-4, upper = aUp,
@@ -44,12 +44,12 @@
fct2 <- function(x){ b.mad*sign(abs(x) - a.mad) }
body(fct2) <- substitute({ b.mad*sign(abs(x) - a.mad) },
list(a.mad = a.mad, b.mad = b.mad))
-
+
return(IC(name = "IC of AnMad type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$objective, asBias = bias, asCov = res$objective - r^2*bias^2),
Infos = matrix(c("rlsOptIC.AnMad", "optimally robust IC for AnMad estimators and 'asMSE'",
"rlsOptIC.AnMad", paste("where a =", round(a, 3))),
ncol=2, byrow = TRUE, dimnames=list(character(0), c("method", "message"))),
- CallL2Fam = call("NormLocationScaleFamily")))
+ CallL2Fam = call("NormLocationScaleFamily")))
}
Modified: pkg/RobLox/R/rlsOptIC_BM.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_BM.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_BM.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,9 +1,9 @@
###############################################################################
-# computation of a
+## computation of a
###############################################################################
.BMrlsGeta <- function(a, bL, bS){
tau1 <- bL/a; tau2 <- (1+bS)/bL; tau3 <- sqrt((1+bS)/a)
-
+
if(tau1 <= tau3)
return(2*(a*(pnorm(tau1)-0.5) - bL*dnorm(tau2)
+ (1+bS)*pnorm(-tau2)) - 1)
@@ -14,11 +14,11 @@
###############################################################################
-# computation of gg
+## computation of gg
###############################################################################
.BMrlsGetgg <- function(a, bL, bS){
tau1 <- bL/a; tau2 <- (1+bS)/bL; tau3 <- sqrt((1+bS)/a)
-
+
if(tau1 <= tau3)
return(1/(2*(-bL*dnorm(tau1) + 3*a*(pnorm(tau1)-0.5)
- 2*bL*dnorm(tau2) + (1+bS)*pnorm(-tau2)) - 1))
@@ -29,11 +29,11 @@
###############################################################################
-# computation of asymptotic variance E|rho|^2
+## computation of asymptotic variance E|rho|^2
###############################################################################
.BMrlsGetvar <- function(bL, bS, a, gg){
tau1 <- bL/a; tau2 <- (1+bS)/bL; tau3 <- sqrt((1+bS)/a)
-
+
if(tau1 <= tau3){
res1 <- 2*(-a*bL*dnorm(tau1) + a^2*(pnorm(tau1)-0.5) + bL^2*(pnorm(tau2)-pnorm(tau1))
+ bL*(1+bS)*dnorm(tau2) - (1+bS)^2*pnorm(-tau2))
@@ -45,13 +45,13 @@
res2 <- 2*(-a^2*tau3^3*dnorm(tau3) + 3*a^2*(-tau3*dnorm(tau3) + pnorm(tau3)-0.5)
+ (1+bS)^2*pnorm(-tau3))
}
-
+
return(res1 + gg^2*(res2 - 1))
}
###############################################################################
-# computation of minimal bL
+## computation of minimal bL
###############################################################################
.BMrlsGetbLmin <- function(bL, bS){
tau <- (1+bS)/bL
@@ -60,17 +60,17 @@
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.BMrlsGetmse <- function(bL.bS, r, MAX){
bL <- bL.bS[1]; bS <- bL.bS[2]
-
+
# constraint for bS
if(bS < 1.0) return(MAX)
-
+
bL.min <- uniroot(.BMrlsGetbLmin, lower = sqrt(pi/2), upper = 10,
tol = .Machine$double.eps^0.5, bS = bS)$root
-
+
# constraint for bL
if(bL < bL.min) return(MAX)
@@ -79,7 +79,7 @@
tau <- (1+bS)/bL
gg <- 1/(4*bL*(1/sqrt(2*pi)-dnorm(tau) + 0.5*tau*pnorm(-tau)) - 1)
bias.2 <- bL^2 + gg^2*bS^2
-
+
Var <- 2*((1+gg^2)*bL^2*(pnorm(tau) - 0.5) + (1-gg^2)*bL*(1+bS)*dnorm(tau)
- (1-gg^2)*(1+bS)^2*pnorm(-tau)) - gg^2
}else{
@@ -87,10 +87,10 @@
tol = .Machine$double.eps^0.5, bL = bL, bS = bS)$root
gg <- .BMrlsGetgg(a = a, bL = bL, bS = bS)
bias.2 <- bL^2 + gg^2*bS^2
-
+
Var <- .BMrlsGetvar(bL = bL, bS = bS, a = a, gg = gg)
}
-
+
return(Var + r^2*bias.2)
}
@@ -102,10 +102,10 @@
res <- optim(c(bL.start, bS.start), .BMrlsGetmse, method = "Nelder-Mead",
control = list(reltol=delta), r = r, MAX = MAX)
bL <- res$par[1]; bS <- res$par[2]
-
+
bL.min <- uniroot(.BMrlsGetbLmin, lower = sqrt(pi/2), upper = 10,
tol = .Machine$double.eps^0.5, bS = bS)$root
-
+
if(abs(bL-bL.min) < 1e-4){
bL <- bL.min
tau <- (1+bS)/bL
@@ -132,7 +132,7 @@
body(fct2) <- substitute({ gg*(abs(x)*pmin(a*abs(x), bL, (1+bS)/abs(x)) - 1) },
list(bL = bL, bS = bS, a = a, gg = gg))
}
-
+
b <- sqrt(bL^2 + gg^2*bS^2)
return(IC(name = "IC of BM type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Modified: pkg/RobLox/R/rlsOptIC_Ha3.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Ha3.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Ha3.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# psi function
+## psi function
###############################################################################
.Ha3rlsGetpsi <- function(x, a, b, c0){
Ind1 <- (abs(x) < a); Ind2 <- (abs(x) < b); Ind3 <- (abs(x) < c0)
@@ -8,7 +8,7 @@
}
###############################################################################
-# chi function
+## chi function
###############################################################################
lsHa3.chi <- function(x, a, b, c0){
A.loc <- 1/(2*(pnorm(a) - 0.5 - a/(c0-b)*(pnorm(c0)-pnorm(b))))
@@ -17,20 +17,20 @@
}
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.Ha3rlsGetbias <- function(x, a, b, c0){
A.loc <- 1/(2*(pnorm(a) - 0.5 - a/(c0-b)*(pnorm(c0)-pnorm(b))))
A.sc <- 1/(2*A.loc*(2*(-a*dnorm(a) + pnorm(a) - 0.5) - a*(dnorm(b)-dnorm(a))
+ a/(c0-b)*(c0*dnorm(c0)+(c0-2*b)*dnorm(b) + 2*(pnorm(b) - pnorm(c0)))))
-
+
return(sqrt(.Ha3rlsGetpsi(x = x, a = a, b = b, c0 = c0)^2*A.loc^2
+ (x*.Ha3rlsGetpsi(x = x, a = a, b = b, c0 = c0)*A.loc - 1)^2*A.sc^2))
}
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Ha3rlsGetvar <- function(a, b, c0){
h1 <- 2*(-a*dnorm(a) + pnorm(a) - 0.5 + a^2*(pnorm(b)-pnorm(a))
@@ -39,7 +39,7 @@
+ pnorm(c0) - pnorm(b)))
A.loc <- 1/(2*(pnorm(a) - 0.5 - a/(c0-b)*(pnorm(c0)-pnorm(b))))
Var.loc <- h1*A.loc^2
-
+
A.sc <- 1/(2*A.loc*(2*(-a*dnorm(a) + pnorm(a) - 0.5) - a*(dnorm(b)-dnorm(a))
+ a/(c0-b)*(c0*dnorm(c0)+(c0-2*b)*dnorm(b)
+ 2*(pnorm(b) - pnorm(c0)))))
@@ -48,21 +48,21 @@
+ (c0^2*b-2*c0*b^2-4*c0+b^3+3*b)*dnorm(b)
+ (c0^2+3)*(pnorm(c0)-pnorm(b))))*A.loc^2 - 1
Var.sc <- h2*A.sc^2
-
+
return(Var.loc+Var.sc)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Ha3rlsGetmse <- function(abc0, r, MAX){
a <- abc0[1]; b <- abc0[2]; c0 <- abc0[3]
-
+
#constraints
if(a < 0 || a > b || b > c0) return(MAX)
Var <- .Ha3rlsGetvar(a = a, b = b, c0 = c0)
-
+
x <- seq(from=0, to=c0, by=0.01)
bias <- sapply(x, .Ha3rlsGetbias, a=a, b=b, c0=c0)
index <- which.max(bias)
@@ -76,12 +76,12 @@
else
b1 <- optimize(f=.Ha3rlsGetbias, lower=x[index-1], upper=x[index+1],
maximum=TRUE, tol=1e-8, a=a, b=b, c0=c0)$objective
-
+
return(Var + r^2*b1^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Ha3 <- function(r, a.start = 0.25, b.start = 2.5, c.start = 5.0,
delta = 1e-6, MAX = 100){
@@ -117,7 +117,7 @@
body(fct2) <- substitute({ Ind1 <- (abs(x) < a); Ind2 <- (abs(x) < b); Ind3 <- (abs(x) < c0)
A.sc*(A.loc*x*(x*Ind1 + a*sign(x)*(Ind2-Ind1) + a*(c0-abs(x))/(c0-b)*sign(x)*(Ind3-Ind2)) - 1) },
list(A.loc = A.loc, A.sc = A.sc, a = a, b = b, c0 = c0))
-
+
return(IC(name = "IC of Ha3 type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$value, asBias = b1, asCov = res$value - r^2*b1^2),
Modified: pkg/RobLox/R/rlsOptIC_Ha4.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Ha4.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Ha4.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.Ha4rlsGetbias <- function(x, k, a, b, c0){
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
@@ -11,7 +11,7 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Ha4rlsGetvar <- function(a, b, c0, k){
h1 <- 2*(-a*dnorm(a) + pnorm(a) - 0.5 + a^2*(pnorm(b)-pnorm(a))
@@ -24,13 +24,13 @@
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
E.psi.4 <- 3*(2*pnorm(k)-1) - 2*(k^3+3*k)*dnorm(k) + 2*k^4*pnorm(-k)
Var.sc <- (E.psi.4 - beta.k^2)/(2*(2*pnorm(k) - 1) - 4*k*dnorm(k))^2
-
+
return(Var.loc+Var.sc)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Ha4rlsGetmse <- function(abc0k, r, MAX){
a <- abc0k[1]; b <- abc0k[2]; c0 <- abc0k[3]; k <- abc0k[4]
@@ -47,12 +47,12 @@
.Ha4rlsGetbias(x = c0, k = k, a = a, b = b, c0 = c0),
.Ha4rlsGetbias(x = sqrt(beta.k), k = k, a = a, b = b, c0 = c0),
.Ha4rlsGetbias(x = k, k = k, a = a, b = b, c0 = c0))
-
+
return(Var + r^2*bias^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Ha4 <- function(r, a.start = 0.25, b.start = 2.5, c.start = 5.0,
k.start = 1.0, delta = 1e-6, MAX = 100){
@@ -79,7 +79,7 @@
fct2 <- function(x){ A.sc*(pmin(x^2, k^2) - beta.k) }
body(fct2) <- substitute({ A.sc*(pmin(x^2, k^2) - beta.k) },
list(k = k, beta.k = beta.k, A.sc = A.sc))
-
+
return(IC(name = "IC of Ha4 type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$value, asBias = bias, asCov = res$value - r^2*bias^2),
Modified: pkg/RobLox/R/rlsOptIC_HaMad.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_HaMad.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_HaMad.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.HaMadrlsGetvar <- function(a, b, c0){
h1 <- 2*(-a*dnorm(a) + pnorm(a) - 0.5 + a^2*(pnorm(b)-pnorm(a))
@@ -9,12 +9,12 @@
A.loc <- 1/(2*(pnorm(a) - 0.5 - a/(c0-b)*(pnorm(c0)-pnorm(b))))
a.mad <- qnorm(0.75)
b.mad <- 1/(4*a.mad*dnorm(a.mad))
-
+
return(h1*A.loc^2 + b.mad^2)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.HaMadrlsGetmse <- function(abc0, r, MAX){
a <- abc0[1]; b <- abc0[2]; c0 <- abc0[3]
@@ -26,14 +26,14 @@
a.mad <- qnorm(0.75)
b.mad <- 1/(4*a.mad*dnorm(a.mad))
bias.2 <- a^2*A.loc^2 + b.mad^2
-
+
Var <- .HaMadrlsGetvar(a = a, b = b, c0 = c0)
-
+
return(Var + r^2*bias.2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.HaMad <- function(r, a.start = 0.25, b.start = 2.5, c.start = 5.0,
delta = 1e-6, MAX = 100){
@@ -62,5 +62,5 @@
"rlsOptIC.HaMad", paste("where a =", round(a, 3), ", b =", round(b, 3),
", c =", round(c0, 3))),
ncol=2, byrow = TRUE, dimnames=list(character(0), c("method", "message"))),
- CallL2Fam = call("NormLocationScaleFamily")))
+ CallL2Fam = call("NormLocationScaleFamily")))
}
Modified: pkg/RobLox/R/rlsOptIC_Hu1.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Hu1.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Hu1.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,13 +1,13 @@
###############################################################################
-# optimal psi function
+## optimal psi function
###############################################################################
lsHu1.chi <- function(x, k){
- beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
+ beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
return(pmin(x^2, k^2) - beta.k)
}
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.Hu1rlsGetbias <- function(x, k){
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
@@ -17,7 +17,7 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Hu1rlsGetvar <- function(k){
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
@@ -31,7 +31,7 @@
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Hu1rlsGetmse <- function(k, r){
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
@@ -48,12 +48,12 @@
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Hu1 <- function(r, kUp = 2.5, delta = 1e-6){
res <- optimize(f = .Hu1rlsGetmse, lower = 1e-4, upper = kUp,
tol = delta, r = r)
-
+
k <- res$minimum
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
A.loc <- 1/(2*pnorm(k)-1)
@@ -70,7 +70,7 @@
fct2 <- function(x){ A.sc*(pmin(x^2, k^2) - beta.k) }
body(fct2) <- substitute({ A.sc*(pmin(x^2, k^2) - beta.k) },
list(k = k, beta.k = beta.k, A.sc = A.sc))
-
+
return(IC(name = "IC of Hu1 type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$objective, asBias = bias, asCov = res$objective - r^2*bias^2),
Modified: pkg/RobLox/R/rlsOptIC_Hu2.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Hu2.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Hu2.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# chi function
+## chi function
###############################################################################
lsHu2.chi <- function(x, c0){
beta.c0 <- 2*pnorm(c0) - 1 - 2*c0*dnorm(c0) + 2*c0^2*pnorm(-c0)
@@ -7,7 +7,7 @@
}
###############################################################################
-# Computation of bias
+## Computation of bias
###############################################################################
.Hu2rlsGetbias <- function(x, k, c0){
beta.c0 <- 2*pnorm(c0) - 1 - 2*c0*dnorm(c0) + 2*c0^2*pnorm(-c0)
@@ -17,46 +17,46 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Hu2rlsGetvar <- function(k, c0){
beta.c0 <- 2*pnorm(c0) - 1 - 2*c0*dnorm(c0) + 2*c0^2*pnorm(-c0)
Var.loc <- (2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*(1-pnorm(k)))/(2*pnorm(k)-1)^2
-
+
hilf <- 3*(2*pnorm(c0)-1) - 2*c0^3*dnorm(c0) - 6*c0*dnorm(c0) + 2*c0^4*pnorm(-c0)
Var.sc <- (hilf - beta.c0^2)/(2*(2*pnorm(c0) - 1) - 4*c0*dnorm(c0))^2
-
+
return(Var.loc + Var.sc)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Hu2rlsGetmse <- function(kc0, r, MAX){
k <- kc0[1]; c0 <- kc0[2]
-
+
# constraints for k, c0
if(k < 0 || c0 < 0) return(MAX)
beta.c0 <- 2*pnorm(c0) - 1 - 2*c0*dnorm(c0) + 2*c0^2*pnorm(-c0)
A.loc <- 1/(2*pnorm(k)-1)
A.sc <- 1/(2*(2*pnorm(c0) - 1) - 4*c0*dnorm(c0))
-
+
bias <- max(.Hu2rlsGetbias(x = 0, k = k, c0 = c0),
.Hu2rlsGetbias(x = k, k = k, c0 = c0),
.Hu2rlsGetbias(x = c0, k = k, c0 = c0),
.Hu2rlsGetbias(x = sqrt(beta.c0), k = k, c0 = c0),
.Hu2rlsGetbias(x = sqrt(max(0, beta.c0-0.5*A.loc^2/A.sc^2)),
k = k, c0 = c0))
-
+
Var <- .Hu2rlsGetvar(k = k, c0 = c0)
-
+
return(Var + r^2*bias^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Hu2 <- function(r, k.start = 1.5, c.start = 1.5, delta = 1e-6, MAX = 100){
res <- optim(c(k.start, c.start), .Hu2rlsGetmse, method = "Nelder-Mead",
@@ -81,7 +81,7 @@
fct2 <- function(x){ A.sc*(pmin(x^2, c0^2) - beta.c0) }
body(fct2) <- substitute({ A.sc*(pmin(x^2, c0^2) - beta.c0) },
list(c0 = c0, beta.c0 = beta.c0, A = A.sc))
-
+
return(IC(name = "IC of Hu2 type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$value, asBias = bias, asCov = res$value - r^2*bias^2),
Modified: pkg/RobLox/R/rlsOptIC_Hu2a.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Hu2a.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Hu2a.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,12 +1,12 @@
###############################################################################
-# psi function
+## psi function
###############################################################################
.Hu2arlsGetpsi <- function(x, k1, k2){
return(sign(x)*pmax(k1,pmin(abs(x),k2)))
}
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.Hu2arlsGetbias <- function(x, k1, k2){
beta.k12 <- (k1^2*(2*pnorm(k1) - 1) - 2*(k2*dnorm(k2)-k1*dnorm(k1))
@@ -18,7 +18,7 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Hu2arlsGetvar <- function(k1, k2){
beta.k12 <- (k1^2*(2*pnorm(k1) - 1) - 2*(k2*dnorm(k2) - k1*dnorm(k1))
@@ -34,7 +34,7 @@
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Hu2arlsGetmse <- function(k12, r, MAX){
k1 <- k12[1]; k2 <- k12[2]
@@ -60,7 +60,7 @@
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Hu2a <- function(r, k1.start = 0.25, k2.start = 2.5, delta = 1e-6, MAX = 100){
res <- optim(c(k1.start, k2.start), .Hu2arlsGetmse, method = "Nelder-Mead",
@@ -86,7 +86,7 @@
fct2 <- function(x){ A.sc*(pmax(k1,pmin(abs(x),k2))^2 - beta.k12) }
body(fct2) <- substitute({ A.sc*(pmax(k1,pmin(abs(x),k2))^2 - beta.k12) },
list(k1 = k1, k2 = k2, beta.k12 = beta.k12, A.sc = A.sc))
-
+
return(IC(name = "IC of Hu2a type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$value, asBias = bias, asCov = res$value - r^2*bias^2),
Modified: pkg/RobLox/R/rlsOptIC_Hu3.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Hu3.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Hu3.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# chi function
+## chi function
###############################################################################
.Hu3rlsGetchi <- function(x, c1, c2){
beta.c12 <- (c1^2*(2*pnorm(c1) - 1) - 2*(c2*dnorm(c2)-c1*dnorm(c1))
@@ -9,7 +9,7 @@
}
###############################################################################
-# Computation of bias
+## Computation of bias
###############################################################################
.Hu3rlsGetbias <- function(x, k, c1, c2){
return(sqrt(pmin(k^2, x^2)/(2*pnorm(k)-1)^2 + .Hu3rlsGetchi(x = x,c1 = c1, c2 = c2)^2/
@@ -18,30 +18,30 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Hu3rlsGetvar <- function(k, c1, c2){
beta.c12 <- (c1^2*(2*pnorm(c1) - 1) - 2*(c2*dnorm(c2) - c1*dnorm(c1))
+ 2*(pnorm(c2)-pnorm(c1)) + 2*c2^2*pnorm(-c2))
Var.loc <- (2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*(1-pnorm(k)))/(2*pnorm(k)-1)^2
-
+
hilf <- (c1^4*(2*pnorm(c1)-1) - 2*(c2^3*dnorm(c2)-c1^3*dnorm(c1))
- 6*(c2*dnorm(c2)-c1*dnorm(c1)) + 6*(pnorm(c2)-pnorm(c1))
+ 2*c2^4*pnorm(-c2))
Var.sc <- (hilf - beta.c12^2)/(4*c1*dnorm(c1)-4*c2*dnorm(c2)+4*(pnorm(c2)-pnorm(c1)))^2
-
+
return(Var.loc + Var.sc)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Hu3rlsGetmse <- function(kc12, r, MAX){
k <- kc12[1]; c1 <- kc12[2]; c2 <- kc12[3]
#constraints
if(k < 0 || c1 < 0 || c2 < 0 || c2 <= c1) return(MAX)
-
+
beta.c12 <- (c1^2*(2*pnorm(c1) - 1) - 2*(c2*dnorm(c2) - c1*dnorm(c1))
+ 2*(pnorm(c2)-pnorm(c1)) + 2*c2^2*pnorm(-c2))
A.loc <- 1/(2*pnorm(k)-1)
@@ -54,15 +54,15 @@
.Hu3rlsGetbias(x = sqrt(beta.c12), k = k, c1 = c1, c2 = c2),
.Hu3rlsGetbias(x = sqrt(max(0, beta.c12-0.5*A.loc^2/A.sc^2)),
k = k, c1 = c1, c2 = c2))
-
+
Var <- .Hu3rlsGetvar(k = k, c1 = c1, c2 = c2)
-
+
return(Var + r^2*bias^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Hu3 <- function(r, k.start = 1.0, c1.start = 0.1, c2.start = 0.5,
delta = 1e-6, MAX = 100){
@@ -92,7 +92,7 @@
body(fct2) <- substitute({ Ind1 <- (abs(x)<c1); Ind2 <- (abs(x)>c2)
A.sc*(c1^2*Ind1 + x^2*(1-Ind1)*(1-Ind2) + c2^2*Ind2 - beta.c12) },
list(c1 = c1, c2 = c2, A.sc = A.sc))
-
+
return(IC(name = "IC of Hu3 type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$value, asBias = bias, asCov = res$value - r^2*bias^2),
Modified: pkg/RobLox/R/rlsOptIC_HuMad.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_HuMad.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_HuMad.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,31 +1,31 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.HuMadrlsGetvar <- function(k){
Var.loc <- (2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k))/(2*pnorm(k)-1)^2
a.mad <- qnorm(0.75)
Var.sc <- 1/(4*a.mad*dnorm(a.mad))^2
-
+
return(Var.loc+Var.sc)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.HuMadrlsGetmse <- function(k, r){
a.mad <- qnorm(0.75)
b.mad <- 1/(4*a.mad*dnorm(a.mad))
bias <- sqrt(k^2/(2*pnorm(k)-1)^2 + b.mad^2)
-
+
Var <- .HuMadrlsGetvar(k = k)
-
+
return(Var + r^2*bias^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.HuMad <- function(r, kUp = 2.5, delta = 1e-6){
res <- optimize(f = .HuMadrlsGetmse, lower = 1e-4, upper = kUp,
@@ -34,7 +34,7 @@
a.mad <- qnorm(0.75)
b.mad <- 1/(4*a.mad*dnorm(a.mad))
bias <- sqrt(k^2/(2*pnorm(k)-1)^2 + b.mad^2)
-
+
A.loc <- 1/(2*pnorm(k)-1)
fct1 <- function(x){ A.loc*sign(x)*pmin(abs(x), k) }
body(fct1) <- substitute({ A.loc*sign(x)*pmin(abs(x), k) },
@@ -42,7 +42,7 @@
fct2 <- function(x){ b.mad*sign(abs(x) - a.mad) }
body(fct2) <- substitute({ b.mad*sign(abs(x) - a.mad) },
list(a.mad = a.mad, b.mad = b.mad))
-
+
return(IC(name = "IC of HuMad type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$objective, asBias = bias, asCov = res$objective - r^2*bias^2),
Modified: pkg/RobLox/R/rlsOptIC_M.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_M.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_M.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# Computation of C_1
+## Computation of C_1
###############################################################################
.MrlsGetC1 <- function(gg){
integrandC1 <- function(x, gg){ dnorm(x)/(1+gg^2*x^2) }
@@ -9,7 +9,7 @@
###############################################################################
-# Computation of C_3
+## Computation of C_3
###############################################################################
.MrlsGetC3 <- function(gg){
integrandC3 <- function(x, gg){ dnorm(x)*x^4/(1+gg^2*x^2) }
@@ -21,20 +21,20 @@
###############################################################################
-# weight function
+## weight function
###############################################################################
.MrlsGetw <- function(x, gg, b, a1, a3){
gvct <- (a1*x + a3*x^3)/sqrt(1 + gg^2*x^2)
bvct <- sqrt(b^2 - gg^2/(1 + gg^2*x^2))
absgvct <- abs(gvct)
ind1 <- (absgvct < bvct)
-
+
return(ind1 + (1-ind1)*bvct/absgvct)
}
###############################################################################
-# psi function
+## psi function
###############################################################################
.MrlsGetpsi <- function(x, gg, b, a1, a3){
gvct <- (a1*x + a3*x^3)/sqrt(1+gg^2*x^2)
@@ -45,7 +45,7 @@
###############################################################################
-# computation of r
+## computation of r
###############################################################################
.MrlsGetr <- function(b, r, gg, a1, a3){
integrandr <- function(x, gg, b, a1, a3){
@@ -65,19 +65,19 @@
###############################################################################
-# computation of a1 und a3
+## computation of a1 und a3
###############################################################################
.MrlsGeta1a3 <- function(gg, b, a1, a3){
C1 <- .MrlsGetC1(gg = gg)
C3 <- .MrlsGetC3(gg = gg)
-
+
integrand1 <- function(x, gg, b, a1, a3){
dnorm(x)*x^2*.MrlsGetw(x, gg, b, a1, a3)/(1+gg^2*x^2)
}
Int1 <- 2*integrate(integrand1, lower = 0, upper = Inf,
rel.tol = .Machine$double.eps^0.5, gg = gg, b = b,
a1 = a1, a3 = a3)$value
-
+
integrand2 <- function(x, gg, b, a1, a3){
dnorm(x)*x^4*.MrlsGetw(x, gg, b, a1, a3)/(1+gg^2*x^2)
}
@@ -91,17 +91,17 @@
Int3 <- 2*integrate(integrand3, lower = 0, upper = Inf,
rel.tol = .Machine$double.eps^0.5, gg = gg, b = b,
a1 = a1, a3 = a3)$value
-
+
D1 <- Int1*Int3 - Int2^2
a1 <- (Int3*C1 + Int2*C3)/D1
a3 <- (Int2*C1 + Int1*C3)/D1
-
+
return(list(a1 = a1, a3 = a3))
}
###############################################################################
-# computation of b, a1 and a3
+## computation of b, a1 and a3
###############################################################################
.MrlsGetba1a3 <- function(r, gg, a1.start, a3.start, bUp, delta, itmax, check){
a1 <- a1.start; a3 <- a3.start
@@ -110,7 +110,7 @@
a1 = a1, a3 = a3)$root, silent = TRUE)
if(!is.numeric(b)) b <- gg
- iter <- 0
+ iter <- 0
repeat{
iter <- iter + 1
if(iter > itmax){
@@ -120,21 +120,20 @@
return(c(NA, NA, NA))
}
a1.old <- a1; a3.old <- a3; b.old <- b
-
+
a1a3 <- .MrlsGeta1a3(gg = gg, b = b, a1 = a1, a3 = a3)
a1 <- a1a3$a1; a3 <- a1a3$a3
-
+
b <- try(uniroot(.MrlsGetr, lower = gg, upper = bUp,
tol = .Machine$double.eps^0.5, r =r, gg = gg,
a1 = a1, a3 = a3)$root, silent = TRUE)
if(!is.numeric(b)) b <- gg
-
+
if(max(abs(b-b.old), abs(a1-a1.old), abs(a3-a3.old)) < delta)
- break
+ break
}
- if(check) # check constraints
- {
+ if(check){ # check constraints
integrand1 <- function(x, gg, b, a1, a3){
x*.MrlsGetpsi(x, gg, b, a1, a3)*dnorm(x)
}
@@ -160,23 +159,23 @@
###############################################################################
-# computation of asymptotic variance \Ew|\rho|^2
+## computation of asymptotic variance \Ew|\rho|^2
###############################################################################
-.MrlsGetvar <- function(gg, b, a1, a3){
+.MrlsGetvar <- function(gg, b, a1, a3){
C1 <- .MrlsGetC1(gg = gg)
-
+
integrandVar <- function(x, gg, b, a1, a3){
dnorm(x)*(a1*x + a3*x^3)^2*.MrlsGetw(x, gg, b, a1, a3)^2/(1+gg^2*x^2)
}
Int <- integrate(integrandVar, lower = 0, upper = Inf,
rel.tol = .Machine$double.eps^0.5, gg = gg, b = b,
a1 = a1, a3 = a3)$value
-
+
return(2*Int + gg^2*C1)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.MrlsGetmse <- function(gg, r, a1.start, a3.start, bUp, delta, itmax, check){
ba1a3 <- .MrlsGetba1a3(r = r, gg = gg, a1.start = a1.start, a3.start = a3.start,
@@ -185,14 +184,14 @@
Var <- .MrlsGetvar(gg = gg, b = b, a1 = a1, a3 = a3)
mse <- Var + r^2*b^2
-
+
# cat("current gamma:\t", gg, "\tcurrent MSE:\t", mse, "\n")
return(mse)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.M <- function(r, ggLo = 0.5, ggUp = 1.5, a1.start = 0.75, a3.start = 0.25,
bUp = 1000, delta = 1e-5, itmax = 100, check = FALSE){
@@ -206,7 +205,7 @@
a1 <- ba1a3$a1; a3 <- ba1a3$a3; b <- ba1a3$b
gg <- res$minimum
-
+
w <- .MrlsGetw
fct1 <- function(x){
wfct <- w
Modified: pkg/RobLox/R/rlsOptIC_MM2.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_MM2.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_MM2.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# psi function
+## psi function
###############################################################################
.MM2rlsGetpsi <- function(x, c0){
Ind1 <- (abs(x) <= 0.8*c0); Ind2 <- (abs(x) <= 1.0*c0)
@@ -9,7 +9,7 @@
}
###############################################################################
-# chi function
+## chi function
###############################################################################
.MM2rlsGetchi <- function(x, d){
Ind1 <- (abs(x) <= d)
@@ -17,7 +17,7 @@
}
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.MM2rlsGetbias <- function(x, c0, d, A.loc, beta.d, A.sc){
return(sqrt(A.loc^2*.MM2rlsGetpsi(x = x, c0 = c0)^2
@@ -26,7 +26,7 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.MM2rlsGetvar <- function(c0, d){
A.loc <- 1/(2*integrate(f = function(x, c0){ x*.MM2rlsGetpsi(x = x, c0 = c0)*dnorm(x) },
@@ -45,16 +45,16 @@
Var.sc <- A.sc^2*(2*integrate(f = function(x, d){ .MM2rlsGetchi(x = x, d = d)^2*dnorm(x) },
lower = 0, upper = Inf, rel.tol = .Machine$double.eps^0.5,
d = d)$value - beta.d^2)
-
+
return(Var.loc + Var.sc)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.MM2rlsGetmse <- function(c0d, r, MAX){
c0 <- c0d[1]; d <- c0d[2]
-
+
# constraints
if(c0 < 0 || d < 0) return(MAX)
@@ -87,12 +87,12 @@
tol = .Machine$double.eps^0.5, c0=c0, d=d, A.loc = A.loc,
beta.d = beta.d, A.sc = A.sc)$objective
}
-
+
return(.MM2rlsGetvar(c0 = c0, d = d) + r^2*b1^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.MM2 <- function(r, c.start = 1.5, d.start = 2.0, delta = 1e-6, MAX = 100){
res <- optim(c(c.start, d.start), .MM2rlsGetmse, method = "Nelder-Mead",
@@ -112,8 +112,8 @@
x <- seq(from=0, to=max(c0,d), by = 0.01)
bias <- sapply(x, .MM2rlsGetbias, c0 = c0, d = d, A.loc = A.loc,
- beta.d = beta.d, A.sc = A.sc)
- index <- which.max(bias)
+ beta.d = beta.d, A.sc = A.sc)
+ index <- which.max(bias)
if(index==length(x))
b1 <- optimize(f=.MM2rlsGetbias, lower=x[index-1], upper=x[index], maximum=TRUE,
Modified: pkg/RobLox/R/rlsOptIC_Tu1.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Tu1.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Tu1.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,18 +1,18 @@
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.Tu1rlsGetbias <- function(x, a){
A.loc <- 1/(2*a*dnorm(a)*(a^2 - 15) + (a^4 - 6*a^2 + 15)*(2*pnorm(a)-1))
A.sc <- 1/(2*A.loc*integrate(f = function(x, a0){ 2*x^2*(a0^2-x^2)*(a0^2-3*x^2)*dnorm(x) },
lower = 0, upper = a, rel.tol = .Machine$double.eps^0.5, a0 = a)$value)
-
+
return(sqrt(x^2*(pmax((a^2-x^2),0))^4*A.loc^2
+ (x^2*pmax((a^2-x^2),0)^2*A.loc - 1)^2*A.sc^2))
}
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Tu1rlsGetvar <- function(a){
A.loc <- 1/(2*a*dnorm(a)*(a^2 - 15) + (a^4 - 6*a^2 + 15)*(2*pnorm(a)-1))
@@ -28,7 +28,7 @@
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Tu1rlsGetmse <- function(a, r){
x <- seq(from = 0, to = a, by = 0.01)
@@ -40,12 +40,12 @@
else
b <- optimize(f=.Tu1rlsGetbias, lower=x[index-1], upper=x[index+1],
maximum=TRUE, tol=.Machine$double.eps^0.5, a=a)$objective
-
+
return(.Tu1rlsGetvar(a = a) + r^2*b^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Tu1 <- function(r, aUp = 10, delta = 1e-6){
res <- optimize(f = .Tu1rlsGetmse, lower = 1e-4, upper = aUp,
@@ -57,8 +57,8 @@
lower = 0, upper = a, rel.tol = .Machine$double.eps^0.5, a0 = a)$value)
x <- seq(from = 0, to = a, by = 0.01)
- bias <- sapply(x, .Tu1rlsGetbias, a = a)
- index <- which.max(bias)
+ bias <- sapply(x, .Tu1rlsGetbias, a = a)
+ index <- which.max(bias)
if(index==length(x))
b <- optimize(f=.Tu1rlsGetbias, lower=x[index-1], upper=x[index],
maximum=TRUE, tol=.Machine$double.eps^0.5, a=a)$objective
@@ -72,7 +72,7 @@
fct2 <- function(x){ A.sc*(A.loc*x^2*(a^2 - x^2)^2*(abs(x) < a) - 1) }
body(fct2) <- substitute({ A.sc*(A.loc*x^2*(a^2 - x^2)^2*(abs(x) < a) - 1) },
list(a = a, A.loc = A.loc, A.sc = A.sc))
-
+
return(IC(name = "IC of Tu1 type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$objective, asBias = b, asCov = res$objective - r^2*b^2),
Modified: pkg/RobLox/R/rlsOptIC_Tu2.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_Tu2.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_Tu2.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# computation of bias
+## computation of bias
###############################################################################
.Tu2rlsGetbias <- function(x, a, k){
beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
@@ -9,7 +9,7 @@
}
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.Tu2rlsGetvar <- function(a, k){
h1 <- 2*integrate(f = function(x, a0){ x^2*(a0^2-x^2)^4*dnorm(x) }, lower = 0,
@@ -24,26 +24,26 @@
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.Tu2rlsGetmse <- function(ak, r, MAX){
a <- ak[1]; k <- ak[2]
-
+
# constraints
if(a < 0 || k < 0) return(MAX)
- beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
+ beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
b <- max(.Tu2rlsGetbias(x = 0, a = a, k = k),
.Tu2rlsGetbias(x = a, a = a, k = k),
.Tu2rlsGetbias(x = a/sqrt(5), a = a, k = k),
.Tu2rlsGetbias(x = k, a = a, k = k),
.Tu2rlsGetbias(x = sqrt(beta.k), a = a, k = k))
-
+
return(.Tu2rlsGetvar(a = a, k = k) + r^2*b^2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.Tu2 <- function(r, a.start = 5, k.start = 1.5, delta = 1e-6, MAX = 100){
res <- optim(c(a.start, k.start), .Tu2rlsGetmse, method = "Nelder-Mead",
@@ -52,13 +52,13 @@
a <- res$par[1]; k <- res$par[2]
A.loc <- 1/(2*a*dnorm(a)*(a^2 - 15) + (a^4 - 6*a^2 + 15)*(2*pnorm(a)-1))
- beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
+ beta.k <- 2*pnorm(k) - 1 - 2*k*dnorm(k) + 2*k^2*pnorm(-k)
bias <- max(.Tu2rlsGetbias(x = 0, a = a, k = k),
.Tu2rlsGetbias(x = a, a = a, k = k),
.Tu2rlsGetbias(x = a/sqrt(5), a = a, k = k),
.Tu2rlsGetbias(x = k, a = a, k = k),
.Tu2rlsGetbias(x = sqrt(beta.k), a = a, k = k))
-
+
A.sc <- 1/(2*(2*pnorm(k) - 1) - 4*k*dnorm(k))
fct1 <- function(x){ A.loc*x*(a^2 - x^2)^2*(abs(x) < a) }
body(fct1) <- substitute({ A.loc*x*(a^2 - x^2)^2*(abs(x) < a) },
Modified: pkg/RobLox/R/rlsOptIC_TuMad.R
===================================================================
--- pkg/RobLox/R/rlsOptIC_TuMad.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rlsOptIC_TuMad.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,5 +1,5 @@
###############################################################################
-# computation of asymptotic variance
+## computation of asymptotic variance
###############################################################################
.TuMadrlsGetvar <- function(a, r){
A.loc <- 1/(2*a*dnorm(a)*(a^2 - 15) + (a^4 - 6*a^2 + 15)*(2*pnorm(a)-1))
@@ -7,12 +7,12 @@
upper = a, rel.tol = .Machine$double.eps^0.5, a0 = a)$value
a.mad <- qnorm(0.75)
b.mad <- 1/(4*a.mad*dnorm(a.mad))
-
+
return(A.loc^2*h1 + b.mad^2)
}
###############################################################################
-# computation of maximum asymptotic MSE
+## computation of maximum asymptotic MSE
###############################################################################
.TuMadrlsGetmse <- function(a, r){
A.loc <- 1/(2*a*dnorm(a)*(a^2 - 15) + (a^4 - 6*a^2 + 15)*(2*pnorm(a)-1))
@@ -20,12 +20,12 @@
b.mad <- 1/(4*a.mad*dnorm(a.mad))
b.2 <- sqrt(A.loc^2*(16/25/sqrt(5)*a^5)^2 + b.mad^2)
-
+
return(.TuMadrlsGetvar(a = a) + r^2*b.2)
}
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rlsOptIC.TuMad <- function(r, aUp = 10, delta = 1e-6){
res <- optimize(f = .TuMadrlsGetmse, lower = 1e-4, upper = aUp,
@@ -44,12 +44,12 @@
fct2 <- function(x){ b.mad*sign(abs(x) - a.mad) }
body(fct2) <- substitute({ b.mad*sign(abs(x) - a.mad) },
list(a.mad = a.mad, b.mad = b.mad))
-
+
return(IC(name = "IC of TuMad type",
Curve = EuclRandVarList(RealRandVariable(Map = list(fct1, fct2), Domain = Reals())),
Risks = list(asMSE = res$objective, asBias = bias, asCov = res$objective - r^2*bias^2),
Infos = matrix(c("rlsOptIC.TuMad", "optimally robust IC for TuMad estimators and 'asMSE'",
"rlsOptIC.TuMad", paste("where a =", round(a, 3))),
ncol=2, byrow = TRUE, dimnames=list(character(0), c("method", "message"))),
- CallL2Fam = call("NormLocationScaleFamily")))
+ CallL2Fam = call("NormLocationScaleFamily")))
}
Modified: pkg/RobLox/R/roblox.R
===================================================================
--- pkg/RobLox/R/roblox.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/roblox.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,14 +1,14 @@
###############################################################################
-# computation of radius-minimax IC
+## computation of radius-minimax IC
###############################################################################
.getlsInterval <- function(r, rlo, rup, mean, sd, delta, A.loc.start,
- a.sc.start, A.sc.start, bUp, itmax){
+ a.sc.start, A.sc.start, bUp, itmax){
Ab <- rlsOptIC.AL(r = r, mean = mean, sd = sd, A.loc.start = A.loc.start,
a.sc.start = a.sc.start, A.sc.start = A.sc.start,
bUp = bUp, delta = delta, itmax = itmax, computeIC = FALSE)
if(rlo == 0){
- efflo <- (sum(diag(Ab$A)) - Ab$b^2*r^2)/(1.5*sd^2)
+ efflo <- (sum(diag(Ab$A)) - Ab$b^2*r^2)/(1.5*sd^2)
}else{
Ablo <- rlsOptIC.AL(r = rlo, mean = mean, sd = sd, A.loc.start = A.loc.start,
a.sc.start = a.sc.start, A.sc.start = A.sc.start,
@@ -20,10 +20,10 @@
a.sc.start = a.sc.start, A.sc.start = A.sc.start,
bUp = bUp, delta = delta, itmax = itmax, computeIC = FALSE)
effup <- (sum(diag(Ab$A)) - Ab$b^2*(r^2 - rup^2))/sum(diag(Abup$A))
-
+
return(effup-efflo)
}
-.getlInterval <- function(r, rlo, rup, mean, sd, bUp){
+.getlInterval <- function(r, rlo, rup, mean, sd, bUp){
Ab <- rlOptIC(r = r, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
if(rlo == 0){
@@ -35,10 +35,10 @@
Abup <- rlOptIC(r = rup, mean = mean, sd = sd, bUp = bUp, computeIC = FALSE)
effup <- (Ab$A - Ab$b^2*(r^2 - rup^2))/Abup$A
-
+
return(effup-efflo)
}
-.getsInterval <- function(r, rlo, rup, mean, sd, delta, bUp, itmax){
+.getsInterval <- function(r, rlo, rup, mean, sd, delta, bUp, itmax){
Ab <- rsOptIC(r = r, mean = mean, sd = sd, bUp = bUp, delta = delta,
itmax = itmax, computeIC = FALSE)
@@ -53,11 +53,11 @@
Abup <- rsOptIC(r = rup, mean = mean, sd = sd, bUp = bUp, delta = delta,
itmax = itmax, computeIC = FALSE)
effup <- (Ab$A - Ab$b^2*(r^2 - rup^2))/Abup$A
-
+
return(effup-efflo)
}
###############################################################################
-# optimally robust estimator for normal location and/or scale
+## optimally robust estimator for normal location and/or scale
###############################################################################
roblox <- function(x, mean, sd, eps, eps.lower, eps.upper, initial.est = "ksMD",
tol = 1e-6, A.loc.start = 1, a.sc.start = 0, A.sc.start = 0.5,
@@ -86,7 +86,7 @@
}
if((initial.est != "ksMD") && (initial.est != "med"))
stop("invalid 'initial.est'")
-
+
if(missing(mean) & missing(sd)){
if(!is.numeric(A.loc.start) || !is.numeric(a.sc.start) || !is.numeric(A.sc.start))
stop("Starting values 'A.loc.start', 'a.sc.start' and 'A.sc.start' have to be numeric")
@@ -133,7 +133,7 @@
A.sc.start = A.sc.start, bUp = bUp, delta = tol,
itmax = itmax)
if(rlo == 0){
- ineff <- (sum(diag(stand(IC1))) - clip(IC1)^2*r^2)/(1.5*sd^2)
+ ineff <- (sum(diag(stand(IC1))) - clip(IC1)^2*r^2)/(1.5*sd^2)
}else{
Ablo <- rlsOptIC.AL(r = rlo, mean = mean, sd = sd, A.loc.start = A.loc.start,
a.sc.start = a.sc.start, A.sc.start = A.sc.start,
@@ -254,4 +254,4 @@
}
}
}
-}
+}
Modified: pkg/RobLox/R/rsOptIC.R
===================================================================
--- pkg/RobLox/R/rsOptIC.R 2008-02-19 05:51:12 UTC (rev 41)
+++ pkg/RobLox/R/rsOptIC.R 2008-02-19 05:59:47 UTC (rev 42)
@@ -1,7 +1,7 @@
###############################################################################
-# computation of centering constang
+## computation of centering constang
###############################################################################
-.ALrsGetz <- function(z, b){
+.ALrsGetz <- function(z, b){
b1 <- sqrt(pmax(z - b,0))
b2 <- sqrt(b + z)
@@ -10,7 +10,7 @@
}
###############################################################################
-# computation of clipping bound
+## computation of clipping bound
###############################################################################
.ALrsGetr <- function(b, r, z){
b1 <- sqrt(pmax(z - b,0))
@@ -22,7 +22,7 @@
###############################################################################
-# optimal IC
+## optimal IC
###############################################################################
rsOptIC <- function(r, mean = 0, sd = 1, bUp = 1000, delta = 1e-6, itmax = 100,
computeIC = TRUE){
@@ -37,7 +37,7 @@
stop("Algorithm did not converge => increase 'itmax'!")
z.old <- z; b.old <- b
-
+
z <- uniroot(.ALrsGetz, lower = 0, upper = 1, tol = .Machine$double.eps^0.5,
b = b)$root
@@ -55,18 +55,18 @@
aa <- dnorm(b2)*(-b2^3 - 3*b2 + 2*z*b2) - dnorm(b1)*(-b1^3 - 3*b1 + 2*z*b1)
aa <- aa + (pnorm(b2) - pnorm(b1))*(z^2 + 3 - 2*z) + b*(1-z)*(1.5 - pnorm(b2) - pnorm(b1))
aa <- aa + b*(b2*dnorm(b2) + b1*dnorm(b1))
-
+
A1 <- 1/(2*aa)
b <- sd*A1*b
a <- sd*A1*(z - 1)
A <- sd^2*A1
-
+
if(computeIC){
return(generateIC(neighbor = ContNeighborhood(radius = r),
L2Fam = NormScaleFamily(mean = mean, sd = sd),
res = list(A = as.matrix(A), a = a, b = b, d = NULL,
risk = list(asMSE = A, asBias = b, asCov = A - r^2*b^2),
- info = c("rsOptIC", "optimally robust IC for AL estimators and 'asMSE'"))))
+ info = c("rsOptIC", "optimally robust IC for AL estimators and 'asMSE'"))))
}else{
return(list(A = A, a = a, b = b))
}
More information about the Robast-commits
mailing list