# R/iprls.R In expectreg: Expectile and Quantile Regression

#### Defines functions acvitpen.lsfitiprls

```iprls <- function(B,DD,yy,pp,lambda,smooth,nb,center,constmat,robust,adaptive=FALSE, LSMaxCores=1)
{
nterms = length(nb)
m = length(yy)
np = length(pp)

#  myapply <- lapply
#  if (.Platform\$OS.type == "unix" && require("multicore"))
#  {
#      if (!multicore:::isChild())
#      {
#          myapply <- mclapply
#      }
#  }

if(length(lambda) < nterms)
lala = rep(lambda[1],nterms)
else
lala = lambda

dummy.reg <- function(pp,lala,smooth,yy,B,DD,nb,nterms,center,constmat,robust,amplitude)
{
print(paste("tau:",pp,sep=" "))

if(smooth == "schall")
{
dc = 1
dw = 1
w <- rep(1,times=m)
it = 1
while((dc >= 0.01 || dw != 0) && it < 100)
{
aa <- itpen.lsfit(yy, B, pp, lala, DD, nb,constmat,robust,amplitude)
vector.a.ma.schall <- aa\$a

w0 <- w
l0 <- lala

for(i in 1:nterms)
{
partbasis = (sum(nb[0:(i-1)])+1):(sum(nb[0:i]))
if(center)
{
partB = B[,-1][,partbasis,drop=FALSE]
partDD = DD[,-1][-1,][,partbasis,drop=FALSE]
partaa = aa\$a[-1][partbasis]
}
else
{
partB = B[,partbasis,drop=FALSE]
partDD = DD[,partbasis,drop=FALSE]
partaa = aa\$a[partbasis]
}

if(any(partDD != 0))
{
v <- partDD %*% partaa
z <- aa\$fitted

w <- as.vector(aa\$weight)

H = solve(t(partB)%*%(w*partB) + lala[i]*t(partDD)%*%partDD)
##H = B%*%H%*%t(B)
#H = diag(H)*aa\$weight
H = apply(sqrt(w)*partB,1,function(x){t(x)%*%H%*%x})

sig2 <- sum(w * (yy - z) ^ 2,na.rm=TRUE) / (m - sum(aa\$diag.hat.ma,na.rm=TRUE))
tau2 <- sum(v ^ 2,na.rm=TRUE) / sum(H,na.rm=TRUE) + 1e-06

lala[i] = max(sig2 / tau2,1e-10,na.rm=TRUE)
}
}

dc <- max(abs(log10(l0 + 1e-6)-log10(lala + 1e-6)))
dw <- sum(w != w0,na.rm=TRUE)
it = it + 1
}
if(it == 100)
warning("Schall algorithm did not converge. Stopping after 100 iterations.")
}
else if(smooth == "acv")
{
acv.min = nlm(acv,p=lala,yy=yy,B=B,asymmetry=pp,DD=DD,nb=nb,constmat=constmat,robust=robust,amplitude=amplitude,iterlim=100)

aa <- itpen.lsfit(yy, B, pp, abs(acv.min\$estimate), DD, nb,constmat,robust,amplitude)
vector.a.ma.schall <- aa\$a
lala <- abs(acv.min\$estimate)
}
else
{
aa <- itpen.lsfit(yy, B, pp, lala, DD, nb,constmat,robust,amplitude)
vector.a.ma.schall <- aa\$a
}

list(vector.a.ma.schall,lala,aa\$weight,aa\$diag.hat.ma)
}

amplitude = NULL

{
for(i in 1:2)
{
#acv.min = nlm(acv,p=lala,yy=yy,B=B,asymmetry=0.5,DD=DD,nb=nb,constmat=constmat,robust=robust,amplitude=amplitude,iterlim=100)
#lala <- abs(acv.min\$estimate)
#aa <- itpen.lsfit(yy, B, 0.5, abs(acv.min\$estimate), DD, nb,constmat,robust,amplitude=amplitude)

coefa = dummy.reg(0.5,lala,smooth,yy,B,DD,nb,nterms,center,constmat,robust,amplitude)[[1]]

res = abs(yy - B %*% coefa)

coefb = dummy.reg(0.5,lala,smooth,res,B,DD,nb,nterms,center,constmat,robust,amplitude=NULL)[[1]]

amplitude = B %*% coefb

#resamp = abs(res - amplitude)

#coefc = dummy.reg(0.5,lala,smooth,resamp,B,DD,nb,nterms,center,constmat,robust,amplitude=NULL)[[1]]

#amp2 = B %*% coefc

#coefd = dummy.reg(0.5,lala,smooth,res,B,DD,nb,nterms,center,constmat,robust,amplitude=amp2)[[1]]

#amplitude = B %*% coefd

#res = abs(yy - B %*% aa\$a)

#acv.min = nlm(acv,p=lala,yy=res,B=B,asymmetry=0.5,DD=DD,nb=nb,constmat=constmat,robust=robust,amplitude=NULL,iterlim=100)

#aa <- itpen.lsfit(res, B, 0.5, abs(acv.min\$estimate), DD, nb,constmat,robust,NULL)

#amplitude = B %*% aa\$a
}
}

if (.Platform\$OS.type == "unix")
coef.vector = mclapply(pp,function(pp) dummy.reg(pp,lala,smooth,yy,B,DD,nb,nterms,center,constmat,robust,amplitude),mc.cores = max(1,min(detectCores()-1, LSMaxCores)))
else if (.Platform\$OS.type == "windows")
coef.vector = mclapply(pp,function(pp) dummy.reg(pp,lala,smooth,yy,B,DD,nb,nterms,center,constmat,robust,amplitude),mc.cores = 1)

lala <- matrix(lambda, nrow=nterms, ncol=np)
vector.a.ma.schall <- matrix(NA, nrow=sum(nb)+(1*center),ncol=np)
ww = matrix(NA,nrow=m,ncol=np)
diag.hat = matrix(NA,nrow=m,ncol=np)

for(i in 1:np)
{
vector.a.ma.schall[,i] = coef.vector[[i]][[1]]
lala[,i] = coef.vector[[i]][[2]]
ww[,i] = coef.vector[[i]][[3]]
diag.hat[,i] = coef.vector[[i]][[4]]
}

return(list(vector.a.ma.schall,lala,amplitude,ww,diag.hat))
}

itpen.lsfit <- function(y, B, p, lambda, DD,nb,constmat,robust,amplitude=NULL)
###### Asymmetric regression with difference penalty
# parameters:
# y - response variable
# B - B-spline basis
# p - asymmetry parameter
# lambda - smoothing parameter
# DD - difference matrix
#
# needed: ncol(B) = ncol(DD)
{
#w1 <- 0 * y + 0.5
n <- ncol(B)
lambda = c(rep(0,times=n - sum(nb)),rep(lambda,times=nb))

P <- sqrt(lambda) * DD
augm <- rep(0, nrow(P))
conpen <- rep(0,nrow(constmat))

diffcon = -1
it.con = 1

my.psi.q<-function(u,q,robust,amplitude)
{
# function for the Huber 2 influence function at quantile q
if(all(is.null(amplitude)))
{
s<-median(abs(u))/0.6745
c<-robust #1.345
}
else
{
s =1/0.6745
c<-abs(amplitude) * robust #1.345
}

u<-u/s
abs(-(1-q)*c*(u< -c)+(1-q)*u*(u<0 & u>=-c)+q*u*(u>=0 & u<c)+q*c*(u>=c))
}

resi = B %*% rep(1,ncol(B))

w1 = my.psi.q(y-resi,p,robust,amplitude)/(abs(y-resi) + 1e-5*(y-resi == 0))

while(any(diffcon < -1e-5) && it.con < 20)
{
it = 1
dw1 = 1

while(dw1 != 0 && it < 100)
{
model <- lsfit(x=rbind(B,P,constmat*conpen), y=c(y, augm,0*conpen), wt=c(w1,(augm+1),1*(conpen>0)), intercept=FALSE)
a1 <- model\$coefficients
z1 <- B %*%a1
w01 <- w1
#w1 <- as.vector(ifelse(y > z1, p, 1 - p))

w1 = my.psi.q(y-z1,p,robust,amplitude)/(abs(y-z1) + 1e-5*(y-z1 == 0))

dw1 <- sum(w1 != w01,na.rm=TRUE)
it = it + 1
}
diffcon =  constmat %*% a1

if(any(diffcon < 0))
{
wc = which(diffcon < 0)
conpen[wc] = conpen[wc] + 100000
}
it.con = it.con + 1
}
diag.hat.ma1 <- hat(model\$qr)[1:length(y)]

if(it == 50)
warning("IWLS weights did not converge after 100 iterations.")

list(a=a1, diag.hat.ma=diag.hat.ma1, weight = w1, fitted = z1)
}

acv <- function(penalty,yy,B,asymmetry,DD,nb,constmat,robust,amplitude=NULL)
# asymmetric cross validation
# computes the acv score for the smoothing of the regression
# score has to be minimized dependant on parameter "penalty"
# therefore a grid search can be applied to this function
# parameters:
# penalty - smoothing parameter lambda
# yy - vector of responses
# B - basis for the approximation
# p - quantile
# DD - penalization matrix
{
aa <- itpen.lsfit(yy, B, asymmetry, abs(penalty), DD, nb, constmat,robust,amplitude)

#H = solve(t(B)%*%(aa\$weight*B) + penalty*t(DD)%*%DD)
##H = B%*%H%*%t(B)
##H = diag(H)*aa\$weight
#H = apply(sqrt(aa\$weight)*B,1,function(x){t(x)%*%H%*%x})

#H = diag(diag(aa\$diag.hat.ma))

score = aa\$weight*(yy - B%*%aa\$a)^2/(1-aa\$diag.hat.ma)^2

mean(score[which(is.finite(score))],na.rm=TRUE)
}
```

## Try the expectreg package in your browser

Any scripts or data that you put into this service are public.

expectreg documentation built on March 18, 2022, 5:57 p.m.