# R/cde.bandwidths.R In hdrcde: Highest Density Regions and Conditional Density Estimation

#### Documented in cde.bandwidths

#' Bandwidth calculation for conditional density estimation
#'
#' Calculates bandwidths for kernel conditional density estimates. Methods
#' described in Bashtannyk and Hyndman (2001) and Hyndman and Yao (2002).
#'
#' Details of the various algorithms are in Bashtannyk and Hyndman (2001) and
#' Hyndman and Yao (2002).
#'
#' @param x Numerical vector: the conditioning variable.
#' @param y Numerical vector: the response variable.
#' @param deg Degree of local polynomial used in estimation.
#' @param link Link function used in estimation. Default "identity". The other
#' possibility is "log" which is recommended if degree > 0.
#' @param method
#'   \describe{
#'     \item{method = 1:}{Hyndman-Yao algorithm if deg>0; Bashtannyk-Hyndman algorithm if deg=0;}
#'     \item{method = 2:}{Normal reference rules;}
#'     \item{method = 3:}{Bashtannyk-Hyndman regression method if deg=0;}
#'     \item{method = 4:}{Bashtannyk-Hyndman bootstrap method if deg=0.}
#' }
#' @param y.margin Values in y-space on which conditional density is
#' calculated. If not specified, an equi-spaced grid of 50 values over the
#' range of y is used.
#' @param passes Number of passes through Bashtannyk-Hyndman algorithm.
#' @param ngrid Number of values of smoothing parameter in grid.
#' @param min.a Smallest value of a to consider if method=1.
#' @param ny Number of values to use for y margin if \code{y.margin} is
#' missing.
#' @param use.sample Used when regression method (3) is chosen.
#' @param GCV Generalized cross-validation. Used only if method=1 and deg>0. If
#' GCV=FALSE, method=1 and deg=0, then the AIC is used instead.  The argument
#' is ignored if deg=0 or method>1.
#' @param b Value of b can be specified only if method=1 and deg>0. For deg=0
#' or method>1, this argument is ignored.
#' @param \dots Other arguments control details for individual methods.
#' @return \item{a}{Window width in \code{x} direction.} \item{b}{Window width
#' in \code{y} direction.}
#' @author Rob J Hyndman
#' @references Hyndman, R.J., Bashtannyk, D.M. and Grunwald, G.K. (1996)
#' "Estimating and visualizing conditional densities". \emph{Journal of
#' Computational and Graphical Statistics}, \bold{5}, 315-336.
#'
#' Bashtannyk, D.M., and Hyndman, R.J. (2001) "Bandwidth selection for kernel
#' conditional density estimation". \emph{Computational statistics and data
#' analysis}, \bold{36}(3), 279-298.
#'
#' Hyndman, R.J. and Yao, Q. (2002) "Nonparametric estimation and symmetry
#' tests for conditional density functions". \emph{Journal of Nonparametric
#' Statistics}, \bold{14}(3), 259-278.
#' @keywords smooth distribution
#' @examples
#' bands <- cde.bandwidths(faithful$waiting,faithful$eruptions,method=2)
#' plot(cde(faithful$waiting,faithful$eruptions,a=bands$a,b=bands$b))
#' @export cde.bandwidths
ngrid=8,min.a=NULL,ny=25,use.sample=FALSE,GCV=TRUE,b=NULL,...)
{
if(deg>0 & method>2)
stop("Method unavailable for degree > 0")

if(missing(y.margin))
y.margin <- seq(min(y),max(y),l=ny)
else
ny <- length(y.margin)

if(deg==0)
return(cde.bandwidths0(x,y,method,y.margin,
ny=ny,passes=passes,use.sample=use.sample,ngrid=ngrid,...))
else if(method == 2)
else # method = 1
{
if(!is.null(b))
bands$b <- b print(paste("Initial values: a=",round(bands$a,3),"  b=",round(bands$b,3)),sep="") if(is.null(min.a)) min.a <- 0.2*bands$a
a.grid <- bands$a*seq(min.a/bands$a,2,l=ngrid)
regout <- cdeband.Mbh(x,y,a.grid,bands$b,y.margin=y.margin,use.sample=use.sample,deg=deg,link=link,GCV=GCV,passes=1,ngrid=ngrid,usequad=FALSE,...) return(list(a=regout$a,b=bands$b,a.grid=regout$a.grid,q=regout$q)) } } "cde.bandwidths0" <- function(x,y,method=1,y.margin,sdlinear=FALSE,xden="normal",penalty=4, ngrid=8,modified=FALSE,k=3, m=25,nx=30,ny=25,passes=2,tol=0.99999,use.sample=FALSE,usequad=TRUE) { if(missing(y.margin)) y.margin <- seq(min(y),max(y),l=ny) else ny <- length(y.margin) if(method==2) return(cdeband.rules0(x,y,sdlinear,xden,k,modified)) else if(method==3) return(cdeband.regress(x,y,use.sample=use.sample,y.margin=y.margin, penalty=penalty,tol=tol,usequad=usequad)) else if(method==4) return(cdeband.bootstrap(x,y,m=m,nx=nx,y.margin=y.margin)) else { firstbands <- bands <- cdeband.rules0(x,y,sdlinear=sdlinear,xden=xden,modified=modified,k=k) print(paste("Initial values: a=",round(bands$a,3),"  b=",round(bands$b,3)),sep="") for(i in 1:passes) { if(i==1) { a.grid <- bands$a*seq(0.1,3,l=ngrid)
b.grid <- bands$b*seq(0.1,3,l=ngrid) } else # Choose values around previous minimums { na <- length(regout$a.grid)
nb <- length(bootout$b.grid) adiff <- abs(bands$a-regout$a.grid) idx <- (1:na)[adiff == min(adiff)] a.grid <- seq(regout$a.grid[max(idx-2,1)],regout$a.grid[min(idx+2,na)],l=ngrid) bdiff <- abs(bands$b-bootout$b.grid) idx <- (1:nb)[bdiff == min(bdiff)] b.grid <- seq(bootout$b.grid[max(idx-2,1)],bootout$b.grid[min(idx+2,nb)],l=min(8,ngrid)) } regout <- cdeband.regress(x,y,a.grid,bands$b,y.margin=y.margin,tol=tol,use.sample=use.sample,usequad=usequad)
regna <- is.na(regout$a) if(regna) bands$a <- firstbands$a else bands$a <- regout$a regna <- is.na(regout$a)
bootout <- cdeband.bootstrap(x,y,bands$a,b.grid,m,nx,y.margin=y.margin) bands$b <- bootout$b } return(list(a=bands$a,b=bands$b,a.grid=regout$a.grid,b.grid=bootout$b.grid,q=regout$q,imse=bootout$imse,regna=regna)) } } "cdeband.rules0" <- function(x,y,sdlinear=FALSE,xden="normal",k=3,modified=FALSE) { if(xden=="normal" & sdlinear) stop("Option not available") if(modified & xden!="normal") stop("Modified estimator requires a normal marginal density") Rk <- 0.5/sqrt(pi) sk <- 1 n <- length(x) junk <- Lsquare(x,y,sdlinear) p <- junk$pc
d <- junk$dx pl <- NA qx <- NA if(xden=="uniform") { u <- max(x) l <- min(x) if(!sdlinear) { a <- ((4*sqrt(pi)*Rk^2*(u-l)*abs(p)^5)/(3*n*sk^4*abs(d)^5))^(1/6) b <- abs(d)*a } else { pl <- junk$pl
qx <- junk$qx z <- ((pl+qx*u)^4 - (pl+qx*l)^4)/((pl+qx*u)^4 *(pl+qx*l)^4) w <- 19*qx^4 + 4*d^4 + 28*d^2*qx^2 a <- ((2^(7.5)*sqrt(pi)*Rk^2*(u-l)^2*qx)/(3*n*sk^4*z*w^(0.75)*(sqrt(w)+2*d^2-3*qx^2)))^(1/6) b <- a*w^(0.25)/sqrt(2) } } else if(xden=="normal") { erf <- 0.954499736 sm <- sqrt(var(x)) if(!modified) { v <- 3*erf*d^2*sm^3*pi - 8*sqrt(2*pi)*p^2*k*exp((-k^2)/2) + 8*pi*p^2*erf anum <- ((16*Rk^2*k*(pi)^(1.25)*p^5*sm^(2.5)) / (n*abs(d)^(5/2)*sk^4))^(1/6) adem <- (((v^5)/(3*(pi)^2*sm^4*erf))^(1/4) + 3*d*((v*erf^(1/3))/3)^(3/4))^(1/6) a <- anum/adem b <- ((d^2*v)/(3*pi*sm*erf))^(1/4) * a } else { u <- 3*d^2*sm^2 + 4*p^2 anum <- 16*pi*Rk^2*sm^2.5*p^5 adem <- n * sk^4 * abs(d)^2.5 * ((u^5/(3*sm^4))^0.25 + (3*d^4*u^3)^0.25) a <- (anum/adem)^(1/6) b <- a *(d^2*u/3/sm^2)^0.25 } } else stop("Unknown distribution") return(list(a=a,b=b,p=p,d=d,pl=pl,q=qx)) } "cdeband.rules" <- function(x,y,deg,link,mean.order,...) { if(missing(mean.order)) { if(deg==1) mean.order <- 2 else mean.order <- 1 } if(mean.order<1 | mean.order>2 | deg<0 | deg>2) stop("Not implemented") if(deg==2 & mean.order!=1) stop("Not implemented") if(deg==0) { if(mean.order!=1 | link!="identity") stop("Not implemented") else return(cdeband.rules0(x,y,...)) } n <- length(x) if(mean.order==1) { model <- lm(y~x) d1 <- abs(model$coef[2])
d2 <- 0
}
else
{
mu <- mean(x)
model <- lm(y~ I(x-mu) + I((x-mu)^2))
d1 <- abs(model$coef[2]) d2 <- model$coef[3]
}
res <- residuals(model)
sigma <- sqrt(sum(res^2)/(length(x)-2))

v <- sqrt(var(x))
gamma <- 3/(64*pi*v*sigma^5)
if(deg==1)
{
eta <- 0.5
tau2 <- 0.25/pi
r1 <- 2
{
alpha <- (3*d1^4+36*d2^2*v^2*(d1^2+v^2*d2^2)+8*d2^2*sigma^2)/(64*pi*sigma^5*v)
beta <- 3*(d1^2+2*d2^2*v^2)/(32*pi*sigma^5*v)
}
else
{
alpha <- (d1^4+12*d2^2*v^2*(d1^2+v^2*d2^2)+2*d2^2*sigma^2)/(16*pi*sigma^5*v)
beta <- (d1^2+2*d2^2*v^2)/(16*pi*sigma^5*v)
}
cc <- (alpha/gamma)^0.25
}
else
{
if(d2!=0)
stop("Not implemented")
{
alpha <- 105*d1^8/(256*pi*sigma^9*v)
beta <- 15*d1^4/(64*pi*sigma^7*v)
cc <- d1^2*sqrt(sqrt(305)+5)/(2*sigma)
eta <- -0.5
tau2 <- 27/(64*pi)
r1 <- 4
}
else
stop("Not implemented")
}
power <- 2/(5*r1+2)
hstar <- (tau2/(n*cc*r1*(2*alpha + beta*cc^2)))^power
bstar <- cc*(hstar)^(r1/2)

return(list(a=hstar,b=bstar,sigma=sigma,d1=d1,d2=d2,v=v))
}

