robeth-package(robeth)
robeth-package()所属R语言包：robeth
Interface for the FORTRAN programs developed at the ETH-Zuerich and then at IUMSP-Lausanne
接口的FORTRAN程序开发的ETH苏黎世和洛桑，在IUMSP，然后
译者：生物统计家园网 机器人LoveR
描述----------Description----------
This package allows the computation of a broad class of procedures based on M-estimation and high breakdown point estimation, including robust regresion, robust testing of linear hypotheses and robust covariances. The reference book quoted below is required for the theoritical background of the statistical and numerical methods
此包允许广泛的一类程序，基于M-估计和高击穿的点估计，包括的强大regresion，功能强大的测试线性假设和强大的协方差的计算。以下引用的参考书所需的背景的统计和数值方法理论
Details
详细信息----------Details----------
</table>
</ TABLE>
（作者）----------Author(s)----------
Alfio Marazzi <Alfio.Marazzi@chuv.ch>
Maintainer: A. Randriamiharisoa <Alex.Randriamiharisoa@chuv.ch>
参考文献----------References----------
Marazzi A., (1993) Algorithm, Routines, and S functions for Robust Statistics. Wadsworth & Brooks/cole, Pacific Grove, California.
实例----------Examples----------
library(robeth)
#[]
# ------------- Examples of Chapter 1: Location problems ------------------------------[-------------第1章：地点问题的例子------------------------------]
#[]
y <- c(6.0,7.0,5.0,10.5,8.5,3.5,6.1,4.0,4.6,4.5,5.9,6.5)
n <- 12
dfvals()
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# M-estimate (tm) of location and confidence interval (tl,tu)[M-估计（TM）的位置和置信区间（TL，涂）]
#[]
dfrpar(as.matrix(y),"huber")
libeth()
s <- lilars(y); t0 <- s$theta; s0 <- s$sigma
s <- lyhalg(y=y,theta=t0,sigmai=s0)
tm <- s$theta; vartm <- s$var
s <- quant(0.975)
tl <- tm-s$x*sqrt(vartm)
tu <- tm+s$x*sqrt(vartm)
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Hodges and Lehmann estimate (th) and confidence interval (zl,zu)[霍奇斯和莱曼估计（TH）和置信区间（ZL，祖冲之）]
# []
m <- n*(n+1)/2 # n even[n为偶数]
k1 <- m/2; k2 <- k1+1
z1 <- lyhdle(y=y,k=k1); z2 <- lyhdle(y=y,k=k2)
th <- (z1$hdle+z2$hdle)/2.
ku <- liindh(0.95,n); kl <- liindh(0.05,n)
zu <- lyhdle(y=y,k=ku$k); zl <- lyhdle(y=y,k=kl$k)
#.......................................................................[.................................................. .....................]
{
cat(" tm, tl, tu : ",round(c(tm,tl,tu),3),"\n")
cat(" th, zl, zu : ",round(c(th,zl$hdle,zu$hdle),3),"\n")
}
# tm, tl, tu : 5.809 4.748 6.87 [TM，TL，TU：5.809 4.748 6.87]
# th, zl, zu : 5.85 5 7 [日，ZL，祖：5.85 5 7]
#=======================================================================[================================================== =====================]
#[]
# Two sample problem[两个样本的问题]
#[]
y <- c(17.9,13.3,10.6,7.6,5.7,5.6,5.4,3.3,3.1,0.9)
n <- 10
x <- c(7.7,5.0,1.7,0.0,-3.0,-3.1,-10.5)
m <- 7
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Estimate (dm) and confidence interval [dl,du] based on Mann-Whitney[估计（DM）和置信区间[DL，杜]的基础上采用Mann-Whitney]
#[]
k1 <- m*n/2; k2 <- k1+1
s1 <- lymnwt(x=x,y=y,k=k1); s2 <- lymnwt(x=x,y=y,k=k2)
dm <- (s1$tmnwt+s2$tmnwt)/2.0
sl <- liindw(0.05,m,n); kl <- sl$k
s <- lymnwt(x=x,y=y,k=kl); dl <- s$tmnwt
s <- lymnwt(x=x,y=y,k=m*n-kl+1); du <- s$tmnwt
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Tau-test . P-value (P)[头测试。 P-值（P）的]
#[]
z <- c(x,y)
dfrpar(as.matrix(z),"huber")
libeth()
s <- lytau2(z=z,m=m,n=n)
P <- s$p
#[]
# estimate (tm) and confidence interval (tl,tu)[估计（TM）和置信区间（TL，TU）]
#[]
tm <- s$deltal
c22<- s$cov[3]
s <- quant(0.975)
tl <- tm-s$x*sqrt(c22)
tu <- tm+s$x*sqrt(c22)
#.......................................................................[.................................................. .....................]
{
cat("dm, dl, du:",round(c(dm,dl,du),3),"\n")
cat(" , tm, tl, tu:",round(c(P,tm,tl,tu),3),"\n")
}
# dm, dl, du: 6.35 2.9 12.9 [DM，DL，杜：6.35 2.9 12.9]
# P, tm, tl, tu: 0.014 6.918 1.562 12.273 [P，TM，TL，TU：0.014 6.918 1.562 12.273]
#[]
# ---- Examples of Chapter 2: M-estimates of coefficients and scale in linear regression ------[----第2章的例子：M-估计的线性回归系数和规模------]
#[]
# Read data; declare the vector wgt; load defaults[读取数据;申报的矢量WGT;载入预设值]
#[]
z <- c(-1, -2, 0, 35, 1, 0, -3, 20,
-1, -2, 0, 30, 1, 0, -3, 39,
-1, -2, 0, 24, 1, 0, -3, 16,
-1, -2, 0, 37, 1, 0, -3, 27,
-1, -2, 0, 28, 1, 0, -3, -12,
-1, -2, 0, 73, 1, 0, -3, 2,
-1, -2, 0, 31, 1, 0, -3, 31,
-1, -2, 0, 21, 1, 0, -1, 26,
-1, -2, 0, -5, 1, 0, -1, 60,
-1, 0, 0, 62, 1, 0, -1, 48,
-1, 0, 0, 67, 1, 0, -1, -8,
-1, 0, 0, 95, 1, 0, -1, 46,
-1, 0, 0, 62, 1, 0, -1, 77,
-1, 0, 0, 54, 1, 0, 1, 57,
-1, 0, 0, 56, 1, 0, 1, 89,
-1, 0, 0, 48, 1, 0, 1, 103,
-1, 0, 0, 70, 1, 0, 1, 129,
-1, 0, 0, 94, 1, 0, 1, 139,
-1, 0, 0, 42, 1, 0, 1, 128,
-1, 2, 0, 116, 1, 0, 1, 89,
-1, 2, 0, 105, 1, 0, 1, 86,
-1, 2, 0, 91, 1, 0, 3, 140,
-1, 2, 0, 94, 1, 0, 3, 133,
-1, 2, 0, 130, 1, 0, 3, 142,
-1, 2, 0, 79, 1, 0, 3, 118,
-1, 2, 0, 120, 1, 0, 3, 137,
-1, 2, 0, 124, 1, 0, 3, 84,
-1, 2, 0, -8, 1, 0, 3, 101)
xx <- matrix(z,ncol=4, byrow=TRUE)
dimnames(xx) <- list(NULL,c("z2","xS","xT","y"))
z2 <- xx[,"z2"]; xS <- xx[,"xS"]; xT <- xx[,"xT"]
x <- cbind(1, z2, xS+xT, xS-xT, xS^2+xT^2, xS^2-xT^2, xT^3)
y <- xx[,"y"]
wgt <- vector("numeric",length(y))
n <- 56; np <- 7
dfvals()
# Set parameters for Huber estimate[胡贝尔估计的参数设置]
dfrpar(x, "huber")
# Compute the constants beta, bet0, epsi2 and epsip[计算常数测试，bet0，epsi2和epsip]
ribeth(wgt)
ribet0(wgt)
s <- liepsh()
epsi2 <- s$epsi2; epsip <- s$epsip
#[]
# Least squares solution (theta0,sigma0)[最小二乘解（theta0，sigma0）]
#[]
z <- riclls(x, y)
theta0<- z$theta; sigma0 <- z$sigma
# Preliminary estimate of the covariance matrix of the coefficients[初步估计的协方差矩阵的系数]
cv <- kiascv(z$xt, fu=epsi2/epsip^2, fb=0.)
cov <- cv$cov
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Solution (theta1,sigma1) by means of RYHALG.[解决方案（THETA1，sigma1）由手段RYHALG。]
#[]
zr <- ryhalg(x,y,theta0,wgt,cov,sigmai=sigma0,ic=0)
theta1<- zr$theta[1:np]; sigma1 <- zr$sigmaf; nit1 <- zr$nit
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Solution (theta2,sigma2) by means of RYWALG (recompute cov) [的溶液（theta2，sigma2）由手段RYWALG（重新计算覆盖）]
#[]
cv <- ktaskv(x, f=epsi2/epsip^2)
zr <- rywalg(x, y, theta0, wgt, cv$cov, sigmai=sigma0)
theta2 <- zr$theta[1:np]; sigma2 <- zr$sigmaf; nit2 <- zr$nit
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Solution (theta3,sigma3) by means of RYSALG with ISIGMA=2.[的的溶液（theta3，sigma3）通过的RYSALG与ISIGMA = 2。]
#[]
zr <- rysalg(x,y, theta0, wgt, cv$cov, sigma0, isigma=2)
theta3 <- zr$theta[1:np]; sigma3 <- zr$sigmaf; nit3 <- zr$nit
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Solution (theta4,sigma4) by means of RYNALG with ICNV=2 and ISIGMA=0.[通过与ICNV = 2和ISIGMA = 0 RYNALG的的溶液（theta4，sigma4）。]
#[]
# Invert cov[反向覆盖]
covm1 <- cv$cov
zc <- mchl(covm1,np)
zc <- minv(zc$a, np)
zc <- mtt1(zc$r,np)
covm1 <- zc$b
zr <- rynalg(x,y, theta0, wgt, covm1, sigmai=sigma3,
iopt=1, isigma=0, icnv=2)
theta4 <- zr$theta[1:np]; sigma4 <- zr$sigmaf; nit4 <- zr$nit
#.......................................................................[.................................................. .....................]
{
cat("theta0 : ",round(theta0[1:np],3),"\n")
cat("sigma0 : ",round(sigma0,3),"\n")
cat("theta1 : ",round(theta1,3),"\n")
cat("sigma1, nit1 : ",round(sigma1,3),nit1,"\n")
cat("theta2 : ",round(theta2,3),"\n")
cat("sigma2, nit2 : ",round(sigma2,3),nit2,"\n")
cat("theta3 : ",round(theta3,3),"\n")
cat("sigma3, nit3 : ",round(sigma3,3),nit3,"\n")
cat("theta4 : ",round(theta4,3),"\n")
cat("sigma4, nit4 : ",round(sigma4,3),nit4,"\n")
}
# theta0 : 68.634 3.634 24.081 -8.053 -0.446 -0.179 -1.634 [theta0：68.634 3.634 24.081 -8.053 -0.446 -0.179 -1.634]
# sigma0 : 26.635 [sigma0：26.635]
# theta1 : 70.006 5.006 24.742 -6.246 -0.079 0.434 -1.487 [THETA1：70.006 5.006 24.742 -6.246 -0.079 0.434 -1.487]
# sigma1, nit1 : 23.564 7 [sigma1，NIT1：23.564]
# theta2 : 70.006 5.006 24.742 -6.245 -0.079 0.434 -1.487 [theta2：70.006 5.006 24.742 -6.245 -0.079 0.434 -1.487]
# sigma2, nit2 : 23.563 7 [sigma2，nit2的：23.563]
# theta3 : 69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48 [theta3：69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48]
# sigma3, nit3 : 22.249 3 [sigma3，nit3：22.249]
# theta4 : 69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48 [theta4：69.993 5.002 24.766 -6.214 -0.055 0.44 -1.48]
# sigma4, nit4 : 22.249 3 [sigma4，nit4：22.249]
#[]
# ---- Examples of Chapter 3: Weights for bounded influence regression ------[----第3章的例子：重量为界的影响力回归------]
#[]
#=======================================================================[================================================== =====================]
rbmost <- function(x,y,cc,usext=userfd) {
n <- nrow(x); np <- ncol(x); dfcomn(xk=np)
.def <- .dFv; .def$itw <- 1; .dFv <<- .def
z <- wimedv(x)
z <- wyfalg(x, z$a, y, exu=usext); nitw <- z$nit
wgt <- 1/z$dist; wgt[wgt>1.e6] <- 1.e6
z <- comval()
bto <- z$bt0; ipso <- z$ipsi; co <- z$c
z <- ribet0(wgt, itype=2, isqw=0)
xt <- x*wgt; yt <- y * wgt
z <- rilars(xt, yt)
theta0 <- z$theta; sigma0 <- z$sigma
rs <- z$rs/wgt; r1 <- rs/sigma0
dfcomn(ipsi=1,c=cc)
z <- liepsh(cc)
den <- z$epsip
g <- Psp(r1)/den # (see Psi in Chpt. 14)[（见幽CHPT。14）]
dfcomn(ipsi=ipso, c=co, bet0=bto)
list(theta=theta0, sigma=sigma0, rs=rs, g=g, nitw=nitw)
}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Mallows-standard estimate (with wyfalg and rywalg)[锦葵标准估计（与wyfalg和rywalg）]
#[]
Mal.Std <- function(x, y, b2=-1, cc=-1, isigma=2) {
n <- length(y); np <- ncol(x)
dfrpar(x, "Mal-Std", b2, cc); .def <- .dFv
if (isigma==1) {dfcomn(d=.dFv$ccc); .def$isg <- 1; .dFv <<- .def}
# Weights [重量]
z <- wimedv(x)
z <- wyfalg(x, z$a, y); nitw <- z$nit
wgt <- Www(z$dist) # See Www in Chpt. 14[请参阅www的CHPT。 14]
# Initial cov. matrix of coefficient estimates[初始覆盖。矩阵的系数估计值]
z <- kiedch(wgt)
cov <- ktaskw(x, z$d, z$e, f=1/n)
# Initial theta and sigma[最初的theta和西格玛]
z <- rbmost(x,y,1.5,userfd)
theta0 <- z$theta; sigma0 <- z$sigma; nitw0 <- z$nitw
# Final theta and sigma[最终的theta和西格玛]
if (.dFv$isg==1) ribeth(wgt) else ribet0(wgt)
z <- rywalg(x, y,theta0,wgt,cov$cov, sigmai=sigma0)
theta1 <- z$theta[1:np]; sigma1 <- z$sigmaf; nit1 <- z$nit
# Covariance matrix of coefficient estimates[系数估计值的协方差矩阵]
z <- kfedcc(wgt, z$rs, sigma=sigma1)
cov <- ktaskw(x, z$d, z$e, f=(sigma1^2)/n)
sd1 <- NULL; for (i in 1:np) { j <- i*(i+1)/2
sd1 <- c(sd1,cov$cov[j]) }
sd1 <- sqrt(sd1)
#.......................................................................[.................................................. .....................]
{
cat("rbmost: theta0 : ",round(theta0[1:np],3),"\n")
cat("rbmost: sigma0, nitw : ",round(sigma0,3),nitw0,"\n")
cat("Mallows-Std: theta1 : ",round(theta1,3),"\n")
cat("Mallows-Std: sd1 : ",round(sd1,3),"\n")
cat("Mallows-Std: sigma1, nitw, nit1 : ",round(sigma1,3),nitw,nit1,"\n")
}
#.......................................................................[.................................................. .....................]
list(theta0=theta0[1:np], sigma0=sigma0, nitw=nitw,
theta1=theta1, sigma1=sigma1, nit1=nit1, sd1=sd1)}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Krasker-Welsch estimate (with wynalg and rynalg)[Krasker韦尔施的估计（与wynalg和rynalg）]
#[]
Kra.Wel <- function(x, y, ckw=-1, isigma=2) {
n <- length(y); np <- ncol(x)
dfrpar(x, "Kra-Wel", ckw); .def <- .dFv
if (isigma==1) {dfcomn(d=.dFv$ccc); .def$isg <- 1; .dFv <<- .def}
# Weights [重量]
z <- wimedv(x)
z <- wynalg(x, z$a); nitw <- z$nit
wgt <- Www(z$dist) # See Www in Chpt. 14[请参阅www的CHPT。 14]
# Initial cov. matrix of coefficient estimates[初始覆盖。矩阵的系数估计值]
z <- kiedch(wgt)
cov <- ktaskw(x, z$d, z$e, f=1/n)
# Initial theta and sigma[最初的theta和西格玛]
z <- rbmost(x, y, cc=1.5)
theta0 <- z$theta; sigma0 <- z$sigma; nitw0 <- z$nitw
# Final theta and sigma[最终的theta和西格玛]
if (.dFv$isg==1) ribeth(wgt) else ribet0(wgt)
z <- rynalg(x, y,theta0,wgt,cov$cov, sigmai=sigma0)
theta2 <- z$theta[1:np]; sigma2 <- z$sigma; nit2 <- z$nit
#[]
# Covariance matrix of coefficient estimates[系数估计值的协方差矩阵]
#[]
z <- kfedcc(wgt, z$rs, sigma=sigma2)
cov <- ktaskw(x, z$d, z$e, f=(sigma2^2)/n)
sd2 <- NULL; for (i in 1:np) { j <- i*(i+1)/2
sd2 <- c(sd2,cov$cov[j]) }
sd2 <- sqrt(sd2)
#.......................................................................[.................................................. .....................]
{
cat("rbmost: theta0 : ",round(theta0[1:np],3),"\n")
cat("rbmost: sigma0, nitw : ",round(sigma0,3),nitw0,"\n")
cat("Krasker-Welsch: theta2 : ",round(theta2,3),"\n")
cat("Krasker-Welsch: sd2 : ",round(sd2,3),"\n")
cat("Krasker-Welsch: sigma2, nitw, nit2 : ",round(sigma2,3),nitw,nit2,"\n")
}
#.......................................................................[.................................................. .....................]
list(theta0=theta0[1:np], sigma0=sigma0, nitw=nitw,
theta2=theta2, sigma2=sigma2, nit2=nit2, sd2=sd2)}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Read data; load defaults[读取数据载入预设值]
# []
z <- c( 8.2, 4, 23.005, 1, 7.6, 5, 23.873, 1,
4.6, 0, 26.417, 1, 4.3, 1, 24.868, 1,
5.9, 2, 29.895, 1, 5.0, 3, 24.200, 1,
6.5, 4, 23.215, 1, 8.3, 5, 21.862, 1,
10.1, 0, 22.274, 1, 13.2, 1, 23.830, 1,
12.6, 2, 25.144, 1, 10.4, 3, 22.430, 1,
10.8, 4, 21.785, 1, 13.1, 5, 22.380, 1,
13.3, 0, 23.927, 1, 10.4, 1, 33.443, 1,
10.5, 2, 24.859, 1, 7.7, 3, 22.686, 1,
10.0, 0, 21.789, 1, 12.0, 1, 22.041, 1,
12.1, 4, 21.033, 1, 13.6, 5, 21.005, 1,
15.0, 0, 25.865, 1, 13.5, 1, 26.290, 1,
11.5, 2, 22.932, 1, 12.0, 3, 21.313, 1,
13.0, 4, 20.769, 1, 14.1, 5, 21.393, 1)
x <- matrix(z, ncol=4, byrow=TRUE)
y <- c( 7.6, 7.7, 4.3, 5.9, 5.0, 6.5, 8.3, 8.2, 13.2, 12.6,
10.4, 10.8, 13.1, 12.3, 10.4, 10.5, 7.7, 9.5, 12.0, 12.6,
13.6, 14.1, 13.5, 11.5, 12.0, 13.0, 14.1, 15.1)
dfvals()
dfcomn(xk=4)
cat("Results\n")
z1 <- Mal.Std(x, y)
z2 <- Kra.Wel(x, y)
#[]
# ---- Examples of Chapter 4: Covariance matrix of the coefficient estimates ------[----第4章：系数估计值的协方差矩阵的例子------]
#[]
#[]
# Read x[1:4] and then set x[,4] <- 1[读取X [1:4]，然后设置X [4] < - 1]
# []
z <- c(80, 27, 89, 1, 80, 27, 88, 1, 75, 25, 90, 1,
62, 24, 87, 1, 62, 22, 87, 1, 62, 23, 87, 1,
62, 24, 93, 1, 62, 24, 93, 1, 58, 23, 87, 1,
58, 18, 80, 1, 58, 18, 89, 1, 58, 17, 88, 1,
58, 18, 82, 1, 58, 19, 93, 1, 50, 18, 89, 1,
50, 18, 86, 1, 50, 19, 72, 1, 50, 19, 79, 1,
50, 20, 80, 1, 56, 20, 82, 1, 70, 20, 91, 1)
x <- matrix(z, ncol=4, byrow=TRUE)
n <- 21; np <- 4; ncov <- np*(np+1)/2
dfvals()
# Cov. matrix of Huber-type estimates[冠状病毒。胡贝尔型估计矩阵]
dfrpar(x, "huber")
s <- liepsh()
epsi2 <- s$epsi2; epsip <- s$epsip
z <- rimtrf(x)
xt <- z$x; sg <- z$sg; ip <- z$ip
zc <- kiascv(xt, fu=epsi2/epsip^2, fb=0.)
covi <- zc$cov # Can be used in ryhalg with ic=0[可以用在与集成电路= 0 ryhalg]
zc <- kfascv(xt, covi, f=1, sg=sg, ip=ip)
covf <- zc$cov
#.......................................................................[.................................................. .....................]
str <- rep(" ", ncov); str[cumsum(1:np)] <- "\n"
{
cat("covf:\n")
cat(round(covf,6),sep=str)
}
#[]
# ---- Examples of Chapter 5: Asymptotic relative efficiency ------[----第5章的例子：渐近相对效率------]
#[]
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Huber[胡贝尔]
#[]
dfcomn(ipsi=1,c=1.345,d=1.345)
.dFv$ite <- 1
z <- airef0(mu=3, ialfa=1, sigmx=1)
#.......................................................................[.................................................. .....................]
{
cat(" airef0 : Huber\n reff, alfa, beta, nit: ")
cat(round(c(z$reff,z$alfa,z$beta,z$nit),3),sep=c(", ",", ",", ","\n"))
}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Schweppe: Krasker-Welsch[施韦普：Krasker韦尔施]
#[]
dfcomn(ipsi=1,c=3.76,iucv=3,ckw=3.76,iwww=1)
.dFv$ite <- 3
z <- airef0(mu=3, ialfa=1, sigmx=1)
#.......................................................................[.................................................. .....................]
{
cat(" airef0 : Krasker-Welsch\n reff, alfa, beta, nit: ")
cat(round(c(z$reff,z$alfa,z$beta,z$nit),3),sep=c(", ",", ",", ","\n"))
}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Mallows-Standard[锦葵标准]
#[]
dfcomn(ipsi=1,c=1.5,iucv=1,a2=0,b2=6.16,iww=3)
.dFv$ite <- 2
z <- airef0(mu=3, ialfa=1, sigmx=1)
#.......................................................................[.................................................. .....................]
{
cat(" airef0 : Mallows-Std \n reff, alfa, beta, nit: ")
cat(round(c(z$reff,z$alfa,z$beta,z$nit),3),sep=c(", ",", ",", ","\n"))
}
#=======================================================================[================================================== =====================]
z <- c(1, 0, 0,
1, 0, 0,
1, 0, 0,
1, 0, 0,
0, 1, 0,
0, 1, 0,
0, 1, 0,
0, 1, 0,
0, 0, 1,
0, 0, 1,
0, 0, 1,
0, 0, 1)
tt <- matrix(z,ncol=3,byrow=TRUE)
n <- nrow(tt); mu <- 2
nu <- ncol(tt)
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Huber[胡贝尔]
#[]
dfrpar(tt,"Huber")
z <- airefq(tt, mu=mu, sigmx=1)
#.......................................................................[.................................................. .....................]
{
cat(" airefq : Huber\n reff, beta, nit: ")
cat(round(c(z$reff,z$beta,z$nit),3),sep=c(", ",", ",", ","\n"))
}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Krasker-Welsch[Krasker - 韦尔施]
#[]
dfrpar(tt,"kra-wel",upar=3.755)
z <- airefq(tt, mu=mu, sigmx=1,init=1)
#.......................................................................[.................................................. .....................]
{
cat(" airefq : Krasker-Welsch\n reff, beta, nit: ")
cat(round(c(z$reff,z$beta,z$nit),3),sep=c(", ",", ",", ","\n"))
}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Mallows Standard[锦葵标准]
#[]
dfrpar(tt,"Mal-Std",1.1*(mu+nu),1.569)
z <- airefq(tt, mu=mu, sigmx=1,init=1)
#.......................................................................[.................................................. .....................]
{
cat(" airefq : Mallows-Std\n reff, beta, nit: ")
cat(round(c(z$reff,z$beta,z$nit),3),sep=c(", ",", ",", ","\n"))
}
#[]
# ---- Examples of Chapter 6: Robust testing in linear models ------[----第6章的例子：线性模型的鲁棒测试------]
#[]
#=======================================================================[================================================== =====================]
tautest <- function(x,y,np,nq) {
# Full model. np variables in x[,1:np][完整模式。 NP变量在x [1：NP]]
n <- nrow(x)
z <- riclls(x[,1:np], y)
theta0 <- z$theta; sigma0 <- z$sigma
z <- liepsh(.dFv$ccc) # ccc is globally defined by dfrpar[CCC是全局定义的dfrpar]
epsi2 <- z$epsi2; epsip <- z$epsip
zc <- ktaskv(x[,1:np],f=epsi2/(epsip^2))
covi <- zc$cov
ribeth(y)
zf <- rywalg(x[,1:np],y,theta0,y,covi,sigmai=sigma0)
thetaf <- zf$theta; sigma <- zf$sigmaf; rf <- zf$rs
f <- kffacv(rf,np=np,sigma=sigma)
zc <- ktaskv(x[,1:np],f=f$fh*(sigma^2)/n)
cov <- zc$cov
# Sub-model: nq variables in x[,1:nq], nq < np[子模型中的变量X [1：NQ NQ，NQ NP]
covi <- cov[1 nq*(nq+1)/2)]
zs <- rywalg(x[,1:nq], y, thetaf, y, covi, sigmai=sigma,isigma=0)
thetas <- zs$theta; rs <- zs$rs
# Compute Tau-test statistic and P-value[计算头检验统计量和P-值]
ztau <- tftaut(rf,rs,y,rho,np,nq,sigma)
ftau <- ztau$ftau
z <- chisq(1,np-nq,ftau*epsip/epsi2)
P <- 1-z$p
#.......................................................................[.................................................. .....................]
{
cat(" F_tau, P, sigma: ")
cat(round(c(ftau,P,sigma),3),sep=c(", ",", ","\n"))
cat(" theta (small model): ",round(thetas[1:nq],3),"\n\n")
}
#.......................................................................[.................................................. .....................]
list(thetas=thetas[1:nq], sigma=sigma, rs=rs,ftau=ftau, P=P)}
#------------------------------------------------------------------------[-------------------------------------------------- ----------------------]
dshift <- function(x, theta, sigma, rs, nq) {
# Shift estimate d and confidence interval[换档估计值d和置信区间]
ncov <- nq*(nq+1)/2
f <- kffacv(rs,np=nq,sigma=sigma)
zc <- ktaskv(x[,1:nq],f=f$fh)
cov <- zc$cov
k11 <- 4.*cov[ncov-nq]
k12 <- 2.*cov[ncov-1]
k22 <- cov[ncov]
za <- quant(0.05/2)
d <- 2*theta[nq-1]/theta[nq]
q <- za$x*sigma/theta[nq]
g <- k22*(q^2)
a <- d-g*k12/k22
b <- abs(q)*sqrt(k11-2*d*k12+d*d*k22-g*(k11-k12*k12/k22))
dL <- (a-b)/(1-g)
dU <- (a+b)/(1-g)
#.......................................................................[.................................................. .....................]
cat(" d, dL, dU: ",round(c(d,dL,dU),3),sep=c("",", ",", ","\n"))
#.......................................................................[.................................................. .....................]
list(d=d, dL=dL, dU=dU)
}
#------------------------------------------------------------------------[-------------------------------------------------- ----------------------]
potcy <- function(m,ml,mu,h,k,d,cs,ct) {
fact <- ((h*k) %% 2) + 1
r <- exp(log(d)*(fact*m+h-k)/2 - log(ct/cs))
rl <- exp(log(d)*(fact*ml+h-k)/2 - log(ct/cs))
ru <- exp(log(d)*(fact*mu+h-k)/2 - log(ct/cs))
list(R=r, Rl=rl, Ru=ru)}
#------------------------------------------------------------------------[-------------------------------------------------- ----------------------]
rbmost <- function(x,y,cc,usext=userfd) {
n <- nrow(x); np <- ncol(x); dfcomn(xk=np)
.def <- .dFv; .def$itw <- 1; .dFv <<- .def
z <- wimedv(x)
z <- wyfalg(x, z$a, y, exu=usext); nitw <- z$nit
wgt <- 1/z$dist; wgt[wgt>1.e6] <- 1.e6
z <- comval()
bto <- z$bt0; ipso <- z$ipsi; co <- z$c
z <- ribet0(wgt, itype=2, isqw=0)
xt <- x*wgt; yt <- y * wgt
z <- rilars(xt, yt)
theta0 <- z$theta; sigma0 <- z$sigma
rs <- z$rs/wgt; r1 <- rs/sigma0
dfcomn(ipsi=1,c=cc)
z <- liepsh(cc)
den <- z$epsip
g <- Psp(r1)/den # (see Psp in Chpt. 14)[（在CHPT。看到PSP 14）]
dfcomn(ipsi=ipso, c=co, bet0=bto)
list(theta=theta0, sigma=sigma0, rs=rs, g=g, nitw=nitw)
}
#=======================================================================[================================================== =====================]
dfvals()
z <- c(-1, -2, 0, 35, 1, 0, -3, 20,
-1, -2, 0, 30, 1, 0, -3, 39,
-1, -2, 0, 24, 1, 0, -3, 16,
-1, -2, 0, 37, 1, 0, -3, 27,
-1, -2, 0, 28, 1, 0, -3, -12,
-1, -2, 0, 73, 1, 0, -3, 2,
-1, -2, 0, 31, 1, 0, -3, 31,
-1, -2, 0, 21, 1, 0, -1, 26,
-1, -2, 0, -5, 1, 0, -1, 60,
-1, 0, 0, 62, 1, 0, -1, 48,
-1, 0, 0, 67, 1, 0, -1, -8,
-1, 0, 0, 95, 1, 0, -1, 46,
-1, 0, 0, 62, 1, 0, -1, 77,
-1, 0, 0, 54, 1, 0, 1, 57,
-1, 0, 0, 56, 1, 0, 1, 89,
-1, 0, 0, 48, 1, 0, 1, 103,
-1, 0, 0, 70, 1, 0, 1, 129,
-1, 0, 0, 94, 1, 0, 1, 139,
-1, 0, 0, 42, 1, 0, 1, 128,
-1, 2, 0, 116, 1, 0, 1, 89,
-1, 2, 0, 105, 1, 0, 1, 86,
-1, 2, 0, 91, 1, 0, 3, 140,
-1, 2, 0, 94, 1, 0, 3, 133,
-1, 2, 0, 130, 1, 0, 3, 142,
-1, 2, 0, 79, 1, 0, 3, 118,
-1, 2, 0, 120, 1, 0, 3, 137,
-1, 2, 0, 124, 1, 0, 3, 84,
-1, 2, 0, -8, 1, 0, 3, 101)
xx <- matrix(z,ncol=4, byrow=TRUE)
dimnames(xx) <- list(NULL,c("z2","xS","xT","y"))
z2 <- xx[,"z2"]; xS <- xx[,"xS"]; xT <- xx[,"xT"]
x <- cbind(1, z2, xS+xT, xS-xT, xS^2+xT^2, xS^2-xT^2, xT^3)
y <- xx[,"y"]
z <- dfrpar(x, "huber",psipar=1.345)
#[]
# Tau-test and shift estimate[头测试和移位估计]
#[]
{
cat("Results (linearity test)\n")
np <- 7; nq <- 4 # Test linearity[测试线性]
z <- tautest(x,y,np,nq)
cat("Results (parallelism test)\n")
np <- 4; nq <- 3 # Test parallelism[测试的并行]
z <- tautest(x,y,np,nq)
z <- dshift(x, z$thetas, z$sigma, z$rs, nq)
}
#------------------------------------------------------------------------[-------------------------------------------------- ----------------------]
# Input data; set defaults[输入数据，默认设置]
# []
z <- c(35.3, 20, 10.98,
29.7, 20, 11.13,
30.8, 23, 12.51,
58.8, 20, 8.40,
61.4, 21, 9.27,
71.3, 22, 8.73,
74.4, 11, 6.36,
76.7, 23, 8.50,
70.7, 21, 7.82,
57.5, 20, 9.14,
46.4, 20, 8.24,
28.9, 21, 12.19,
28.1, 21, 11.88,
39.1, 19, 9.57,
46.8, 23, 10.94,
48.5, 20, 9.58,
59.3, 22, 10.09,
70.0, 22, 8.11,
70.0, 11, 6.83,
74.5, 23, 8.88,
72.1, 20, 7.68,
58.1, 21, 8.47,
44.6, 20, 8.86,
33.4, 20, 10.36,
28.6, 22, 11.08)
x <- matrix(z, ncol=3, byrow=TRUE)
y <- x[,3]; x[,2:3] <- x[,1:2]; x[,1] <- 1
n <- length(y); np <- ncol(x); nq <- np - 1
#[]
# Optimal tau-test based on Schweppe-type estimates[最优头试验的基础上施韦普型估计]
#[]
z <- tauare(itype=3,mu=1,cpsi=2.665,bb=0,sigmax=1)
dfrpar(x, "Sch-Tau",upar=2.67); .dFv$isg <- 1; dfcomn(d=.dFv$ccc)
# Full model: initial estimates of theta, sigma and weights[完整模式：THETA，标准差和重量的初步估计]
dfcomn(xk=np) # Needed for userfd [所需的userfd]
zr <- rbmost(x,y,cc=1.5)
theta0 <- zr$theta; sigma0 <- zr$sigma; nitw0 <- zr$nitw
#[]
# Initial and final values of weights [的权重的初始值和最终值]
#[]
z <- wimedv(x)
z <- wyfalg(x, z$a, zr$g, nvarq=nq, igwt=1)
wgt <- 1/z$dist ; wgt[wgt>1.e6] <- 1.e6
# Full model: covariance matrix of coefficients and inverse[完整模式：系数和协方差矩阵的逆]
z <- kiedch(wgt)
zc <- ktaskw(x, z$d, z$e, f=1/n, iainv=1)
cov <- zc$cov; ainv <- zc$ainv
# Full model: Final estimate of theta and sigma[完整模式：最后估计的θ波和sigma]
ribeth(wgt)
zf <- rywalg(x, y,theta0,wgt,cov, sigmai=sigma0)
thetaf <- zf$theta[1:np]; sigma <- zf$sigmaf; nitf <- zf$nit
# Small model: Final estimate of theta and sigma[小模型：最后估计的θ波和sigma]
covi <- cov[1nq*(nq+1)/2)]
xt <- x[,1:nq,drop=FALSE]
zs <- rywalg(xt, y, theta0, wgt, covi, sigmai=sigma, isigma=0)
thetas <- zs$theta[1:nq]; nits <- zs$nit
# Compute Tau-test statistic[计算头检验统计量]
ft <- tftaut(zf$rs,zs$rs,wgt,rho,np,nq,sigma)
ftau <- ft$ftau
# P-value[P-值]
z <- ttaskt(cov, ainv, np, nq, fact=n)
z <- tteign(z$covtau,nq)
xlmbda <- z$xlmbda[1np-nq)]
mult <- rep(1, length=np)
delta <- rep(0, length=np)
z <- ruben(xlmbda, delta, mult,ftau, xmode=-1, maxit=100, eps=0.0001)
P <- 1-z$cumdf
#.......................................................................[.................................................. .....................]
{
cat(" Optimal Tau-test : Schweppe-type estimates\n")
cat(" theta0: ",round(theta0[1:np],3),"\n sigma0, nitw: ")
cat(round(c(sigma0,nitw0),3),sep=c(", ","\n"))
cat(" thetaf: ",round(thetaf,3),"\n sigma, nit: ")
cat(round(c(sigma,nitf),3),sep=c(", ","\n"))
cat(" thetas: ",round(thetas,3),"\n sigma, nit: ")
cat(round(c(sigma,nits),3),sep=c(", ","\n"))
cat(" F_tau =",round(ftau,3),", P =",P,"\n")
}
#=======================================================================[================================================== =====================]
rn2mal <- function(x,y,b2,c,nq) {
n <- length(y); np <- ncol(x)
#[]
# Rn2-test on Mallows estimators [RN2测试锦葵估计]
# ==============================[==============================]
dfrpar(x, "mal-std", b2, c)
.def <- .dFv
.def$isg <- 1; .dFv <<- .def; dfcomn(d=.dFv$ccc)
#[]
# Initial and final values of weights[的权重的初始值和最终值]
#[]
z <- wimedv(x)
z <- wyfalg(x, z$a, wgt)
wgt <- Www(z$dist); nitw <- z$nit
#[]
# Initial theta and sigma (using weighted LAR)[最初的theta和西格玛（采用加权LAR）]
#[]
ribet0(wgt, isqw=0)
xt <- x*wgt
yt <- y * wgt
z <- rilars(xt, yt)
theta0 <- z$theta; sigma0 <- z$sigma
#[]
# Initial value of COV[初始值的COV]
#[]
z <- kiedch(wgt)
zc <- ktaskw(x, z$d, z$e, f=1/n)
covi <- zc$cov
#[]
# Solution by means of RYWALG. [解决手段RYWALG。]
#[]
z <- ribeth(wgt)
beta <- z$bta
zw <- rywalg(x, y, theta0, wgt, covi, sigmai=sigma0)
theta1 <- zw$theta[1:np]; sigma1 <- zw$sigmaf; nit1 <- zw$nit
#[]
# Unscaled covariance matrix of coefficients[系数的未缩放的协方差矩阵]
#[]
zc <- kfedcb(wgt, zw$rs, sigma=sigma1)
z <- ktaskw(x, zc$d, zc$e, f=1/n)
cov1 <- z$cov
#[]
# Rn2-test statistic and significance[RN2检验统计量和意义]
#[]
z <- tfrn2t(cov1,theta1,n,nq)
rn2m <- z$rn2t/(n*sigma1^2)
z <- chisq(1,np-nq,rn2m)
p1 <- 1.-z$p
list(theta1=theta1, sigma1=sigma1, wgt=wgt, nitw=nitw, nit1=nit1,
rn2m=rn2m, p1=p1)}
#------------------------------------------------------------------------[-------------------------------------------------- ----------------------]
#[]
# Read data[读取数据]
# []
z <- c(35.3, 20, 10.98,
29.7, 20, 11.13,
30.8, 23, 12.51,
58.8, 20, 8.40,
61.4, 21, 9.27,
71.3, 22, 8.73,
74.4, 11, 6.36,
76.7, 23, 8.50,
70.7, 21, 7.82,
57.5, 20, 9.14,
46.4, 20, 8.24,
28.9, 21, 12.19,
28.1, 21, 11.88,
39.1, 19, 9.57,
46.8, 23, 10.94,
48.5, 20, 9.58,
59.3, 22, 10.09,
70.0, 22, 8.11,
70.0, 11, 6.83,
74.5, 23, 8.88,
72.1, 20, 7.68,
58.1, 21, 8.47,
44.6, 20, 8.86,
33.4, 20, 10.36,
28.6, 22, 11.08)
x <- matrix(z, ncol=3, byrow=TRUE)
y <- x[,3]; x[,2:3] <- x[,1:2]; x[,1] <- 1
n <- length(y); np <- ncol(x); nq <- np - 1
wgt <- vector("numeric",length(y))
z <- rn2mal(x, y, 4, 1.5, nq)
#.......................................................................[.................................................. .....................]
{
cat("Rn2-test on Mallows estimators\n")
cat(" theta1: ",round(z$theta1,3),"\n sigma1, nitw, nit1: ")
cat(round(c(z$sigma1,z$nitw,z$nit1),3),sep=c(", ",", ","\n"))
cat(" Rn2 =",round(z$rn2m,3),", P =",z$p1,"\n")
}
#[]
# ---- Examples of Chapter 7: Breakdown point regression ------[----第7章的例子：故障点回归------]
#[]
#[]
# Read data; load defaults[读取数据载入预设值]
# []
z <- c(80, 27, 89, 42,
80, 27, 88, 37,
75, 25, 90, 37,
62, 24, 87, 28,
62, 22, 87, 18,
62, 23, 87, 18,
62, 24, 93, 19,
62, 24, 93, 20,
58, 23, 87, 15,
58, 18, 80, 14,
58, 18, 89, 14,
58, 17, 88, 13,
58, 18, 82, 11,
58, 19, 93, 12,
50, 18, 89, 8,
50, 18, 86, 7,
50, 19, 72, 8,
50, 19, 79, 8,
50, 20, 80, 9,
56, 20, 82, 15,
70, 20, 91, 15)
x <- matrix(z, ncol=4, byrow=TRUE)
y <- x[,4]; x[,4] <- 1
n <- length(y); np <- ncol(x)
nq <- np+1
dfvals()
#[]
# Least median of squares[至少中位数的平方]
#[]
zr <- hylmse(x,y, nq, ik=1, iopt=3, intch=1)
theta1 <- zr$theta; xmin1 <- zr$xmin
zr <- hylmse(x,y, nq, ik=2, iopt=3, intch=1)
theta2 <- zr$theta; xmin2 <- zr$xmin
zr <- hylmse(x,y, nq, ik=1, iopt=1, intch=1)
theta3 <- zr$theta; xmin3 <- zr$xmin
#[]
# Least trimmed squares[最不修剪广场]
#[]
zr <- hyltse(x,y, nq, ik=1, iopt=3, intch=1)
theta4 <- zr$theta; xmin4 <- zr$smin
zr <- hyltse(x,y, nq, ik=2, iopt=3, intch=1)
theta5 <- zr$theta; xmin5 <- zr$smin
zr <- hyltse(x,y, nq, ik=1, iopt=1, intch=1)
theta6 <- zr$theta; xmin6 <- zr$smin
#[]
# S-estimate[S-估计]
#[]
z <- dfrpar(x,'S')
z <- ribetu(y)
zr <- hysest(x,y, nq, iopt=3, intch=1)
theta7 <- zr$theta[1:np]; xmin7 <- zr$smin
zr <- hysest(x,y, nq, iopt=1, intch=1)
theta8 <- zr$theta[1:np]; xmin8 <- zr$smin
#.......................................................................[.................................................. .....................]
{
cat("Results\n theta1 = (")
cat(round(theta1,3),sep=", ")
cat("), xmin1 ="); cat(round(xmin1,3))
cat("\n theta2 = ("); cat(round(theta2,3),sep=", ")
cat("), xmin2 ="); cat(round(xmin2,3))
cat("\n theta3 = ("); cat(round(theta3,3),sep=", ")
cat("), xmin3 ="); cat(round(xmin3,3))
cat("\n theta4 = ("); cat(round(theta4,3),sep=", ")
cat("), xmin4 ="); cat(round(xmin4,3))
cat("\n theta5 = ("); cat(round(theta5,3),sep=", ")
cat("), xmin5 ="); cat(round(xmin5,3))
cat("\n theta6 = ("); cat(round(theta6,3),sep=", ")
cat("), xmin6 ="); cat(round(xmin6,3))
cat("\n theta7 = ("); cat(round(theta7,3),sep=", ")
cat("), xmin7 ="); cat(round(xmin7,3))
cat("\n theta8 = ("); cat(round(theta8,3),sep=", ")
cat("), xmin8 ="); cat(round(xmin8,3),"\n")
}
#[]
# ---- Examples of Chapter 8: M-estimates of covariance matrices ------[----第8章的例子：M-估计的协方差矩阵------]
#[]
#[]
# Read data; set defaults[读取数据，默认设置]
# []
z <- c(4.37, 5.23, 4.38, 5.02,
4.56, 5.74, 4.42, 4.66,
4.26, 4.93, 4.29, 4.66,
4.56, 5.74, 4.38, 4.90,
4.30, 5.19, 4.22, 4.39,
4.46, 5.46, 3.48, 6.05,
3.84, 4.65, 4.38, 4.42,
4.57, 5.27, 4.56, 5.10,
4.26, 5.57, 4.45, 5.22,
4.37, 5.12, 3.49, 6.29,
3.49, 5.73, 4.23, 4.34,
4.43, 5.45, 4.62, 5.62,
4.48, 5.42, 4.53, 5.10,
4.01, 4.05, 4.45, 5.22,
4.29, 4.26, 4.53, 5.18,
4.42, 4.58, 4.43, 5.57,
4.23, 3.94, 4.38, 4.62,
4.42, 4.18, 4.45, 5.06,
4.23, 4.18, 4.50, 5.34,
3.49, 5.89, 4.45, 5.34,
4.29, 4.38, 4.55, 5.54,
4.29, 4.22, 4.45, 4.98,
4.42, 4.42, 4.42, 4.50,
4.49, 4.85)
cx <- matrix(z, ncol=2, byrow=TRUE)
n <- nrow(cx); np <- ncol(cx)
dst0 <- vector("numeric",n)
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Classical covariance[古典协方差]
#[]
t0 <- apply(cx, 2, mean)
xmb <- sweep(cx, 2, t0)
cv0 <- crossprod(xmb)/n
# Mahalanobis distances[马氏距离]
cvm1 <- solve(cv0)
for (i in 1:n) {
z <- xmb[i,,drop=FALSE]; dst0[i] <- sqrt(z %*% cvm1 %*% t(z))}
#=======================================================================[================================================== =====================]
# M-estimate of covariance[M-估计的协方差]
#[]
zc <- cicloc()
za <- cia2b2(nvar=np)
a2 <- za$a2; b2 <- za$b2
zd <- cibeat(a2, b2, np)
cw <- zc$c; dv <- zd$d
dfcomn(iucv=1, a2=a2, b2=b2, bt=dv, cw=cw)
# zf <- cifact(a2,b2,np); fc <- zf$fc[采埃孚 - cifact（A2，B2，NP）; FC < - ZF $ FC]
z <- cimedv(cx)
ai <- z$a; ti <- z$t; fc <- 1
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# With prescription F0[随着处方F0]
zd <- cyfalg(cx,ai,ti)
zc <- cfrcov(zd$a,np,fc)
cv1 <- zc$cov; t1 <- zd$t; dst1 <- zd$dist; nt1 <- zd$nit
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# With prescription NH[随着处方NH]
zd <- cynalg(cx,ai,ti)
zc <- cfrcov(zd$a,np,fc)
cv2 <- zc$cov; t2 <- zd$t; dst2 <- zd$dist; nt2 <- zd$nit
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# With prescription CG[随着处方CG]
zd <- cygalg(cx,ai,ti)
zc <- cfrcov(zd$a,np,fc)
cv3 <- zc$cov; t3 <- zd$t; dst3 <- zd$dist; nt3 <- zd$nit
#.......................................................................[.................................................. .....................]
{
cat("Results\n\n cv0[1,1],cv0[2,1],cv0[2,2] = (")
cat(round(as.vector(cv0)[-2],3),sep=", ")
cat(")\n t0 = ("); cat(round(t0,3),sep=", ")
cat(")\n dist0 :\n ")
cat(round(dst0,3),sep=c(rep(", ",9),",\n "))
cat("\n cv1[1,1],cv1[2,1],cv1[2,2] = (")
cat(round(cv1,3),sep=", ")
cat(")\n t1 = ("); cat(round(t1,3),sep=", ")
cat("), nit1 =",nt1); cat("\n dist1 :\n ")
cat(round(dst1,3),sep=c(rep(", ",9),",\n "))
cat("\n cv2[1,1],cv2[2,1],cv2[2,2] = (")
cat(round(cv2,3),sep=", ")
cat(")\n t2 = ("); cat(round(t2,3),sep=", ")
cat("), nit2 =",nt2); cat("\n dist2 :\n ")
cat(round(dst2,3),sep=c(rep(", ",9),",\n "))
cat("\n cv3[1,1],cv3[2,1],cv3[2,2] = (")
cat(round(cv3,3),sep=", ")
cat(")\n t3 = ("); cat(round(t3,3),sep=", ")
cat("), nit3 =",nt3); cat("\n dist3 :\n ")
cat(round(dst3,3),sep=c(rep(", ",9),",\n "))
}
#[]
# ----------- Examples of Chapter 9: Mixed procedures --------------[-----------第9章：混合程序的例子--------------]
#[]
bindec <- function(np,ind,cpc,cpr) {
n <- length(ind)
ccar <- matrix("-",ncol=np, nrow=n)
for (i in 1:n) {
j <- 0
num <- abs(ind[i])
while (num != 0 & j < np) {
j <- j+1
if (num %% 2 == 1) ccar[i,j] <- "X"
num <- num %/% 2}}
data.frame(Cp=round(cpc,3),Cp.r=round(cpr,3),ipr=ind,i=ccar)
}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Read data[读取数据]
# []
z <- c(-1, -2, 0, 35, 1, 0, -3, 20,
-1, -2, 0, 30, 1, 0, -3, 39,
-1, -2, 0, 24, 1, 0, -3, 16,
-1, -2, 0, 37, 1, 0, -3, 27,
-1, -2, 0, 28, 1, 0, -3, -12,
-1, -2, 0, 73, 1, 0, -3, 2,
-1, -2, 0, 31, 1, 0, -3, 31,
-1, -2, 0, 21, 1, 0, -1, 26,
-1, -2, 0, -5, 1, 0, -1, 60,
-1, 0, 0, 62, 1, 0, -1, 48,
-1, 0, 0, 67, 1, 0, -1, -8,
-1, 0, 0, 95, 1, 0, -1, 46,
-1, 0, 0, 62, 1, 0, -1, 77,
-1, 0, 0, 54, 1, 0, 1, 57,
-1, 0, 0, 56, 1, 0, 1, 89,
-1, 0, 0, 48, 1, 0, 1, 103,
-1, 0, 0, 70, 1, 0, 1, 129,
-1, 0, 0, 94, 1, 0, 1, 139,
-1, 0, 0, 42, 1, 0, 1, 128,
-1, 2, 0, 116, 1, 0, 1, 89,
-1, 2, 0, 105, 1, 0, 1, 86,
-1, 2, 0, 91, 1, 0, 3, 140,
-1, 2, 0, 94, 1, 0, 3, 133,
-1, 2, 0, 130, 1, 0, 3, 142,
-1, 2, 0, 79, 1, 0, 3, 118,
-1, 2, 0, 120, 1, 0, 3, 137,
-1, 2, 0, 124, 1, 0, 3, 84,
-1, 2, 0, -8, 1, 0, 3, 101)
xx <- matrix(z,ncol=4, byrow=TRUE)
dimnames(xx) <- list(NULL,c("z2","xS","xT","y"))
z2 <- xx[,"z2"]; xS <- xx[,"xS"]; xT <- xx[,"xT"]
x <- cbind(1, z2, xS+xT, xS-xT, xS^2+xT^2, xS^2-xT^2, xT^3)
y <- xx[,"y"]
wgt <- vector("numeric",length(y))
n <- 56; np <- 7
dfvals()
# Compute classical sigma and the t-statistics[计算经典的标准差和t-统计]
dfrpar(x,"ols",-1,-1)
z <- mirtsr(x,y,.dFv$ite)
sigmc <- z$sigma; tstac <- z$t
# Compute robust sigma and the t-statistics[计算强大的标准差和t-统计]
dfrpar(x,"huber",-1,-1)
z <- mirtsr(x,y,.dFv$ite)
sigmr <- z$sigma; tstar <- z$t
#[]
# All possible regressions including the constant and linear terms[包括所有可能回归的常数及线性]
#[]
vp <- rep(-0.5, length=np)
vp[1] <- 3; vp[3] <- 2; vp[4] <- 1
za <- mfragr(x, y, vp, nc=18, .dFv$ite, sigmac=sigmc, sigmar=sigmr)
#[]
# Priorites by means of t-directed search[叔定向搜索由装置的优先权的]
#[]
zt <- mfragr(x, y, tstar, nc=7, .dFv$ite, sigmac=sigmc, sigmar=sigmr)
#.......................................................................[.................................................. .....................]
{
cat(" Estimates of sigma\n ")
cat(" sigmc =",round(sigmc,3),", sigmr =",round(sigmr,3),"\n")
cat(" Regressions on subset of variables:\n")
cat(" C{p} C{p,@} ipr 1 2 3 4 5 6 7\n")
cat(t(bindec(np,za$ipr,za$cpc,za$cpr)),sep=c(rep(" ",9),"\n"))
cat("\n t-directed search\n")
cat(" tstar[1:7]=(", round(tstar,3),sep=c("",rep(", ",6)))
cat(")\n C_p C{p,@} ipr 1 2 3 4 5 6 7\n")
cat(t(bindec(np,zt$ipr,zt$cpc,zt$cpr)),sep=c(rep(" ",9),"\n"))
}
#=======================================================================[================================================== =====================]
#[]
# Read data; set defaults[读取数据，默认设置]
# []
z <- c(4.37, 5.23, 4.48, 5.42, 4.38, 5.02, 4.53, 5.10,
4.56, 5.74, 4.01, 4.05, 4.42, 4.66, 4.45, 5.22,
4.26, 4.93, 4.29, 4.26, 4.29, 4.66, 4.53, 5.18,
4.56, 5.74, 4.42, 4.58, 4.38, 4.90, 4.43, 5.57,
4.30, 5.19, 4.23, 3.94, 4.22, 4.39, 4.38, 4.62,
4.46, 5.46, 4.42, 4.18, 3.48, 6.05, 4.45, 5.06,
3.84, 4.65, 4.23, 4.18, 4.38, 4.42, 4.50, 5.34,
4.57, 5.27, 3.49, 5.89, 4.56, 5.10, 4.45, 5.34,
4.26, 5.57, 4.29, 4.38, 4.45, 5.22, 4.55, 5.54,
4.37, 5.12, 4.29, 4.22, 3.49, 6.29, 4.45, 4.98,
3.49, 5.73, 4.42, 4.42, 4.23, 4.34, 4.42, 4.50,
4.43, 5.45, 4.49, 4.85, 4.62, 5.62)
cx <- matrix(z, ncol=2, byrow=TRUE)
n <- nrow(cx); np <- ncol(cx)
y <- vector("numeric",length=n)
#[]
# Minimum Volume Ellipsoid covariances [最小体积椭球的协方差]
#[]
dfvals()
z <- mymvlm(cx,y,ilms=0,iopt=3,iseed=5321)
dst <- z$d; cv <- z$cov; xvol <- z$xvol
#.......................................................................[.................................................. .....................]
{
cat("Minimum Volume Ellipsoid covariances\n cv = (")
cat(round(cv,3),sep=c(", ",", "))
cat("), Objective function value =",round(xvol,3),"\ndistances:\n")
cat(round(dst,3),sep=c(rep(", ",9),",\n"))
}
#======================================================================= [================================================== =====================]
#[]
# Read data; load defaults[读取数据载入预设值]
# []
z <- c(80, 27, 89, 42,
80, 27, 88, 37,
75, 25, 90, 37,
62, 24, 87, 28,
62, 22, 87, 18,
62, 23, 87, 18,
62, 24, 93, 19,
62, 24, 93, 20,
58, 23, 87, 15,
58, 18, 80, 14,
58, 18, 89, 14,
58, 17, 88, 13,
58, 18, 82, 11,
58, 19, 93, 12,
50, 18, 89, 8,
50, 18, 86, 7,
50, 19, 72, 8,
50, 19, 79, 8,
50, 20, 80, 9,
56, 20, 82, 15,
70, 20, 91, 15)
x <- matrix(z, ncol=4, byrow=TRUE)
y <- x[,4]; x[,4] <- 1
n <- length(y); np <- ncol(x)
nq <- np+1
dfvals()
#[]
# High breakdown point & high efficiency regression[高击穿点和高效率的回归]
#[]
dfrpar(x,"S",-1,-1)
z <- myhbhe(x, y, iseed=5431)
#.......................................................................[.................................................. .....................]
{
cat("High breakdown point & high efficiency regression\n")
cat(" theta0 = ("); cat(round(z$theta0,3),sep=", ")
cat("), sigma0 =",round(z$sigm0,3))
cat("\n theta1 = ("); cat(round(z$theta1,3),sep=", ")
cat("), sigma1 = ",round(z$sigm1,3),", tbias =",sep="")
cat(round(z$tbias,3),"\n")
}
#[]
# -------- Examples of Chapter 10: M-estimates for discrete GLM ---------[--------第10章的例子：M-估计的离散GLM ---------]
#[]
glmb <- function(x,y,n,np,upar) {
#[]
# BI estimates of theta, A, ci and wa: Bernouilli responses, b=upar[BI估计θ，A，CI和WA：伯努利反应，B = UPAR]
#[]
# Initial theta, A (A0) and c (c0)[初步θ，A（A0）和c（C0）]
#[]
ni <- vector("integer",n)
z <- gintac(x,y,ni,icase=1,b=upar,c=1.5)
theta0 <- z$theta[1:np]; A0 <- z$a; c0 <- z$ci
# Initial distances |Ax_i| and cut off points a_i (wa)[初始距离| Ax_i并切断，点A_I（WA）]
# for (i in 1:n)[（1：N）]
# {ax <- mlyd(A0,x[i,],np) # See Chpt. 12 for mlyd [{斧头< - mlyd（A0，X [I]，NP），＃看看CHPT。 12 mlyd]
# znr <- ax$y %*% t(ax$y)[氧化锌非线性电阻< - 斧头$ Y％％T（AX $ Y）]
# wa[i] <- upar/sqrt(znr)}[WA [我] < - 的UPAR / SQRT（氧化锌非线性电阻）}]
wa <- upar/sqrt(z$dist) # ==> Modif. in gintac[==> Modif。在gintac]
vtheta <- x %*% theta0
z <- gfedca(vtheta, c0, wa, ni, icase=1)
zc <- ktaskw(x, z$d, z$e, f=1/n) # See Chpt. 4[请参阅CHPT。 4]
covi <- zc$cov
# Final theta, A, c (ci) and a (wa)[最终θ，A，C（CI）和（WA）]
z <- gymain(x, y, ni, covi, A0, theta0, b = upar)
theta <- z$theta; A <- z$a; ci <- z$ci; wa <- z$wa; nit <- z$nit
# Final cov. matrix and std. dev's of coeff. estimates[最后的冠状病毒。矩阵和性传播疾病。开发的COEFF。估计]
z <- gfedca(z$vtheta, ci, wa, ni, icase=1)
sdev <- NULL
zc <- ktaskw(x, z$d, z$e, f=1/n)
for (i in 1:np) {ii <- i*(i+1)/2; sdev <- c(sdev,zc$cov[ii])}
sdev <- sqrt(sdev)
list(theta=theta, A=A, ci=ci, wa=wa, nit=nit, sdev=sdev)}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
# Read data; load defaults[读取数据载入预设值]
# []
z <- c(3.70, 0.825, 1, 3.50, 1.090, 1,
1.25, 2.500, 1, 0.75, 1.500, 1,
0.80, 3.200, 1, 0.70, 3.500, 1,
0.60, 0.750, 0, 1.10, 1.700, 0,
0.90, 0.750, 0, 0.90, 0.450, 0,
0.80, 0.570, 0, 0.55, 2.750, 0,
0.60, 3.000, 0, 1.40, 2.330, 1,
0.75, 3.750, 1, 2.30, 1.640, 1,
3.20, 1.600, 1, 0.85, 1.415, 1,
1.70, 1.060, 0, 1.80, 1.800, 1,
0.40, 2.000, 0, 0.95, 1.360, 0,
1.35, 1.350, 0, 1.50, 1.360, 0,
1.60, 1.780, 1, 0.60, 1.500, 0,
1.80, 1.500, 1, 0.95, 1.900, 0,
1.90, 0.950, 1, 1.60, 0.400, 0,
2.70, 0.750, 1, 2.35, 0.300, 0,
1.10, 1.830, 0, 1.10, 2.200, 1,
1.20, 2.000, 1, 0.80, 3.330, 1,
0.95, 1.900, 0, 0.75, 1.900, 0,
1.30, 1.625, 1)
x <- matrix(z, ncol=3, byrow=TRUE)
y <- x[,3]; x[,3] <- log(x[,2]); x[,2] <- log(x[,1]) ; x[,1] <- 1
n <- length(y); np <- ncol(x)
dfvals()
upar <- 3.2*sqrt(np)
z1 <- glmb(x,y,n,np,upar)
upar <- 3.7*sqrt(np)
z2 <- glmb(x,y,n,np,upar)
z <- glmb(x,y,n,np,300) # Classical estimates[古典估计]
#.......................................................................[.................................................. .....................]
{
cat("\n Robust estimates : upar=5.5426, nit =",z1$nit,"\n")
cat(" {theta[i] (sdev[i]), i=1:3}\n ")
for (i in 1:3) cat(round(z1$theta[i],3)," (",round(z1$sdev[i],3),") ",sep="")
cat("\n\n Robust estimates : upar=6.4086, nit =",z2$nit,"\n")
cat(" {theta[i] (sdev[i]), i=1:3}\n ")
for (i in 1:3) cat(round(z2$theta[i],3)," (",round(z2$sdev[i],3),") ",sep="")
cat("\n\n Classical estimates : upar=300, nit =",z$nit,"\n")
cat(" {theta[i] (sdev[i]), i=1:3}\n ")
for (i in 1:3) cat(round(z$theta[i],3)," (",round(z$sdev[i],3),") ",sep="")
cat("\n")
}
#-----------------------------------------------------------------------[-------------------------------------------------- ---------------------]
转载请注明:出自 生物统计家园网(http://www.biostatistic.net)。
注：
注1：为了方便大家学习，本文档为生物统计家园网机器人LoveR翻译而成，仅供个人R语言学习参考使用，生物统计家园保留版权。
注2：由于是机器人自动翻译，难免有不准确之处，使用时仔细对照中、英文内容进行反复理解，可以帮助R语言的学习。
注3：如遇到不准确之处，请在本贴的后面进行回帖，我们会逐渐进行修订。
|
发表于 2012-9-27 20:53:34
