version 0.2.0

cran · Feb 13, 2019 · 43b4969 · 43b4969
1 parent 8de09b2
commit 43b4969
Show file tree

Hide file tree

Showing 10 changed files with 720 additions and 260 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,17 +1,16 @@
 Package: rlmDataDriven
 Type: Package
 Title: Robust Regression with Data Driven Tuning Parameter
-Version: 0.1.0
+Version: 0.2.0
 Author: You-Gan Wang
 Maintainer: The package maintainer <you-gan.wang@qut.edu.au>
-Imports: stats, MASS
-Description: Data driven approach for robust regression estimation.
-              See Wang et al. (2007), <doi:10.1198/106186007X180156>.
+Imports: stats, MASS, tseries
+Description: Data driven approach for robust regression estimation in homoscedastic and heteroscedastic context. See Wang et al. (2007), <doi:10.1198/106186007X180156> regarding homoscedastic framework.  
 License: GPL (>= 2.0)
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 6.1.0
 NeedsCompilation: no
-Packaged: 2018-08-10 04:37:08 UTC; wangy
+Packaged: 2019-02-13 04:45:42 UTC; liquetwe
 Repository: CRAN
-Date/Publication: 2018-09-17 15:10:03 UTC
+Date/Publication: 2019-02-13 06:50:04 UTC
diff --git a/MD5 b/MD5
@@ -1,8 +1,12 @@
-3f3e1c776ddd3574b88183e68b2d686a *DESCRIPTION
-e684857efafa7b194d6d0459581b4b22 *NAMESPACE
-2bae9720c36a6d5ea17f668691892b0f *R/DD-internal.R
+8db8b586ba39d089c62d7346ca292afc *DESCRIPTION
+83af3ef5dc9d34f8c8636e908f2db9af *NAMESPACE
+a856d9caba372bc66122240d916e9754 *R/DD-internal.R
 2801b315c3d0a5a1cb7b2af36ff8b95a *R/rlmDD.R
+98b0f6d67265eb93b6cdda6045faffc5 *R/rlmDD_het.R
+611920c4014bf0220079dab79a9f8d38 *R/whm.R
 86db0b2c4c3a8a7f44ae40bbaf3ab29e *data/plasma.RData
-4f1e7315c723f966f354edb580e1ef61 *man/DD-internal.Rd
+1cbee3fec615e7b13578f68d8ffc1f75 *man/DD-internal.Rd
 887918ce9fb1ef10e70d2bcaff33e50e *man/plasma.Rd
-b09ddfb1b3755e554a3d4070a3dca96a *man/rlmDD.Rd
+ba7fcd8d357187f8d5fc09425658376f *man/rlmDD.Rd
+0ed97475764a9ede18e52152941f939a *man/rlmDD_het.Rd
+505e7f3a60fe4534f1195efd8c708c34 *man/whm.Rd
diff --git a/NAMESPACE b/NAMESPACE
@@ -1,14 +1,15 @@
 exportPattern("^[[:alpha:]]+")
-importFrom("utils", "data")
+importFrom("utils", "data", "setTxtProgressBar", "txtProgressBar")
 importFrom("MASS","rlm")
 importFrom("MASS","psi.huber")
 importFrom("MASS","psi.bisquare")
 importFrom("stats", "lm", "median")
 importFrom("graphics", "axis", "legend")
-importFrom("stats", "mad")
+importFrom("stats", "mad", "pacf", "uniroot", "var")
+import("tseries")
 export(eff,ESL_O)
 export(rho.h,psi.Huber,dpsi.Huber)
 export(rho.b,psi.Tukey,dpsi.Tukey)
 export(rho.e,psi.Exp,dpsi.Exp)
-
+export(chi,create_lag)
 