cdeband.Mbh <- function(x,y,a.grid,b,ny=50,use.sample=FALSE,nx=100,y.margin,passes=2,deg=deg,
{
## Initialization

if(missing(a.grid))
a.grid <- bands$a*seq(0.1,5,l=ngrid) if(missing(b)) b <- bands$b
a.grid <- a.grid[a.grid>0]
na <- length(a.grid)

n <- length(x)
if(use.sample & n > nx)
idx <- sample(1:n,nx,replace=F)
else
idx <- 1:n
xx <- x[idx]
yy <- y[idx]
if(missing(y.margin))
y.margin <- seq(min(y),max(y),l=ny)

## First pass

if(is.na(first$a)) # Expand grid { a.grid <- seq(a.grid[1]/20,1.1*a.grid[na],l=ngrid) firstb <- CDEband.Mbh(xx,yy,a.grid,b,y.margin,deg,link,GCV=GCV) first$a.grid <- c(first$a.grid,firstb$a.grid)
idx <- order(first$a.grid) first$q <- c(first$q,firstb$q)[idx]
first$a.grid <- first$a.grid[idx]
first$a <- firstb$a
}

## Second pass

if(passes>1 & !is.na(first$a)) { na <- length(first$a.grid)
idx <- c(1:na)[first$a.grid==first$a]
if(na>4)
idx <- idx + (-2:2)
else
idx <- idx + (-1:1)
idx <- idx[idx>=1 & idx<=na]
newa.grid <- seq(0.95*first$a.grid[idx[1]],1.05*first$a.grid[idx[length(idx)]],l=ngrid)
fulla.grid <- c(first$a.grid,second$a.grid)
idx <- order(fulla.grid)
fullq <- c(first$q,second$q)[idx]
fulla.grid <- fulla.grid[idx]
a <- second$a } else { a <- first$a
fulla.grid <- first$a.grid fullq <- first$q
}

{
na <- length(fulla.grid)
idx <- (c(1:na)[fulla.grid==a])[1]
if(na>4)
idx <- idx + (-2:2)
else
idx <- idx + (-1:1)
idx <- idx[idx>=1 & idx<=na & fullq[idx] != Inf]
if(length(idx)>3)
{
fit <- lm(q ~ a+I(a^2),data=data.frame(q=fullq[idx],a=fulla.grid[idx]))
a <- -0.5*fit$coef[2]/fit$coef[3]
}
}
else
{
a <- fulla.grid[fullq==min(fullq,na.rm=T)]
a <- a[!is.na(a)]
}

return(list(a=a,b=b,a.grid=fulla.grid,q=fullq))
}

