## Test disabled
quit()
## test kappa functions
require(robustlmm)
## Some session / library diagnostics for "us" (the package developers):
sessionInfo()
.libPaths()
packageDescription("robustbase")
packageDescription("lme4")
packageDescription("robustlmm")
## --------------------------------------------------------------------
calcKappaTau <- robustlmm:::calcKappaTau
calcTau <- robustlmm:::calcTau
E <- function(fun) {
## sum(na.omit(fun(rfm@pp$ghz)*rfm@pp$ghw))
int <- function(x) fun(x)*dnorm(x)
integrate(int, -Inf, Inf)$value
}
testKappaTau <- function(rho, wExp) {
rho2 <- if (wExp == 2) psi2propII(rho) else rho
kappa <- calcKappaTau(rho2, 1)
fun <- switch(wExp,
function(x) rho@psi(x)*x - kappa*rho@wgt(x),
function(x) rho@psi(x)^2 - kappa*rho@wgt(x)^2)
E(fun)
}
testKappa <- function(rho) {
stopifnot(all.equal(0, testKappaTau(rho, 1), tolerance = 1e-4),
all.equal(0, testKappaTau(rho, 2), tolerance = 1e-4))
}
testKappa(huberPsiRcpp)
testKappa(smoothPsi)
E2 <- function(fun, ...) {
int1 <- function(x, y, ...) fun(x, y, ...)*dnorm(x)*dnorm(y)
int2 <- Vectorize(function(y, ...) integrate(int1, -Inf, Inf, y=y, ...)$value)
integrate(int2, -Inf, Inf, ...)$value
}
rfm <- rlmer(Yield ~ (1|Batch), Dyestuff)
## FIXME make this an argument of the original call
rfm@pp$initGH(50)
testTau <- function(rho, rho.sigma, wExp, a, s) {
psi <- rho@psi
rho.e <- rho.sigma
rho.sigma2 <- if (wExp == 2) psi2propII(rho.sigma) else rho.sigma
kappa <- calcKappaTau(rho.sigma2, 1)
i <- 1
fun <- switch(wExp + 1,
{ ## wExp.e == 0:
wgt <- function(x) ifelse(x == 0, rho.e$Dpsi(0)/2, rho@rho(x)/(x*x))
function(w, v, tau) {
t <- (v-a[i]*psi(v)+w*s[i])/tau
rho.e@rho(t) - kappa*wgt(t)
} }, ## wExp.e == 1:
function(w, v, tau) {
t <- (v-a[i]*psi(v)+w*s[i])/tau
rho.e@psi(t)*t - kappa*rho.e@wgt(t)
}, ## wExp.e == 2:
function(w, v, tau) {
t <- (v-a[i]*psi(v)+w*s[i])/tau
rho.e@psi(t)^2 - kappa*rho.e@wgt(t)^2
})
tau <- calcTau(a, s, rho, rho.sigma2, rfm@pp, kappa)
print(tau)
ret <- tau
for (i in seq_along(a)) ret[i] <- E2(fun, tau = tau[i])
print(ret)
ret
}
a <- c(0.1, 0.4, 1, 1, 0.8)
s <- c(0.4, 0.5, 0.8, 0.1, 0.1)
## FIXME: increase accuracy of calcTau.
stopifnot(all.equal(rep(0, length(a)), testTau(huberPsiRcpp, huberPsiRcpp, 1, a, s), tolerance = 1e-2),
all.equal(rep(0, length(a)), testTau(huberPsiRcpp, huberPsiRcpp, 2, a, s), tolerance = 1e-2))
stopifnot(all.equal(rep(0, length(a)), testTau(smoothPsi, huberPsiRcpp, 1, a, s), tolerance = 1e-2),
all.equal(rep(0, length(a)), testTau(smoothPsi, huberPsiRcpp, 2, a, s), tolerance = 1e-1))
stopifnot(all.equal(rep(0, length(a)), testTau(smoothPsi, smoothPsi, 1, a, s), tolerance = 1e-2),
all.equal(rep(0, length(a)), testTau(smoothPsi, smoothPsi, 2, a, s), tolerance = 1e-1))
sPsi <- chgDefaults(smoothPsi, k = 1, s = 10)
stopifnot(all.equal(rep(0, length(a)), testTau(smoothPsi, sPsi, 1, a, s), tolerance = 1e-2),
all.equal(rep(0, length(a)), testTau(smoothPsi, sPsi, 2, a, s), tolerance = 1e-2))
## to save time in checks, also test init argument here:
initList <- list(fixef=fixef(rfm), u=getME(rfm, "u"),
sigma=sigma(rfm), theta=getME(rfm, "theta"))
system.time(rfm2 <- rlmer(Yield ~ (1|Batch), Dyestuff,
init=initList))
rfm2@call <- rfm@call
rfm2@optinfo <- rfm@optinfo
rfm2@pp$.setTau_e <- rfm@pp$.setTau_e
rfm2@pp$.tau_e <- rfm@pp$.tau_e
rfm2@pp$.setTbk <- rfm@pp$.setTbk
rfm2@pp$.Tbk <- rfm@pp$.Tbk
rfm2@pp$groupsA <- rfm@pp$groupsA
rfm@pp$initGH()
stopifnot(all.equal(rfm, rfm2, tolerance = 1e-6))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.