Here is a function that I wrote to help me visualize equations of 2 variables (i.e. a 3D surface) that kept cropping up in my Stats 630 class. The core component is fairly simple: I create a 2D grid (a matrix), and for each point on the grid, evaluate the incoming equation or function. The X, Y, and Z coordinates are then shipped off to the RGL package, along with a couple of decorations and RGL does the rest. I’ve also included the empirical marginal distributions for the function on 2 “faces” of the enclosing box. I’ve included a short video of this the code in action at the bottom of this page (below the actual code).
# equationViewer - a native R tool for visualizing equations in 3D
#
# f - a function such as: function(x,y) sin(sqrt(x^2 + y^2))
# or function(x,y) 3*x^1.5 + 2*y
# - in the absence of an actual equation, you can just ship in a data frame or matrix.
# If it has 2 columns, it will think of them as 'x' and 'y' coordinates, each
# representing a single point of data. The code will then convert that to a
# 2D kernel density estimate and plot out the results.
# - Altenatively, if your dataframe or matrix has more than two columns, it will
# think of it as a matrix of 'z' values and plot out their heights.
# xRange - range over which to evaluate the function. Default is 0 to 1.
# yRange - range over which to evaluate the function. Default is 0 to 1.
# gridSize - I find that just 20 or 30 points per axis is plenty. Too many and it will eat
# up your memory and available compute power
# dropFunction - an extra bit. You can supply something like function(x,y) { x < y } and the
# equationViewer will drop out all points on your grid where x is less than y.
# Thus it allows you remove slices from your picture.
# append - add another surface to an already existing surface
#
# Examples:
# equationViewer(f=function(x,y) { sin(x) + sin(y) }, xRange=c(-2,2), yRange=c(-2,2))
# equationViewer(f=data.frame(x=rnorm(1000), y=rt(1000, df=5)))
# equationViewer(f=volcano) # see ?volcano for details
# TO DO:
# - change xRange to xlim, yRange to ylim,
# - args are z, x, y (account for function in different parts (how to do a 3D spiral)
#
equationViewer <- function(f, xRange=c(0,1), yRange=c(0,1), gridSize=21,
dropFunction, append=FALSE, xlab=NULL, ylab=NULL) {
if (!(require(rgl))) {
cat("\n\n Looks like you don't have the package 'rgl' installed\n",
" on this machine. I'll install it now... \n\n")
install.packages("rgl")
}
xRangeTrue <- NULL
yRangeTrue <- NULL
# Cute default function to plot
if (missing(f)) {
f <- function(x,y) sin(sqrt(x^2 + y^2))
xRange = c(-pi,pi); yRange=c(-pi,pi);
}
if (is.function(f)) {
x <- seq(xRange[1], xRange[2], length.out=gridSize) # set up a 20x20 grid of x's and y's
y <- seq(yRange[1], yRange[2], length.out=gridSize)
g <- expand.grid(x=x, y=y)
g$z <- f(g$x, g$y)
mat <- matrix(g$z, nrow=gridSize, ncol=gridSize, byrow=FALSE) # put in z values in a matrix
# example of cutting off the p.d.f. via an "equation"
if (!missing(dropFunction)) mat[dropFunction(g$x, g$y)] <- NA
# Otherwise the user didn't supply a function, rather they supplied a set
# of data consisting of X and Y coordinates that will be slotted into a kde2d
# to give the appearance of a density.
} else if (ncol(f) == 2) {
if (missing(xlab)) xlab = colnames(data.frame(f))[1]
if (missing(ylab)) ylab = colnames(data.frame(f))[2]
dens <- kde2d(f[,1], f[,2], n=gridSize)
x <- dens$x
y <- dens$y
mat <- dens$z
g <- expand.grid(x=x, y=y)
if (!missing(dropFunction)) mat[dropFunction(g$x, g$y)] <- NA
# The following few steps are for graphical purposes
# Artificially inflate the density estimate to make it look nice.
mat <- 0.3 * diff(range(c(x,y))) * mat / max(mat, na.rm=TRUE)
f <- function(x,y) kde2d(x,y)
xRange=range(x); yRange=range(y);
# Last possibility is that they supplied an actual data matrix with the
# cell values as the 'height' values
} else {
mat <- as.matrix(f)
x <- 1:nrow(mat)
y <- 1:ncol(mat)
g <- expand.grid(x=x, y=y)
if (!missing(xRange)) xRangeTrue <- xRange # use for the axes in a minute
if (!missing(yRange)) yRangeTrue <- yRange # use for the axes in a minute
xRange=range(x); yRange=range(y);
f <- function(x,y) Input_Matrix # formalism to keep up with requirements below
if (!missing(dropFunction)) mat[dropFunction(g$x, g$y)] <- NA
}
# as.Table <- function(x) t(x[nrow(x):1,])
mat <- mat[,ncol(mat):1]
heightScale = diff(range(mat, na.rm=TRUE))
colValues = colSums(mat, na.rm=TRUE); # marginal sums of mat
colValues = (colValues - min(colValues)) / diff(range(colValues)) # values from 0..1
if (!all(is.finite(colValues))) colValues <- rep(1,length(colValues))
colValues = colValues * heightScale + min(mat, na.rm=TRUE)
rowValues = rowSums(mat, na.rm=TRUE); # marginal sums of mat
rowValues = (rowValues - min(rowValues)) / diff(range(rowValues)) # values from 0..1
if (!all(is.finite(rowValues))) rowValues <- rep(1,length(rowValues))
rowValues = rowValues * heightScale + min(mat, na.rm=TRUE)
# A couple of steps to assign colors to the heights in the "z matrix"
numColors <- 20
colorLUT <- terrain.colors(numColors+1, alpha=0)
colorIndex <- round(numColors*as.vector(mat)/diff(range(mat, na.rm=TRUE)))
colorIndex <- colorIndex - min(colorIndex, na.rm=TRUE) + 1
colorValues <- colorLUT[colorIndex]
if (!append) open3d(userMatrix=rotationMatrix(pi/3, pi/3, pi/3, 0))
# N.B. rgl thinks of the X and Z axes differently than me.
# It takes a *lot* of horsepower to render partially transparent objects...
if (length(mat) > 25000) {
rgl.surface(x=x, z=y, y=mat, color=colorValues)
# ... but partial alpha transparency is good if you can do it
} else {
rgl.surface(x=x, z=y, y=mat, color=colorValues, alpha=0.7, back="lines", specular="grey")
# surface3d(x=x, y=y, z=mat, color=colorValues, alpha=0.7, back="lines", specular="grey")
}
if (!append) {
# N.B. rgl thinks of the X and Z axes differently than me
xMin <- xRange[1] - 0.05 * diff(xRange)
xSeq <- signif(seq(from=xRange[1], to=xRange[2], length.out=5), 2)
yMin <- yRange[1] - 0.05 * diff(yRange)
ySeq <- signif(seq(from=yRange[1], to=yRange[2], length.out=5), 2)
zRange <- range(mat, na.rm=TRUE)
zMin <- zRange[1] - 0.05 * diff(zRange)
zSeq <- signif(seq(from=zRange[1], to=zRange[2], length.out=5), 2)
# Some gymnastics are required if we want the x- and y- axes to be on a different
# scale of numbers than each other or the z- axis. That is to say, after an
# hour or two of messing around, I can't see a way in RGL to define the x-, y- and
# z- axes to be on different scales without forcing the resulting object into a
# 3D cube (e.g. persp3D).
if (is.null(xRangeTrue)) {
axis3d('x', pos=c(xMin, zMin, yMin), nticks=5, labels=xSeq, at=xSeq, col='red', alpha=0.7)
} else {
xSeqTrue <- signif(seq(from=xRangeTrue[1], to=xRangeTrue[2], length.out=5), 2)
axis3d('x', pos=c(xMin, zMin, yMin), nticks=5, labels=xSeqTrue, at=xSeq, col='red', alpha=0.7)
}
axis3d('y', pos=c(xMin, zMin, yMin), nticks=5, labels=zSeq, at=zSeq, col='red', alpha=0.7)
if (is.null(yRangeTrue)) {
axis3d('z', pos=c(xMin, zMin, yMin), nticks=5, labels=rev(ySeq), at=ySeq, col='red', alpha=0.7)
} else {
ySeqTrue <- signif(seq(from=yRangeTrue[1], to=yRangeTrue[2], length.out=5), 2)
axis3d('z', pos=c(xMin, zMin, yMin), nticks=5, labels=rev(ySeqTrue), at=ySeq, col='red', alpha=0.7)
}
equationName = paste("f(x,y) =", paste(deparse(f)[-1], collapse=""), sep="")
equationName = gsub("[{}]", "", equationName, perl=TRUE)
title3d(main=equationName, family="serif", font=1, color="black", cex=0.8,
xlab=ifelse(is.null(xlab), 'X axis', xlab),
ylab='Z axis',
zlab=ifelse(is.null(ylab), 'Y axis', ylab))
# Add on the marginal density for x (ie f(x))
segments3d(x=rep(x, each=2),
y=unlist(lapply(rowValues, function(a) c(a,zMin))), # interdigitate zeros trick
z=rep(yMin, 2*length(x)), col='red', alpha=0.3, lwd=0.5)
lines3d(x=x, y=rowValues, z=rep(yMin, length(x)), col='red', alpha=0.3, lwd=2)
text3d(x=median(x, na.rm=TRUE), y=median(rowValues, na.rm=TRUE), z=yMin,
alpha=0.8, adj=c(0.5,0.5), family="serif", font=4, cex=0.9, col='red', text="f(x)")
# Add on the marginal density for y (ie f(y))
segments3d(z=rep(y, each=2),
y=unlist(lapply(colValues, function(a) c(a,zMin))), # interdigitate zeros trick
x=rep(xMin, 2*length(x)), col='red', alpha=0.3, lwd=0.5)
lines3d(x=rep(xMin, length(y)), y=colValues, z=y, col='red', alpha=0.3, lwd=2)
text3d(z=median(y, na.rm=TRUE), x=xMin, y=median(colValues, na.rm=TRUE),
alpha=0.8, adj=c(0.5,0.5), family="serif", font=4, cex=0.9, col='red', text="f(y)")
play3d(spin3d(axis=c(0,1,0), rpm=-9), duration=1) # rpm=-5, duration=3
}
rgl.bringtotop()
invisible(list(x=x, y=y, z=mat)) # return the matrix, mat, without printing it
}
Leave a Reply