{
na <- length(a.grid)
q <- numeric(na)
cat("\n Trying a=")
for(i in 1:na)
{
cat(round(a.grid[i],3)," ")
if(GCV)
q[i] <- junk$GCV else q[i] <- junk$AIC
}
cat("\n")
idx <- (1:na)[q==min(q,na.rm=TRUE)]
idx <- idx[!is.na(idx)]
if(idx==1 | idx==na)
{
warning("No minimum found")
a <- NA
}
else
a <- a.grid[idx]
return(list(a=a,a.grid=a.grid,q=q))
}

"cdeband.bootstrap" <- function(x,y,a.grid,b.grid,m=25,nx=30,ny=25,y.margin)
{
## Fit parametric model and calculate parametric conditional density
df <- data.frame(x=x,y=y)
aic <- numeric(4)
n <- length(x)
fit1 <- lm(y ~ 1,data=df)
aic[1] <- deviance(fit1) + 2*(n-fit1$df.residual)*deviance(fit1)/fit1$df.resid
fit2 <- lm(y ~ x,data=df)
aic[2] <- deviance(fit2) + 2*(n-fit2$df.residual)*deviance(fit2)/fit2$df.resid
fit3 <- lm(y ~ x + I(x^2),data=df)
aic[3] <- deviance(fit3) + 2*(n-fit3$df.residual)*deviance(fit3)/fit3$df.resid
fit4 <- lm(y~ x+I(x^2)+I(x^3),data=df)
aic[4] <- deviance(fit4) + 2*(n-fit4$df.residual)*deviance(fit4)/fit4$df.resid
fit <- switch((1:4)[aic==min(aic)],fit1,fit2,fit3,fit4)
rse <- summary(fit)$sigma fits <- fitted(fit) n <- length(x) if(missing(y.margin)) y.margin <- seq(min(y),max(y),l=ny) else ny <- length(y.margin) if(length(x)<nx) { x.margin <- sort(x) nx <- length(x) } else x.margin <- sort(sample(x,nx)) fit.grid <- approx(x,fits,xout=x.margin,rule=2)$y  # Quicker than predict and works when constant model used.
truecde <- list(x=x.margin,y=y.margin,z=matrix(NA,nx,ny))
for(i in 1:nx)
truecde$z[i,] <- dnorm(y.margin,fit.grid[i],rse) ## Get bandwidth grids bands <- cdeband.rules0(x,y) if(missing(a.grid)) a.grid <- bands$a*seq(0.4,2,0.2)
if(missing(b.grid))
b.grid <- bands$b*seq(0.4,2,0.2) ## Simulate bootstrap samples m <- max(m,5) bootsamples <- matrix(NA,nrow=n,ncol=m) for(i in 1:m) bootsamples[,i] <- fits + rnorm(n,0,rse) ## Calculate IMSE diff.cde <- numeric(m) na <- length(a.grid) nb <- length(b.grid) delta <- y.margin[2]-y.margin[1] imse <- matrix(NA,na,nb) cat("\n Trying (a,b)=") ## First pass for(i in 1:na) { a <- a.grid[i] for(j in 1:nb) { b <- b.grid[j] cat("(",round(a,3),",",round(b,3),") ",sep="") for(k in 1:3) { bootcde <- cde(x,bootsamples[,k],deg=0,link="identity",a=a,b=b,x.margin=x.margin,y.margin=y.margin) diff.cde[k] <- sum((bootcde$z - truecde$z)^2)/(m*nx) } imse[i,j] <- mean(diff.cde[1:3]) } } cat("\n") ## Second pass near minimum idx <- (imse==min(imse)) idx.a <- max(3,(1:na)[apply(idx,1,sum)==1]) + c(-2:2) idx.b <- max(3,(1:nb)[apply(idx,2,sum)==1]) + c(-2:2) idx.a <- idx.a[idx.a>0 & idx.a <=na] idx.b <- idx.b[idx.b>0 & idx.b <=nb] for(i in idx.a) { a <- a.grid[i] for(j in idx.b) { b <- b.grid[j] cat("(",round(a,3),",",round(b,3),") ",sep="") for(k in 4:m) { bootcde <- cde(x,bootsamples[,k],deg=0,link="identity",a=a,b=b,x.margin=x.margin,y.margin=y.margin) diff.cde[k] <- sum((bootcde$z - truecde$z)^2)/(m*nx) } imse[i,j] <- (3*imse[i,j] + (m-3)*mean(diff.cde[4:m]))/m } } cat("\n") # Find minimum on grid idx <- (imse==min(imse)) idx.a <- (1:na)[apply(idx,1,sum)>0] idx.b <- (1:nb)[apply(idx,2,sum)>0] best.a <- a.grid[idx.a[1]] best.b <- b.grid[idx.b[1]] bestmse <- imse[idx.a[1],idx.b[1]] ## Fit quadratic surface to points near minimum and find minimum of surface na <- length(idx.a) nb <- length(idx.b) if(na>2 & nb>2) { fit <- lm(imse ~ a + b + I(a^2) + I(b^2) + I(a*b), data=data.frame(imse=c(imse[idx.a,idx.b]),a=rep(a.grid[idx.a],nb),b=rep(b.grid[idx.b],rep(na,nb)))) best.b <- (fit$coef[2]*fit$coef[6] - 2*fit$coef[3]*fit$coef[4])/(4*fit$coef[4]*fit$coef[5] - fit$coef[6]*fit$coef[6]) best.a <- (best.b*fit$coef[6] + fit$coef[2])/(-2*fit$coef[4])
}
else if(nb>2) # na=1 or 2
{
for(i in 1:na)
{
fit <- lm(imse ~ b+I(b^2),data=data.frame(imse=imse[idx.a[i],idx.b],b=b.grid[idx.b]))
bestb <- -0.5*fit$coef[2]/fit$coef[3]
minimse <- predict(fit,newdata=data.frame(b=bestb))
if(minimse<bestmse & bestb >= min(b.grid) & bestb <= max(b.grid))
{
best.b <- bestb
best.a <- a.grid[idx.a[i]]
bestmse <- minimse
}
}
}
else if(na>2) # nb=1 or 2
{
for(i in 1:nb)
{
fit <- lm(imse ~ a+I(a^2),data=data.frame(imse=imse[idx.a,idx.b[i]],a=a.grid[idx.a]))
besta <- -0.5*fit$coef[2]/fit$coef[3]
minimse <- predict(fit,newdata=data.frame(a=besta))
if(minimse<bestmse  & besta >= min(a.grid) & besta <= max(a.grid))
{
best.a <- besta
best.b <- b.grid[idx.b[i]]
bestmse <- minimse
}
}
}

