Nothing
#library(colorspace) is take off so it can go to gamlss
scattersmooth <- function(x, y,
nbin = 100,
lambda = 1,
ndot = 500,
csize = 0.3,
ticks = TRUE,
xlim = c(min(x), max(x)),
ylim = c(min(y), max(y)),
show = TRUE,
save = FALSE,
data = NULL,
xlab = NULL,
ylab = NULL,
cols = heat.colors(10:200),
col.points = "blue",
...)
{
#-------------------------------------------------------------------------------
#-------------------------------------------------------------------------------
smooth2d = function (H, lambda) {
# 2-D exponential smoother
if (length(lambda) == 1) lambda = c(lambda, lambda)
# Build penalty matrices
m <- nrow(H)
n <- ncol(H)
E1 <- diag(m)
E2 <- diag(n)
D1x <- diff(E1)
D1y <- diff(E2)
D2x <- diff(D1x)
D2y <- diff(D1y)
Qx <- E1 + lambda[1]^2 * t(D2x) %*% D2x + 2 * lambda[1] * t(D1x) %*% D1x
Qy <- E2 + lambda[2]^2 * t(D2y) %*% D2y + 2 * lambda[2] * t(D1y) %*% D1y
# Smooth
H <- t(solve(Qy, t(solve(Qx, H))))
return(H)
}
#-------------------------------------------------------------------------------
fillhist <- function(xb, yb, nb)
{
H <- matrix(rep(0, prod(nb)), nb[1], nb[2])
for (i in 1:length(xb))
{
H[xb[i],yb[i]] <- H[xb[i],yb[i]] + 1
}
H
}
#-------------------------------------------------------------------------------
# main function starts here
# check correct input
xlab <- if (is.null(xlab)) deparse(substitute(x)) else xlab
ylab <- if (is.null(ylab)) deparse(substitute(y)) else ylab
y <- if (!is.null(data)) get(deparse(substitute(y)), envir=as.environment(data)) else y
x <- if (!is.null(data)) get(deparse(substitute(x)), envir=as.environment(data)) else x
# ## --------------
# if (!is.null(data))
# {
# YYY <- paste(paste(deparse(substitute(data)),"$", sep=""),deparse(substitute(y)), sep="" )
# y <- eval(parse(text=YYY))
# # XXX <- paste(paste(deparse(substitute(data)),"$", sep=""),deparse(substitute(x)), sep="" )
# x <- eval(parse(text=XXX))
# }
if (length(x) != length(y))
stop("lengths of x and y do not match")
if ( (length(x) < 2) | (length(y) < 2) )
stop("x and y should be vectors")
if ( !is.numeric(x) | !is.numeric(y) )
stop("x and y should contain numeric values")
if ( all(!is.numeric(lambda)) | all(lambda < 0) | (length(lambda) > 2) )
stop("lambda should be numeric and positive")
if ( length(lambda) == 1 )
lambda <- c(lambda, lambda)
if (all(!is.numeric(nbin)) | all(nbin < 1) | (length(nbin) > 2) )
stop("nbin should be a strictly positive integer")
if ( length(nbin) == 1 )
nbin <- c(nbin, nbin)
if (!is.numeric(ndot) | (ndot < 0) | (length(ndot) > 1) )
stop("ndot should be a non-negative integer")
ndot <- floor(ndot)
m <- length(x)
# Put the x-values into bins
xmin <- xlim[1]
xmax <- xlim[2]
dx <- (xmax - xmin) / (nbin[1] - 1)
xbin <- floor(1 + (x - xmin) / dx)
xgrid <- xmin + (1:nbin[1] - 0.5) * dx
# Put the y-values into bins
ymin <- ylim[1]
ymax <- ylim[2]
dy <- (ymax - ymin) / (nbin[2]-1)
ybin <- floor(1 + (y - ymin) / dy)
ygrid <- ymin + (1:nbin[2] - 0.5) * dy
# Select point within boundaries
sel <- 1 <= xbin & xbin <= nbin[1] & 1 <= ybin & ybin <= nbin[2]
xbin <- xbin[sel]
ybin <- ybin[sel]
nmiss <- length(x) - sum(sel)
# Create and smooth the unsmoothed histogram
Hraw <- fillhist(xbin, ybin, nbin)
Hsmooth <- smooth2d(Hraw, lambda)
if (ticks == FALSE) par(xaxt = 'n', yaxt = 'n')
if (show)
{
# cols = #gray(0:n.col/n.col)
#heat_hcl(100:200)
# heat.colors(10:200)
image(x = xgrid, y = ygrid, z = -Hsmooth, xlab = xlab, ylab = ylab,
col = cols)
# Plot selection of dots
if (ndot > 0)
{
ndot <- min(m, ndot)
sel <- sort.list(rnorm(m))[1:ndot]
points(x[sel], y[sel], cex = csize, col = col.points, pch = 15)
}
box()
}
sg <- list()
if (save==TRUE)
{
# Return results
sg <- list(Hraw = Hraw, Hsmooth = Hsmooth, xgrid = xgrid, ygrid = ygrid,
xbin = xbin, ybin = ybin, nmiss = nmiss, seldots = sel)
sg
}
else
{
return(invisible(sg))
}
}
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.