Nothing
## VAS & MEM #########################################################
## non-parametric regression for weighted data -- an alternative way
## to handle measurement errors. Remarks: (1) a single valued 'bw' is
## used to avoid heavy computational burdens. (2) LSCV is used to
## select the optimal bandwidth.
.vasbin2 <- function(x,y){
x <- as.matrix(x)
y <- as.numeric(y)
nx <- apply(!is.na(x),1,sum)
wx <- matrix(1/nx, nrow=nrow(x),ncol=ncol(x))
y0 <- matrix(y,nrow=length(y), ncol=ncol(x))
wx[is.na(x)] <- 0
x <- as.numeric(x)
w <- as.numeric(wx)
y1 <- as.numeric(y0)
sele <- w > 0
x <- x[sele]
w <- w[sele]
y <- y1[sele]
list(x=x,w=w,y=y)
}
VAS.npr <- function(y,x,bw,s.x,from,to,gridsize,conf.level=0.95)
{
## y: dependent variable
## x: predictor
## wlpsmooth if x is a matrix
## lpsmooth if x is a vector
if(missing(s.x)){
s.x <- 0
}else{
stopifnot(is.numeric(s.x))
stopifnot(length(s.x)==1)
stopifnot(s.x>=0)
}
y <- as.matrix(y)
x <- as.matrix(x)
stopifnot(nrow(x)==nrow(y))
stopifnot(ncol(y)==1)
if(ncol(x)==1){
y <- as.numeric(y)
x <- as.numeric(x)
sele <- is.na(x) | is.na(y)
if(any(sele)){
if(sum(!sele)<5)
stop("too many missing values")
x <- x[!sele]
y <- y[!sele]
}
out <- lpsmooth(y=y,x=x,bw=bw,
from=from,to=to,gridsize=gridsize,
conf.level=conf.level)
}else{
res <- .vasbin2(x,y)
out <- wlpsmooth(y=res$y,x=res$x,w=res$w,s.x=s.x,bw=bw,
from=from,to=to,gridsize=gridsize,
conf.level=conf.level)
}
out
}
## VAS.pdf estimate the density function using nonparametric
## regression (local polynomial with df=1) via binning. 'x' will be
## binned automatically. An existing R function lpsommoth will be
## called with the following settings: (1) lscv=TRUE (to refine the
## bandwidth estimate); (2) adaptive=FALSE (enhance the speed).
VAS.pdf <- function(x,y,bw,s.x, gridsize = 401L,type="absolute",
range.x, truncate = FALSE){
type <- match.arg(tolower(type),
c("absolute","percent","relative"))
if(missing(s.x)){
s.x <- 0
}else{
stopifnot(is.numeric(s.x))
stopifnot(length(s.x)==1)
stopifnot(s.x>0)
}
wx <- .vasbin(x=x,y=y,type=type)
x <- wx$x
w <- wx$w
n <- sum(w)
w <- w/n
mu <- sum(w*x)
vx <- sum(w*(x-mu)^2)*n/(n-1)
## density estimation with RR
if(wx$ME){
if(s.x > 0){
out <- wdekde(x=wx$x,w=wx$w,bandwidth=bw,s.x=s.x)
}else{
out <- wkde(x=wx$x,w=wx$w,bandwidth=bw)
}
}else{
out <- density(wx$x)
}
x0 <- out$x; x1 <- x0
f0 <- out$y; f1 <- f0
if(!missing(range.x)){
stopifnot(is.numeric(range.x))
stopifnot(length(range.x)==2)
a <- range.x[[1]]
b <- range.x[[2]]
## apply RR to left end
sele <- x0 < a
if(any(sele)){
k <- sum(sele)
x1[sele] <- NA
f1[(k+1):(2*k)] <- f1[(k+1):(2*k)]+rev(f1[sele])
f1[sele] <- 0
}
## apply RR to right end
sele <- x0 > b
if(any(sele)){
k <- sum(sele)
n <- length(x1)
x1[sele] <- NA
f1[(n-2*k+1):(n-k)] <- f1[(n-2*k+1):(n-k)]+rev(f1[sele])
f1[sele] <- 0
}
sele <- f1 > 0
x1 <- x1[sele]
f1 <- f1[sele]
}
list(x=x1,y=f1,x0=x0,y0=f0,
data=wx$x, wt=wx$w,
mean=mu,s=sqrt(vx),bw=out$bw)
}
.pc <- function(x,y) (x-y)/x
.vaschange <- function(x,nx,ny,type){
x0 <- x[1:nx]
x1 <- x[(nx+1):(nx+ny)]
if(type=="absolute"){
out <- as.numeric(outer(x0,x1,"-"))
}else{
out <- as.numeric(outer(x0,x1,FUN=.pc))
}
out
}
.vasbin <- function(x,y,type){
x <- as.matrix(x)
if(ncol(x)==1) ME <- FALSE
else ME <- TRUE
if(missing(y)){
nx <- apply(!is.na(x),1,sum)
wx <- matrix(1/nx, nrow=nrow(x),ncol=ncol(x))
wx[is.na(x)] <- 0
x <- as.numeric(x)
w <- as.numeric(wx)
sele <- w > 0
x <- x[sele]
w <- w[sele]
}else{
y <- as.matrix(y)
if(ncol(x)==1&&ncol(y)==1) ME <- FALSE
else ME <- TRUE
stopifnot(nrow(x)==nrow(y))
nx <- ncol(x); ny <- ncol(y)
xy <- cbind(x,y)
x <- apply(xy, 1, .vaschange,nx=nx,ny=ny,type=type)
x <- as.matrix(x)
nx <- apply(!is.na(x),1,sum)
wx <- matrix(1/nx, nrow=nrow(x),ncol=ncol(x))
wx[is.na(x)] <- 0
x <- as.numeric(x)
w <- as.numeric(wx)
sele <- w > 0
x <- x[sele]
w <- w[sele]
}
list(x=x,w=w,ME=ME)
}
.vasdenest <- function(x,nclass=20,bw){
## 'x' can be a matrix if repeated measures are available
if(any(is.na(x))){
x <- x[!is.na(x)]
warning("missing value(s) removed")
}
stopifnot(is.numeric(nclass))
stopifnot(length(nclass)==1)
nclass <- round(nclass)
stopifnot(nclass>5)
k <- length(x)/3
if(nclass > k) nclass <- k
a <- min(x); b <- max(x)
## determine the number of classes
breaks <- seq(a,b,length=nclass)
xhist <- hist(x, breaks=breaks,plot=FALSE)
y0 <- sqrt(0.25+xhist$counts)
d <- diff(xhist$breaks)
x0 <- xhist$breaks[-1] - d/2
x1 <- (x0 - a)/(b-a)
out <- lpsmooth(y=y0, x=x1, lscv=FALSE,bw=bw)
fhat <- (out$y)^2
x0 <- out$x
d <- diff(x0)
fh <- (fhat[-1]+fhat[-length(fhat)])/2
fTotal <- sum(fh*d)
f0 <- fhat/fTotal/(b-a)
x0 <- x0 * (b-a) + a
list(y=f0, x=x0, pars=out$pars)
}
#.rootunroot <- function(x,breaks,conf.level=0.95){
# name <- deparse(substitute(x))
# if(any(is.na(x)))
# stop("missing value(s) not allowed in 'x'")
# if(any(x<0))
# stop("negative value(s) not allowed in 'x'")
# n <- length(x)
# if(length(breaks)-n != 1)
# stop("invalid length of 'breaks'")
# if(any(diff(breaks)<=0))
# stop("invalid values in 'breaks'")
#}
VAS.ecdf <- function(x,w,alpha=0.05){
if(missing(w)){
if(any(is.na(x)))
x <- x[!is.na(x)]
stopifnot(alpha>0 && alpha<1)
n <- length(x)
xord <- sort(x)
delta <- 0.05*diff(range(xord))
xord <- c(xord[1]-delta,xord)
Fn <- (0:n)/n
en <- sqrt(0.5 * log(2/alpha))/n
Ux <- Fn + en
Ux[Ux>1] <- 1
Lx <- Fn - en
Lx[Lx<0] <- 0
}else{
n <- length(x)
stopifnot(length(w)==n)
stopifnot(all(w>=0))
sele <- is.na(x) | is.na(w)
if(any(sele)){
x <- x[!sele]
w <- w[!sele]
n <- length(x)
}
n0 <- sum(w)
xord <- order(x)
x.ord <- x[xord]
w.ord <- w[xord]
delta <- 0.05*diff(range(x.ord))
Fn <- cumsum(w.ord)/n0
xord <- c(x.ord[1]-delta, x.ord)
Fn <- c(0,Fn)
en <- sqrt(0.5 * log(2/alpha))/n0
Ux <- Fn + en
Ux[Ux>1] <- 1
Lx <- Fn - en
Lx[Lx<0] <- 0
}
structure(list(x=xord,
y=Fn,
lb=Lx,
ub=Ux,
alpha=alpha,
data=x,
ecdf=TRUE),
class="VAS")
}
print.VAS <- function(x,...){
print(summary(x$data))
}
.drawECDF <- function(x,y,ecdf=TRUE,...){
if(ecdf){
k <- length(x)
xn <- 2*x[k] - x[k-1]
x <- sort(c(x[1],x[-1],x[-1],xn))
y <- sort(c(y,y))
}
lines(x,y,...)
}
plot.VAS <- function(x,cb=FALSE,type='n',...){
type <- 'n'
plot(x$y~x$x,type=type,...)
lines(x,cb=cb,...)
}
lines.VAS <- function(x,cb=FALSE,...){
.drawECDF(x$x, x$y,ecdf=x$ecdf,...)
if(cb){
.drawECDF(x$x, x$lb,col='gray',lty=2)
.drawECDF(x$x, x$ub,col='gray',lty=2)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.