return(list(a=best.a,b=best.b,a.grid=a.grid,b.grid=b.grid,imse=imse*delta))
}

## This function called by cdeband.regress.

CDEband.regress <- function(x, y, a.grid, b, y.margin, penalty=4, tol=0.999)
{
# Calculate v
n <- length(x)
ny <- length(y.margin)
if(ny>1)
delta <- y.margin[2]-y.margin[1]
else
delta <- 1
v <- Kernel(y,y.margin,b)
na <- length(a.grid)
q <- numeric(na)
cat("\n Trying a=")
for(i in 1:na)
{
cat(round(a.grid[i],3)," ")
w <- Kernel(x,x,a.grid[i])
row.sum <- apply(w,1,sum)
w <- w  / matrix(rep(row.sum,n),ncol=n)
junk <- (v - (t(w) %*% v))^2
diag.w <- diag(w)
pen <- switch(penalty,1+2*diag.w,1/(1-diag.w)^2,exp(2*diag.w),(1+diag.w)/(1-diag.w),1/(1-2*diag.w))
tmp <- apply(junk,1,sum)
q[i] <- mean(tmp*pen,na.rm=TRUE)
#        if(sum(diag.w>tol)>0)
#            q[i] <- Inf
#        else
#        {
#            pen <- switch(penalty,1+2*diag.w,1/(1-diag.w)^2,exp(2*diag.w),(1+diag.w)/(1-diag.w),1/(1-2*diag.w))
#            q[i] <- mean(apply(junk,1,sum)*pen)
#        }
}
cat("\n")
idx <- (1:na)[q==min(q)]
if(idx==1 | idx==na)
{
#        warning("No minimum found")
a <- NA
}
else
a <- a.grid[idx]
#    browser()
return(list(a=a,a.grid=a.grid,q=q*delta))
}