diff --git a/R/DD-internal.R b/R/DD-internal.R
@@ -1,127 +1,142 @@
-##huber's function
-psi.Huber<-function(r,c)
-{
-  newr<-pmin(c,pmax(r,-c))
-  newr
-}
-
-##derivative of huber's function
-dpsi.Huber<-function(r,c)
-{
-  newr<-as.numeric(abs(r)<=c)
-  newr
-}
-
-
-##bisquare's function
-psi.Tukey<-function(r,c)
-{
-  r<-pmin(c,pmax(r,-c))
-  newr<-r*(1-(r/c)^2)^2
-  newr
-}
-
-##derivative of bisquare's function
-dpsi.Tukey<-function(r,c)
-{
-  r<-pmin(c,pmax(r,-c))
-  newr<-(1-(r/c)^2)*(1-5*(r/c)^2)
-  newr
-}
-
-##modified huber's function
-psi.Exp<-function(r,c)
-{
-  newr<-r*exp(-(r/c)^2)
-  newr
-}
-
-##derivative of modified huber's function
-dpsi.Exp<-function(r,c)
-{
-  newr<-(1-(2*r^2)/c^2)*exp(-(r/c)^2)
-  newr
-}
-
-
-##find the most efficient tuning constant
-eff<-function(r,method,plot)
-{
-  nn<-  length(r)
-  if (method =="Huber")
-  {
-    tau<-(1:30)/10
-    new<-unlist(lapply(tau,FUN=function(x){(sum(dpsi.Huber(r,x)))^2/(nn*sum(psi.Huber(r,x)^2))}))
-  }
-  else if (method=="Bisquare")
-  {
-    tau<-(15.48:59.48)/10
-    new<-unlist(lapply(tau,FUN=function(x){(sum(dpsi.Tukey(r,x)))^2/(nn*sum(psi.Tukey(r,x)^2))}))
-  }
-  else if (method=="Exponential")
-  {
-    tau<-(5:100)/10
-    new<-unlist(lapply(tau,FUN=function(x){(sum(dpsi.Exp(r,x)))^2/(nn*sum(psi.Exp(r,x)^2))}))
-  }
-
-  if (plot =="Y")
-  {
-    plot(new, type="o", xaxt = "n", xlab="Tunning parameter",
-        ylab="Efficiency",ylim=c(0,max(new)*1.3))
-    axis(1, at=seq_along(tau), labels=tau)
-    if (method =="Exponential") {
-      legend(length(tau)*0.8,max(new)*1.3,c("Exp"),bty = "n")
-    } else {
-      legend(length(tau)*0.8,max(new)*1.3,c(method),bty = "n")
-    }
-  }
-
-  xx<- order(new)
-  return(tau[xx[length(tau)]])
-
-}
-
-
-rho.h<-function(u,c){
-  phi<-0;
-  x1<-(abs(u)<=c)
-  phi<-x1* (u^2/2) + (1-x1)*(abs(u)*c-c^2/2)
-  phi
-}
-
-rho.b<-function(u,c){
-  phi<-0;
-  x1<-(abs(u)<=c)
-  phi<-x1*(1-(1-(u/c)^2)^3) + (1-x1)
-  phi
-}
-
-rho.e<-function(u,c){
-  1-exp(-(u/c)^2)
-}
-
-
-ESL_O<-function(x,xx,y,beta,newc,maxit=500, toler=1e-6)
-  {
-  it=0;delta=1;
-  while (delta>toler && it<maxit){
-    it<-it+1;
-    beta0<-beta;
-    w<-0;
-    e<-y-x%*%beta0;
-    S_n<-max(mad(e),1e-6)
-    u<-e/S_n
-    for(i in seq_along(u)){
-      w[i]<-exp(-(u[i]/newc)^2);
-    }
-
-    W<-diag(w);
-    #beta<-solve(t(x)%*%W%*%x)%*%t(x)%*%W%*%y;
-    rlm1<-rlm(y~xx,psi=psi.huber,weights=w,k=newc)
-    beta<-rlm1$coef
-    r2<-summary(lm(y~xx,weights=w))$r.squared
-    delta<-sqrt(sum((beta-beta0)^2));
-  }
-  list(esti=rlm1, Std.Error=summary(rlm1)$coef[,2], weights=w, tunning=newc, R2=r2)
-}
-
+##huber's function
+psi.Huber<-function(r,c)
+{
+  newr<-pmin(c,pmax(r,-c))
+  newr
+}
+
+##derivative of huber's function
+dpsi.Huber<-function(r,c)
+{
+  newr<-as.numeric(abs(r)<=c)
+  newr
+}
+
+
+##bisquare's function
+psi.Tukey<-function(r,c)
+{
+  r<-pmin(c,pmax(r,-c))
+  newr<-r*(1-(r/c)^2)^2
+  newr
+}
+
+##derivative of bisquare's function
+dpsi.Tukey<-function(r,c)
+{
+  r<-pmin(c,pmax(r,-c))
+  newr<-(1-(r/c)^2)*(1-5*(r/c)^2)
+  newr
+}
+
+##modified huber's function
+psi.Exp<-function(r,c)
+{
+  newr<-r*exp(-(r/c)^2)
+  newr
+}
+
+##derivative of modified huber's function
+dpsi.Exp<-function(r,c)
+{
+  newr<-(1-(2*r^2)/c^2)*exp(-(r/c)^2)
+  newr
+}
+
+
+##find the most efficient tuning constant
+eff<-function(r,method,plot)
+{
+  nn<-  length(r)
+  if (method =="Huber")
+  {
+    tau<-(1:30)/10
+    new<-unlist(lapply(tau,FUN=function(x){(sum(dpsi.Huber(r,x)))^2/(nn*sum(psi.Huber(r,x)^2))}))
+  }
+  else if (method=="Bisquare")
+  {
+    tau<-(15.48:59.48)/10
+    new<-unlist(lapply(tau,FUN=function(x){(sum(dpsi.Tukey(r,x)))^2/(nn*sum(psi.Tukey(r,x)^2))}))
+  }
+  else if (method=="Exponential")
+  {
+    tau<-(5:100)/10
+    new<-unlist(lapply(tau,FUN=function(x){(sum(dpsi.Exp(r,x)))^2/(nn*sum(psi.Exp(r,x)^2))}))
+  }
+
+  if (plot =="Y")
+  {
+    plot(new, type="o", xaxt = "n", xlab="Tunning parameter",
+        ylab="Efficiency",ylim=c(0,max(new)*1.3))
+    axis(1, at=seq_along(tau), labels=tau)
+    if (method =="Exponential") {
+      legend(length(tau)*0.8,max(new)*1.3,c("Exp"),bty = "n")
+    } else {
+      legend(length(tau)*0.8,max(new)*1.3,c(method),bty = "n")
+    }
+  }
+
+  xx<- order(new)
+  return(tau[xx[length(tau)]])
+
+}
+
+
+rho.h<-function(u,c){
+  phi<-0;
+  x1<-(abs(u)<=c)
+  phi<-x1* (u^2/2) + (1-x1)*(abs(u)*c-c^2/2)
+  phi
+}
+
+rho.b<-function(u,c){
+  phi<-0;
+  x1<-(abs(u)<=c)
+  phi<-x1*(1-(1-(u/c)^2)^3) + (1-x1)
+  phi
+}
+
+rho.e<-function(u,c){
+  1-exp(-(u/c)^2)
+}
+
+
+ESL_O<-function(x,xx,y,beta,newc,maxit=500, toler=1e-6)
+  {
+  it=0;delta=1;
+  while (delta>toler && it<maxit){
+    it<-it+1;
+    beta0<-beta;
+    w<-0;
+    e<-y-x%*%beta0;
+    S_n<-max(mad(e),1e-6)
+    u<-e/S_n
+    for(i in seq_along(u)){
+      w[i]<-exp(-(u[i]/newc)^2);
+    }
+
+    W<-diag(w);
+    #beta<-solve(t(x)%*%W%*%x)%*%t(x)%*%W%*%y;
+    rlm1<-rlm(y~xx,psi=psi.huber,weights=w,k=newc)
+    beta<-rlm1$coef
+    r2<-summary(lm(y~xx,weights=w))$r.squared
+    delta<-sqrt(sum((beta-beta0)^2));
+  }
+  list(esti=rlm1, Std.Error=summary(rlm1)$coef[,2], weights=w, tunning=newc, R2=r2)
+}
+
+
+chi<- function(obj){
+  pmin(obj^2/1.041^2,1)-0.5
+}
+
+create_lag <- function(vec, n.lag){
+  list.lag<-list()
+  for(i in 1: n.lag){
+    lagged.vec <- c(rep(NA,i),vec)[1:length(vec)]
+    list.lag[[paste0("lag", i)]] <- lagged.vec
+  }
+  return(list.lag)
+
+}
+