{
## Initialization

bands <- cdeband.rules0(x,y)
if(missing(a.grid))
a.grid <- bands$a*seq(0.01,3,l=10) if(missing(b)) b <- bands$b
a.grid <- a.grid[a.grid>0]
na <- length(a.grid)

n <- length(x)
if(use.sample & n > nx)
idx <- sample(1:n,nx,replace=FALSE)
else
idx <- 1:n
xx <- x[idx]
yy <- y[idx]
if(missing(y.margin))
y.margin <- seq(min(y),max(y),l=ny)

## First pass

first <- CDEband.regress(xx,yy,a.grid,b,y.margin,penalty=penalty,tol=tol)
if(is.na(first$a)) # Expand grid { a.grid <- seq(a.grid[1]/20,1.1*a.grid[na],l=50) firstb <- CDEband.regress(xx,yy,a.grid,b,y.margin,penalty=penalty,tol=tol) first$a.grid <- c(first$a.grid,firstb$a.grid)
idx <- order(first$a.grid) first$q <- c(first$q,firstb$q)[idx]
first$a.grid <- first$a.grid[idx]
first$a <- firstb$a
}

## Second pass

if(passes>1 & !is.na(first$a)) { na <- length(first$a.grid)
idx <- c(1:na)[first$a.grid==first$a]
if(na>4)
idx <- idx + (-2:2)
else
idx <- idx + (-1:1)
idx <- idx[idx>=1 & idx<=na]
newa.grid <- seq(0.95*first$a.grid[idx[1]],1.05*first$a.grid[idx[length(idx)]],l=length(a.grid))
second <- CDEband.regress(xx,yy,newa.grid,b,y.margin,tol=tol,penalty=penalty)
fulla.grid <- c(first$a.grid,second$a.grid)
idx <- order(fulla.grid)
fullq <- c(first$q,second$q)[idx]
fulla.grid <- fulla.grid[idx]
a <- second$a } else { a <- first$a
fulla.grid <- first$a.grid fullq <- first$q
}

{
na <- length(fulla.grid)
idx <- (c(1:na)[fulla.grid==a])[1]
if(na>4)
idx <- idx + (-2:2)
else
idx <- idx + (-1:1)
idx <- idx[idx>=1 & idx<=na & fullq[idx] != Inf]
if(length(idx)>3)
{
fit <- lm(q ~ a+I(a^2),data=data.frame(q=fullq[idx],a=fulla.grid[idx]))
besta <- -0.5*fit$coef[2]/fit$coef[3]
if(besta >0)
a <- besta
}
}

return(list(a=a,b=b,a.grid=fulla.grid,q=fullq))
}

"Lsquare" <- function(x,y,sdlinear=TRUE)
{
# Fits linear model.
# If sdlinear=TRUE, allows heteroskedastic errors.

x <- c(x)
y <- c(y)
model <- lm(y~x)
res <- residuals(model)
pc <- sqrt(sum(res^2)/(length(x)-2))
if(!sdlinear)
return(list(pc=pc,dx=model$coef[2])) ## Initial estimate of heteroskedasticity model.res <- lm(abs(res)~x) ## Recalculate linear fit weights <- fitted(model.res)^(-2) model <- lm(y~x,weights=weights) ## Recalculate heteroskedasticity res <- residuals(model) f1 <- function(co,res,x) { sum((res * res - (co[1] + co[2] * x)^2)^2) } junk <- nlm(f1,model.res$coef,res=res,x=x)
return(list(pc=pc,pl=junk$estimate[1],qx=junk$estimate[2],dx=model\$coef[2]))
}


## Try the hdrcde package in your browser

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

hdrcde documentation built on Jan. 18, 2021, 9:05 a.m.