| Title: | Dynamic Trees for Learning and Design |
|---|---|
| Description: | Inference by sequential Monte Carlo for dynamic tree regression and classification models with hooks provided for sequential design and optimization, fully online learning with drift, variable selection, and sensitivity analysis of inputs. Illustrative examples from the original dynamic trees paper (Gramacy, Taddy & Polson (2011); <doi:10.1198/jasa.2011.ap09769>) are facilitated by demos in the package; see demo(package="dynaTree"). |
| Authors: | Robert B. Gramacy [aut, cre], Matt A. Taddy [aut], Christoforos Anagnostopoulos [aut] |
| Maintainer: | Robert B. Gramacy <[email protected]> |
| License: | LGPL |
| Version: | 1.2-17 |
| Built: | 2024-11-01 11:26:31 UTC |
| Source: | CRAN |
Inference by sequential Monte Carlo for dynamic tree regression and classification models with hooks provided for sequential design and optimization, fully online learning with drift, variable selection, and sensitivity analysis of inputs. Illustrative examples from the original dynamic trees paper are facilitated by demos in the package; see demo(package="dynaTree")
For a fuller overview including a complete list of functions, and
demos, please use help(package="dynaTree").
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
https://bobby.gramacy.com/r_packages/dynaTree/
plgp, tgp
Uses analytic integration (at the leaves) to calculate the (regression) ALC statistic, or calculates the predictive (class) entropy at the input (X) locations; or calculate ALC at new predictive locations either analytically or numerically
## S3 method for class 'dynaTree' alcX(object, rect = NULL, categ = NULL, approx = FALSE, verb = 0) ## S3 method for class 'dynaTree' entropyX(object, verb = 0) ## S3 method for class 'dynaTree' alc(object, XX, rect = NULL, categ = NULL, approx = FALSE, Xref = NULL, probs = NULL, verb = 0)## S3 method for class 'dynaTree' alcX(object, rect = NULL, categ = NULL, approx = FALSE, verb = 0) ## S3 method for class 'dynaTree' entropyX(object, verb = 0) ## S3 method for class 'dynaTree' alc(object, XX, rect = NULL, categ = NULL, approx = FALSE, Xref = NULL, probs = NULL, verb = 0)
object |
a |
rect |
for |
categ |
A vector of logicals of length |
approx |
a scalar logical that, when |
XX |
a design |
Xref |
|
probs |
weights for the reference locations to be used in a Monte Carlo approximation; usually these weights are class probabilities for response surfaces under constraints |
verb |
a positive scalar integer indicating how many predictive locations
(iterations) after which a progress statement should be
printed to the console; a (default) value of |
This function is most useful for selecting object$X
locations to remove from the analysis, perhaps in an online inference
setting. See retire.dynaTree for more details. The
output is the same as using predict.dynaTree
using XX = object$X, alc = "rect", and Xref =
rect
entropyX only apples to classification models
(object$model != "class"), and alcX applies (only)
to the other, regression, models
The alc function is more generic and allows ALC calculations
at new, predictive, XX locations. This functionality used
to be part of the predict.dynaTree function, but were
separated out for computational reasons. The previous version was
Monte Carlo-based (using Xref) whereas the new version
also allows analytic calculation (now the default, via rect)
The entire object is returned with a new entry called
alcX containing a vector of length nrow(X) with
the ALC values, or entropyX containing the entropy values,
or alc if general ALC calculations at new XX locations
Robert B. Gramacy rbg@vt,
Matt Taddy, and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
Anagnostopoulos, C., Gramacy, R.B. (2013) “Information-Theoretic Data Discarding for Dynamic Trees on Data Streams.” Entropy, 15(12), 5510-5535; arXiv:1201.5568
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, predict.dynaTree, and
retire.dynaTree
## fit the model to the parabola data n <- 100 Xp <- runif(n,-3,3) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) rect <- c(-3,3) out <- dynaTree(Xp, Yp, model="linear", icept="augmented") ## calculate the alcX out <- alcX(out, rect=rect) ## to compare to analytic out <- alc(out, XX=out$X[,-1], rect=rect) ## plot comparison between alcX and predict-ALC plot(out$X[,-1], out$alcX) o <- order(out$X[,2]) lines(out$X[o,-1], out$alc[o], col=2, lty=2) ## now compare to approximate analytic ## (which writes over out$alc) out <- alc(out, XX=out$X[,-1], rect=rect, approx=TRUE) lines(out$X[o,-1], out$alc[o], col=3, lty=3) ## clean up deletecloud(out) ## similarly with entropyX for classification models ## see demo("design") for more iterations and ## design under other active learning heuristics ## like ALC, and EI for optimization; also see ## demo("online") for an online learning example where ## ALC is used for retirement## fit the model to the parabola data n <- 100 Xp <- runif(n,-3,3) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) rect <- c(-3,3) out <- dynaTree(Xp, Yp, model="linear", icept="augmented") ## calculate the alcX out <- alcX(out, rect=rect) ## to compare to analytic out <- alc(out, XX=out$X[,-1], rect=rect) ## plot comparison between alcX and predict-ALC plot(out$X[,-1], out$alcX) o <- order(out$X[,2]) lines(out$X[o,-1], out$alc[o], col=2, lty=2) ## now compare to approximate analytic ## (which writes over out$alc) out <- alc(out, XX=out$X[,-1], rect=rect, approx=TRUE) lines(out$X[o,-1], out$alc[o], col=3, lty=3) ## clean up deletecloud(out) ## similarly with entropyX for classification models ## see demo("design") for more iterations and ## design under other active learning heuristics ## like ALC, and EI for optimization; also see ## demo("online") for an online learning example where ## ALC is used for retirement
A stub for class dynaTree and its custom generic methods
This is just a stub file. See sens.dynaTree
and retire.dynaTree for more information on
the generic methods used in this package
A virtual Class: No objects may be created from it.
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
signature(object = "dynaTree"): ...
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
Gramacy, R.B., Taddy, M.A., and S. Wild (2011). “Variable Selection and Sensitivity Analysis via Dynamic Trees with an Application to Computer Code Performance Tuning” arXiv:1108.4739
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, predict.dynaTree,
update.dynaTree, retire.dynaTree,
sens.dynaTree, alcX.dynaTree
showClass("dynaTree")showClass("dynaTree")
A function to initialize and fit dynamic tree models to regression and classification data by the sequential Monte Carlo (SMC) method of particle learning (PL)
dynaTree(X, y, N = 1000, model = c("constant", "linear", "class", "prior"), nu0s20 = c(0,0), ab = c(0.95, 2), minp = NULL, sb = NULL, nstart = minp, icept = c("implicit", "augmented", "none"), rprop = c("luvar", "luall", "reject"), verb = round(length(y)/10)) dynaTrees(X, y, N = 1000, R = 10, sub = length(y), model = c("constant", "linear", "class", "prior"), nu0s20 = c(0,0), ab=c(0.95, 2), minp = NULL, sb = NULL, nstart = minp, icept = c("implicit", "augmented", "none"), rprop = c("luvar", "luall", "reject"), XX = NULL, yy = NULL, varstats = FALSE, lhs = NULL, plotit = FALSE, proj = 1, rorder = TRUE, verb = round(sub/10), pverb=round(N/10), ...)dynaTree(X, y, N = 1000, model = c("constant", "linear", "class", "prior"), nu0s20 = c(0,0), ab = c(0.95, 2), minp = NULL, sb = NULL, nstart = minp, icept = c("implicit", "augmented", "none"), rprop = c("luvar", "luall", "reject"), verb = round(length(y)/10)) dynaTrees(X, y, N = 1000, R = 10, sub = length(y), model = c("constant", "linear", "class", "prior"), nu0s20 = c(0,0), ab=c(0.95, 2), minp = NULL, sb = NULL, nstart = minp, icept = c("implicit", "augmented", "none"), rprop = c("luvar", "luall", "reject"), XX = NULL, yy = NULL, varstats = FALSE, lhs = NULL, plotit = FALSE, proj = 1, rorder = TRUE, verb = round(sub/10), pverb=round(N/10), ...)
X |
A design |
y |
A vector of length |
N |
a positive scalar integer indicating the number of particles to be used |
R |
a scalar integer |
sub |
Optional argument allowing only a subset of the |
model |
indicates the type of model to be used at the leaves of the tree;
|
nu0s20 |
a two-vector indicating Inverse Gamma prior parameters
|
ab |
tree prior parameter |
minp |
a positive scalar integer describing the smallest allowable
region in the treed partition; if |
sb |
an optional two-vector of positive integers indicating
|
nstart |
a positive scalar integer |
icept |
indicates the type of intertcept term used (only applies to
|
XX |
a design |
yy |
an optional vector of “true” responses at the |
varstats |
if |
lhs |
an optional |
plotit |
a scalar |
proj |
when |
rorder |
a scalar |
rprop |
indicates the scheme used to construct a grow proposal.
The best setting, |
verb |
a positive scalar integer indicating how many time steps
(iterations) should pass before a progress statement is
printed to the console; a value of |
pverb |
a positive scalar integer indicating after many particles
should be processed for prediction before a progress statement is
printed to the console; a value of |
... |
extra arguments to |
The dynaTree function processes the X and y
pairs serially via PL. It builds up a particle cloud
which is stored as an object in C. A “pointer” to that
object is the primary return value. The dynaTrees function
fits several (R) different dynamic tree models on different
time-orderings of the data indices and also
obtains samples from the posterior predictive distribution at
new XX locations. These predictions can be averaged
over each repeat, or used to assess the Monte Carlo predictive
error.
Three different leaf models are supported: two for
regression and one for classification. If model == "class"
then the y values must contain representatives from
every class (1:max(y)). For details of these models and
the complete description of their use at the leaves of the dynamic
trees, see the Taddy, et al., (2009) reference, below.
The tree prior is specified by ab=c(alpha, beta)
via the and minp.
It was originally described by Chipman et al., (1998, 2002)
and subsequently augmented to enforce a minimum number of points
(minp) in each region.
Once a "dynaTree"-class object has been built (by
dynaTree), predictions and estimates of sequential design and
optimization criteria can be obtained via
predict.dynaTree, a generic prediction method.
These values can be used to augment the design, and the
update.dynaTree function can be used to quickly
update the fit with the augmenting data
Both functions return an object of class "dynaTree",
which is a list containing the following fields
m |
|
T |
|
N |
the number of particles used |
X |
a copy of the design matrix |
y |
a copy of the responses |
model |
a copy of the specified leaf model |
params |
a vector containing |
verb |
a copy of the verbosity argument |
lpred |
a vector of |
icept |
a copy of the intercept argument |
time |
the total computing time used to build the particle cloud |
num |
a “pointer” to the |
-
The dynaTrees function can obtain predictive samples
(via predict.dynaTree) at each of the R
repeats. Therefore, the "dynaTree" object returned contains
extra fields collecting these predictive samples, primarily
comprising of R columns of information for each of the fields
returned by predict.dynaTree; see that function for
more details. Likewise, when varstats = TRUE the returned
object also contains vpu, vpt and parde[ fields
whose columns contain the varpropuse and
varproptotal outputs.
Likewise, dynaTrees, can provide variable usage summaries
if varstats = TRUE, in which case the output includes
vpu and vpt fields; See varpropuse
and varproptotal for more details
The dynaTrees function does not return num since
it does not leave any allocated particle clouds on the C-side
As mentioned in the details section, above, the
dynaTree function returns a pointer to a particle
cloud allocated in C. This pointer is used
for prediction, via predict.dynaTree and for
later updating/augmentation of data, via
update.dynaTree.
This information will not be “freed” unless
the user specifically calls deletecloud(num)
or deleteclouds(). Failing to call one
of these functions (when done with the corresponding
object(s)) could result in a memory leak;
see their documentation for more details.
The C-side memory cannot be saved in the workspace,
so they cannot persist across R sessions
To copy a "dynaTree"-class object, use
copy.dynaTree, which will also copy the C-side
memory allocated to the object
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
Gramacy, R.B., Taddy, M.A., and S. Wild (2011). “Variable Selection and Sensitivity Analysis via Dynamic Trees with an Application to Computer Code Performance Tuning” arXiv:1108.4739
Carvalho, C., Johannes, M., Lopes, H., and Polson, N. (2008). “Particle Learning and Smoothing”. Discussion Paper 2008-32, Duke University Dept. of Statistical Science.
Chipman, H., George, E., & McCulloch, R. (1998). Bayesian CART model search (with discussion). Journal of the American Statistical Association, 93, 935–960.
Chipman, H., George, E., & McCulloch, R. (2002). Bayesian treed models. Machine Learning, 48, 303–324.
https://bobby.gramacy.com/r_packages/dynaTree/
predict.dynaTree, update.dynaTree,
plot.dynaTree, deletecloud,
copy.dynaTree, getBF,
varpropuse, varproptotal,
sens.dynaTree, relevance.dynaTree
## simple parabolic data n <- 100 Xp <- sort(runif(n,-3,3)) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) ## fit a piece-wise linear model parab.fit <- dynaTree(Xp, Yp, model="linear") ## obtain predictions at a new set of locations ## and plot parab.fit <- predict(parab.fit, XX=seq(-3, 3, length=100)) plot(parab.fit) ## try duplicating the object parab.fit.copy <- copy(parab.fit) ## must delete the cloud or memory may leak deletecloud(parab.fit); parab.fit$num <- NULL ## to delete all clouds, do: deleteclouds() ## for more examples of dynaTree see update.dynaTree ## Motorcycle accident data if(require("MASS")) { data(mcycle) Xm <- mcycle[,1] Ym <- mcycle[,2] XXm <- seq(min(mcycle[,1]), max(mcycle[,1]), length=100) R <- 2 ## use R >= 10 for better results ## small R is for faster CRAN checks ## fit constant model with R=2 repeats and predictions moto.fit <- dynaTrees(Xm, Ym, XX=XXm, R=R, plotit=TRUE) ## plot the averages plot(moto.fit, ptype="mean") ## clouds automatically deleted by dynaTrees } ## Not run: ## 2-d/3-class classification data library(plgp) library(tgp) xx <- seq(-2, 2, length=20) XX <- expand.grid(xx, xx) X <- dopt.gp(125, Xcand=XX)$XX C <- exp2d.C(X) ## fit a classification model with R=10 repeats, class.fit <- dynaTrees(X, C, XX=XX, model="class") ## for plot the output (no generic plotting available) cols <- c(gray(0.85), gray(0.625), gray(0.4)) par(mfrow=c(1,2)) library(interp) ## plot R-averaged predicted class mclass <- apply(class.fit$p, 1, which.max) image(interp(XX[,1], XX[,2], mclass), col=cols, xlab="x1", ylab="x2", main="repeated class mean") points(X) ## plot R-averaged entropy ment <- apply(class.fit$entropy, 1, mean) image(interp(XX[,1], XX[,2], ment), xlab="x1", ylab="x2", main="repeated entropy mean") ## End(Not run)## simple parabolic data n <- 100 Xp <- sort(runif(n,-3,3)) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) ## fit a piece-wise linear model parab.fit <- dynaTree(Xp, Yp, model="linear") ## obtain predictions at a new set of locations ## and plot parab.fit <- predict(parab.fit, XX=seq(-3, 3, length=100)) plot(parab.fit) ## try duplicating the object parab.fit.copy <- copy(parab.fit) ## must delete the cloud or memory may leak deletecloud(parab.fit); parab.fit$num <- NULL ## to delete all clouds, do: deleteclouds() ## for more examples of dynaTree see update.dynaTree ## Motorcycle accident data if(require("MASS")) { data(mcycle) Xm <- mcycle[,1] Ym <- mcycle[,2] XXm <- seq(min(mcycle[,1]), max(mcycle[,1]), length=100) R <- 2 ## use R >= 10 for better results ## small R is for faster CRAN checks ## fit constant model with R=2 repeats and predictions moto.fit <- dynaTrees(Xm, Ym, XX=XXm, R=R, plotit=TRUE) ## plot the averages plot(moto.fit, ptype="mean") ## clouds automatically deleted by dynaTrees } ## Not run: ## 2-d/3-class classification data library(plgp) library(tgp) xx <- seq(-2, 2, length=20) XX <- expand.grid(xx, xx) X <- dopt.gp(125, Xcand=XX)$XX C <- exp2d.C(X) ## fit a classification model with R=10 repeats, class.fit <- dynaTrees(X, C, XX=XX, model="class") ## for plot the output (no generic plotting available) cols <- c(gray(0.85), gray(0.625), gray(0.4)) par(mfrow=c(1,2)) library(interp) ## plot R-averaged predicted class mclass <- apply(class.fit$p, 1, which.max) image(interp(XX[,1], XX[,2], mclass), col=cols, xlab="x1", ylab="x2", main="repeated class mean") points(X) ## plot R-averaged entropy ment <- apply(class.fit$entropy, 1, mean) image(interp(XX[,1], XX[,2], ment), xlab="x1", ylab="x2", main="repeated entropy mean") ## End(Not run)
Electricity Pricing Data Set Exhibiting Concept Drift
data(elec2)data(elec2)
A data frame with 27552 observations on the following 5 variables.
x1a numeric vector
x2a numeric vector
x3a numeric vector
x4a numeric vector
yclass label
This data has become a benchmark of sorts in streaming classification. It was first described by Harries (1999) and used thereafter for several performance comparisons [e.g., Baena-Garcia et al. (2006); Kuncheva and Plumpton, (2008)]. It holds information for the Australian New South Wales (NSW) Electricity Market, containing 27552 records dated from May 1996 to December 1998, each referring to a period of 30 minutes subsampled as the completely observed portion of 45312 total records with missing values. These records have seven fields: a binary class label, two time stamp indicators (day of week, time), and four covariates capturing aspects of electricity demand and supply.
An appealing property of this dataset is that it is expected to contain drifting data distributions since, during the recording period, the electricity market was expanded to include adjacent areas. This allowed for the production surplus of one region to be sold in the adjacent region, which in turn dampened price levels.
M. Harries. “Splice-2 Comparative Evaluation: Electricity Pricing”. University of New South Wales, School of Computer Science and Engineering technical report (1999)
Anagnostopoulos, C., Gramacy, R.B. (2013) “Information-Theoretic Data Discarding for Dynamic Trees on Data Streams.” Entropy, 15(12), 5510-5535; arXiv:1201.5568
M. Baena-Garcia, J. del Campo-Avila, R., Fidalgo, A. Bifet, R. Gavalda and R. Morales-Bueno. “Early drift detection method”. ECML PKDD 2006 Workshop on Knowledge Discovery from Data Streams, pp. 77-86 (2006)
L.I. Kuncheva C.O. and Plumpton. “Adaptive Learning Rate for Online Linear Discriminant Classifiers”. SSPR and SPR 2008, Lecture Notes in Computer Science (LNCS), 5342, pp. 510-519 (2008)
## this is a snipet from the "elec2" demo; see that demo ## for a full comparison to dynaTree models which can ## cope with drifting concepts ## set up data data(elec2) X <- elec2[,1:4] y <- drop(elec2[,5]) ## predictive likelihood for repated trials T <- 200 ## use nrow(X) for a longer version, ## short T is for faster CRAN checks hits <- rep(NA, T) ## fit the initial model n <- 25; N <- 1000 fit <- dynaTree(X[1:n,], y[1:n], N=N, model="class") w <- 1 for(t in (n+1):T) { ## predict the next data point ## full model fit <- predict(fit, XX=X[t,], yy=y[t]) hits[t] <- which.max(fit$p) == y[t] ## sanity check retiring index if(any(fit$X[w,] != X[t-n,])) stop("bad retiring") ## retire fit <- retire(fit, w) ## update retiring index w <- w + 1; if(w >= n) w <- 1 ## update with new point fit <- update(fit, X[t,], y[t], verb=100) } ## free C-side memory deleteclouds() ## plotting a moving window of hit rates over time rhits <- rep(0, length(hits)) for(i in (n+1):length(hits)) { rhits[i] <- 0.05*as.numeric(hits[i]) + 0.95*rhits[i-1] } ## plot moving window of hit rates plot(rhits, type="l", main="moving window of hit rates", ylab="hit rates", xlab="t")## this is a snipet from the "elec2" demo; see that demo ## for a full comparison to dynaTree models which can ## cope with drifting concepts ## set up data data(elec2) X <- elec2[,1:4] y <- drop(elec2[,5]) ## predictive likelihood for repated trials T <- 200 ## use nrow(X) for a longer version, ## short T is for faster CRAN checks hits <- rep(NA, T) ## fit the initial model n <- 25; N <- 1000 fit <- dynaTree(X[1:n,], y[1:n], N=N, model="class") w <- 1 for(t in (n+1):T) { ## predict the next data point ## full model fit <- predict(fit, XX=X[t,], yy=y[t]) hits[t] <- which.max(fit$p) == y[t] ## sanity check retiring index if(any(fit$X[w,] != X[t-n,])) stop("bad retiring") ## retire fit <- retire(fit, w) ## update retiring index w <- w + 1; if(w >= n) w <- 1 ## update with new point fit <- update(fit, X[t,], y[t], verb=100) } ## free C-side memory deleteclouds() ## plotting a moving window of hit rates over time rhits <- rep(0, length(hits)) for(i in (n+1):length(hits)) { rhits[i] <- 0.05*as.numeric(hits[i]) + 0.95*rhits[i-1] } ## plot moving window of hit rates plot(rhits, type="l", main="moving window of hit rates", ylab="hit rates", xlab="t")
Extract a path (log) Bayes factors (BFs) from the log marginal posterior
probabilities of two "dynaTree"-class objects
getBF(obj1, obj2)getBF(obj1, obj2)
obj1 |
a |
obj2 |
another |
Simply calculates a difference in log marginal posterior
probabilities, setting BFs to zero for initial elements of the
path where one of the objects has more zero marginal probabilities
than the other. The BF is for the model in obj1 over
obj2. If the objects are the output of repeated
fits as obtained from dynaTrees, then multiple
traces are returned
Returns a vector or matrix of a trace(s) of Bayes factors that
can be plotted; see examples below
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
Gramacy, R.B., Taddy, M.A., and S. Wild (2011). “Variable Selection and Sensitivity Analysis via Dynamic Trees with an Application to Computer Code Performance Tuning” arXiv:1108.4739
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, update.dynaTree,
link{logpost}
## parabola data n <- 100 Xp <- sort(runif(n,-3,3)) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) XXp <- seq(-3,3,length=100) ## comparison by log Bayes Factor R <- 2 ## use R >= 10 for better results ## small R is for faster CRAN checks o <- apply(matrix(runif(n*(R-1)), ncol=R-1), 2, order) lpc.p <- dynaTrees(Xp, Yp, R=R, rorder=o, verb=0) lpl.p <- dynaTrees(Xp, Yp, model="linear", R=R, rorder=o, verb=0) bf.p <- getBF(lpl.p, lpc.p) ## plot the log Bayes factors matplot(bf.p, type="l", lty=1, col="gray", main="parabola", xlab="time", ylab="log Bayes factor") ## see demo("reg1d") for further examples## parabola data n <- 100 Xp <- sort(runif(n,-3,3)) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) XXp <- seq(-3,3,length=100) ## comparison by log Bayes Factor R <- 2 ## use R >= 10 for better results ## small R is for faster CRAN checks o <- apply(matrix(runif(n*(R-1)), ncol=R-1), 2, order) lpc.p <- dynaTrees(Xp, Yp, R=R, rorder=o, verb=0) lpl.p <- dynaTrees(Xp, Yp, model="linear", R=R, rorder=o, verb=0) bf.p <- getBF(lpl.p, lpc.p) ## plot the log Bayes factors matplot(bf.p, type="l", lty=1, col="gray", main="parabola", xlab="time", ylab="log Bayes factor") ## see demo("reg1d") for further examples
Plotting predictive distributions constructed from dynamic tree (regression) models for 1-d data – provided primarily for use in our 1-d examples and for illustrative purposes
## S3 method for class 'dynaTree' plot(x, proj = 1, add = FALSE, ylim = NULL, col = 2, lwd = 1, ptype = c("each", "mean"), ...)## S3 method for class 'dynaTree' plot(x, proj = 1, add = FALSE, ylim = NULL, col = 2, lwd = 1, ptype = c("each", "mean"), ...)
x |
a |
add |
a scalar |
proj |
when |
ylim |
user-specified y-axis limits values; see |
col |
user-specified color value; see |
lwd |
user-specified line-width value; see |
ptype |
type of plot used to visualize several predictive
samples obtained from |
... |
other arguments to the generic |
This plotting function only handles the predictive output from
1-dimensional regression dynaTree models as obtained by
first calling dynaTree and then
predict.dynaTree on the resulting output at new
XX locations. It is provided to help make the illustration
of our 1-d examples easier and to serve as an aid in a user's
development of custom plotting functions in higher dimensions
The only output of this function is a pretty plot
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
https://bobby.gramacy.com/r_packages/dynaTree/
predict.dynaTree, dynaTree,
update.dynaTree
## see dynaTree, dynaTrees and update.dynaTree for examples ## which use this plot function## see dynaTree, dynaTrees and update.dynaTree for examples ## which use this plot function
Predicting and calculating sequential design and optimization statistics at new design points (i.e., active learning heuristics) for dynamic tree models
## S3 method for class 'dynaTree' predict(object, XX, yy = NULL, quants = TRUE, ei = FALSE, verb = 0, ...) ## S3 method for class 'dynaTree' coef(object, XX, verb = 0, ...)## S3 method for class 'dynaTree' predict(object, XX, yy = NULL, quants = TRUE, ei = FALSE, verb = 0, ...) ## S3 method for class 'dynaTree' coef(object, XX, verb = 0, ...)
object |
a |
XX |
a design |
yy |
an optional vector of “true” responses at the |
quants |
a scalar |
ei |
a scalar |
verb |
a positive scalar integer indicating how many predictive locations
(iterations) after which a progress statement should be
printed to the console; a (default) value of |
... |
to comply with the generic |
predict returns predictive summary statistics by averaging over the
samples from the posterior predictive distribution obtained
from each of the particles in the cloud pointed to by the
object (object)
coef returns a matrix of regression coefficients used in linear
model leaves (model = "linear") leaves, averaged over all particles,
for each XX location. For other models it prints a warning and
defaults to predict.
The value(s) calculated are appended to object; the new
fields are described below
Note that ALC calculations have been moved to the alc.dynaTree
function(s)
The object returned is of class "dynaTree", which includes a
copy of the list elements from the object passed in,
with the following (predictive)
additions depending on whether object$model is for
regression ("constant" or "linear") or classification
("class").
For regression:
mean |
a vector containing an estimate of the predictive mean
at the |
vmean |
a vector containing an estimate of the variance of predictive mean
at the |
var |
a vector containing an estimate of the predictive
variance (average variance plus variance of mean) at the |
df |
a vector containing the average degrees of freedom at the |
q1 |
a vector containing an estimate of the 5% quantile of
the predictive distribution at the |
q2 |
a vector containing an estimate of the 95% quantile of
the predictive distribution at the |
yypred |
if |
ei |
a vector containing an estimate of the EI statistic,
unless |
;
For classification:
p |
a |
entropy |
a |
;
For coef a new XXc field is created so as not to trample
on XXs that may have been used in a previous predict,
plus
coef |
a |
matrix of particle- averaged regression coefficients.
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, update.dynaTree,
plot.dynaTree, alc.dynaTree,
entropyX.dynaTree
## see the example(s) section(s) of dynaTree and ## update.dynaTree and the demos (demo(package=dynaTree))## see the example(s) section(s) of dynaTree and ## update.dynaTree and the demos (demo(package=dynaTree))
Re-pass the X-y pairs in the object
in a random (or specified) order to temporarily double the
size of the particle set
## S3 method for class 'dynaTree' rejuvenate(object, odr = order(runif(length(object$y))), verb = round(length(object$y)/10))## S3 method for class 'dynaTree' rejuvenate(object, odr = order(runif(length(object$y))), verb = round(length(object$y)/10))
object |
a |
odr |
an integer vector of |
verb |
a positive scalar integer indicating how many time steps
(iterations) should pass before a progress statement is
printed to the console; a value of |
The rejuvenate function causes the particle set to
temporarily double, to have size 2 * object$N. The new
object$N particles represent a discrete approximation
to the dynaTree posterior under the ordering specified
by odr, which may be random. Subsequent calls to
update.dynaTree cause the particle set to revert back
to object$N particles as only that many are obtained from
the particle learning resample step.
This function can be particularly useful in online learning contexts,
where retire.dynaTree is used to retain information
on discarded data, especially when the data is discarded historically
to deal with drifting concepts. Since the new, rejuvenated, particles
are based only on the active data, object$X-object$y
pairs (and not the retired data via informative leaf priors),
subsequent update.dynaTree steps allow the data
to dictate if old (informative prior) or new (default prior) particles
are best for the new concept
The returned list is the same as dynaTree –
i.e., a "dynaTree"-class object but with 2 * object$N
particles. Note that object$N is not updated to reflect this
fact, but the C-side object will indeed have a double particle set.
Repeated calls to rejuvenate will cause the particle set to
double again.
The object (object) must contain a pointer to a particle
cloud (object$num) which has not been deleted by
deletecloud. In particular, it cannot be
an object returned from dynaTrees
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
Anagnostopoulos, C., Gramacy, R.B. (2013) “Information-Theoretic Data Discarding for Dynamic Trees on Data Streams.” Entropy, 15(12), 5510-5535; arXiv:1201.5568
Carvalho, C., Johannes, M., Lopes, H., and Polson, N. (2008). “Particle Learning and Smoothing”. Discussion Paper 2008-32, Duke University Dept. of Statistical Science.
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, alcX.dynaTree,
entropyX.dynaTree, update.dynaTree,
retire.dynaTree
## see retire.dynaTree for a combined example ## illustrating rejuvenation## see retire.dynaTree for a combined example ## illustrating rejuvenation
Computes relevance statistics for each input coordinate by calculating their particle-averaged mean reduction in variance each time that coordinate is used as a splitting variable in (an internal node of) the tree(s)
relevance.dynaTree(object, rect = NULL, categ = NULL, approx = FALSE, verb = 0)relevance.dynaTree(object, rect = NULL, categ = NULL, approx = FALSE, verb = 0)
object |
a |
rect |
an optional |
categ |
A vector of logicals of length |
approx |
a scalar logical indicating if the count of the number of data points in the leaf should be used in place of its area; this can help with numerical accuracy in high dimensional input spaces |
verb |
a positive scalar integer indicating how many particles should
be processed (iterations) before a progress statement should be
printed to the console; a (default) value of |
Each binary split in the tree (in each particle) emits a reduction in variance (for regression models) or a reduction in entropy (for classification). This function calculates these reductions and attributes them to the variable(s) involved in the split(s). Those with the largest relevances are the most useful for prediction. A sensible variable selection rule based on these relevances is to discard those variables whose median relevance is not positive. See the Gramacy, Taddy, & Wild (2011) reference below for more details.
The new set of particles is appended to the old set. However
after a subsequent update.dynaTree call the total
number of particles reverts to the original amount.
Note that this does not work well with dynaTree objects
which were built with model="linear". Rather, a full
sensitivity analysis (sens.dynaTree) is needed. Usually
it is best to first do model="constant" and then use
relevance.dynaTree. Bayes factors (getBF)
can be used to back up any variable selections implied by the
relevance. Then, if desired, one can re-fit on the new (possibly
reduced) set of predictors with model="linear".
There are no caveats with model="class"
The entire object is returned with a new entry called
relevance containing a matrix with ncol(X)
columns. Each row contains the sample from the relevance of
each input, and there is a row for each particle
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Gramacy, R.B., Taddy, M.A., and S. Wild (2011). “Variable Selection and Sensitivity Analysis via Dynamic Trees with an Application to Computer Code Performance Tuning” arXiv:1108.4739
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, sens.dynaTree,
predict.dynaTree
varpropuse, varproptotal
## see the examples in sens.dynaTree for the relevances; ## Also see varpropuse and the class2d demo via ## demo("class2d")## see the examples in sens.dynaTree for the relevances; ## Also see varpropuse and the class2d demo via ## demo("class2d")
Allows the removal (or “retireing”
of X-y pairs from a
"dynaTree"-class object to facilitate online
learning; “retireed” pairs ar absorbed into
the leaf prior(s)
## S3 method for class 'dynaTree' retire(object, indices, lambda = 1, verb = 0)## S3 method for class 'dynaTree' retire(object, indices, lambda = 1, verb = 0)
object |
a |
indices |
a vector of positive integers in |
lambda |
a scalar proportion (forgetting factor) used to downweight the previous prior summary statistics |
verb |
a nonzero scalar causes info about the “retireed” indices,
i.e., their |
Primarily for use in online learning contexts. After
“retireing” the predictive distribution remains unchanged,
because the sufficient statistics of the removed pairs enters
the prior in the leaves of the tree of each particle. Further
update.dynaTree calls (adding data) may cause
changes to the posterior predictive as grow moves cannot keep
the “retires”; see a forthcoming paper for more
details. In many ways, retire.dynaTree is the
opposite of update.dynaTree except that the loss of
information upon “retireing” is not complete.
Drifting regression or classification relationships may be modeled
with a forgetting factor lambda < 1
The alcX.dynaTree provides a good, and computationally
efficient, heuristic for choosing which points to “retire” for
regression models, and likewise link{entropyX.dynaTree} for
classification models.
Note that classification models (model = "class") are
not supported, and implicit intercepts (icept = "implicit")
with linear models (model = "linear") are not supported
at this time
returns a "dynaTree"-class object with updated attributes
In order to use model = "linear" with
dynaTree and retirement one must also specify
icept = "augmented" which automatically augments an
extra column of ones onto the input X design matrix/matrices.
The retire function only supports this icept case
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Anagnostopoulos, C., Gramacy, R.B. (2013) “Information-Theoretic Data Discarding for Dynamic Trees on Data Streams.” Entropy, 15(12), 5510-5535; arXiv:1201.5568
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, alcX.dynaTree,
entropyX.dynaTree, update.dynaTree,
rejuvenate.dynaTree
n <- 100 Xp <- runif(n,-3,3) XX <- seq(-3,3, length=200) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) rect <- c(-3,3) out <- dynaTree(Xp, Yp, model="linear", icept="augmented") ## predict and plot out <- predict(out, XX) plot(out, main="parabola data", lwd=2) ## randomly remove half of the data points out <- retire(out, sample(1:n, n/2, replace=FALSE)) ## predict and add to plot -- shouldn't change anything out <- predict(out, XX) plot(out, add=TRUE, col=3) points(out$X[,-1], out$y, col=3) ## now illustrating rejuvenation, which should result ## in a change to the predictive surface out <- rejuvenate(out) out <- predict(out, XX) plot(out, add=TRUE, col=4) legend("top", c("original", "retired", "rejuvenated"), col=2:4, lty=1) ## clean up deletecloud(out) ## see demo("online") for an online learning example ## where ALC is used for retirementn <- 100 Xp <- runif(n,-3,3) XX <- seq(-3,3, length=200) Yp <- Xp + Xp^2 + rnorm(n, 0, .2) rect <- c(-3,3) out <- dynaTree(Xp, Yp, model="linear", icept="augmented") ## predict and plot out <- predict(out, XX) plot(out, main="parabola data", lwd=2) ## randomly remove half of the data points out <- retire(out, sample(1:n, n/2, replace=FALSE)) ## predict and add to plot -- shouldn't change anything out <- predict(out, XX) plot(out, add=TRUE, col=3) points(out$X[,-1], out$y, col=3) ## now illustrating rejuvenation, which should result ## in a change to the predictive surface out <- rejuvenate(out) out <- predict(out, XX) plot(out, add=TRUE, col=4) legend("top", c("original", "retired", "rejuvenated"), col=2:4, lty=1) ## clean up deletecloud(out) ## see demo("online") for an online learning example ## where ALC is used for retirement
A Monte Carlo sensitivity analysis using random Latin hypercube samples (LHSs) or bootstrap resamples for each particle to estimate main effects as well as 1st order and total sensitivity indices
## S3 method for class 'dynaTree' sens(object, class = NULL, nns = 1000, nME = 100, span = 0.3, method = c("lhs", "boot"), lhs = NULL, categ = NULL, verb = 0)## S3 method for class 'dynaTree' sens(object, class = NULL, nns = 1000, nME = 100, span = 0.3, method = c("lhs", "boot"), lhs = NULL, categ = NULL, verb = 0)
object |
a |
class |
only valid for |
nns |
A positive integer scalar indicating the size of each
LHS or bootstrap drawn for use in the Monte Carlo
integration scheme underlying the sensitivity analysis;
the total number of locations is |
nME |
A positive integer scalar indicating number of grid points, in each input dimension, upon which main effects will be estimated |
span |
A positive real-valued scalar giving the smoothing parameter
for main effects integration: the fraction of |
method |
indicates whether LHS or bootstrap should be used |
lhs |
if
|
categ |
A vector of logicals of length |
verb |
a positive scalar integer indicating how many predictive locations
(iterations) after which a progress statement should be
printed to the console; a (default) value of |
Saltelli (2002) describes a Latin Hypercube sampling based method for estimation of the 'Sobol' sensitivity indices:
1st Order for input ,
where is the -th input.
Total Effect for input ,
where is all inputs except for the -th.
All moments are with respect to the appropriate marginals of the
uncertainty distribution – that is, the probability
distribution on the inputs with respect to which sensitivity is being
investigated.
Under this approach, the integrals involved are approximated through
averages over properly chosen samples based on two LH samples
proportional to U. If nns is the sample size for the
Monte Carlo estimate, this scheme requires nns*(ncol(X)+2)
function evaluations.
The sens.dynaTree function implements the method for unknown functions
, through prediction via one of the tgp regression
models conditional on an observed set of X locations.
For each particle, treated as sample from the dynaTree
model posterior,
the nns*(ncol(X)+2) locations are drawn randomly from the
LHS scheme and realizations of the sensitivity indices are
calculated. Thus we obtain a posterior sample of the indices,
incorporating variability from both the Monte Carlo estimation and
uncertainty about the function output. Since a subset of the
predictive locations are actually an LHS proportional to the
uncertainty distribution, we can also estimate the main effects
through simple non-parametric regression (a moving average).
See the Gramacy, Taddy, & Wild (2011) reference below for more details.
If method = "boot" is used then simply replace LHS above
with a bootstrap resample of the object$X locations.
As with prediction, the dynaTrees function enables
repeated calls to sens.dynaTree
The object returned is of class "dynaTree", which includes a
copy of the list elements from the object passed in,
with the following (sensitivity-analysis specific)
additions.
MEgrid |
An |
MEmean |
A |
MEq1 |
same as |
MEq2 |
same as |
S |
An |
T |
same as |
In the case of object$model = "class" the entries
listed above will themselves be lists with an entry for each
class specified on input, or all classes as is the
default
The quality of sensitivity analysis is dependent on the size of
the LHSs used for integral approximation; as with any Monte
Carlo integration scheme, the sample size (nns) must
increase with the dimensionality of the problem. The total
sensitivity indices are forced non-negative,
and if negative values occur it is necessary to increase
nnd. Postprocessing replaces negative values with NA
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Saltelli, A. (2002) Making best use of model evaluations to compute sensitivity indices. Computer Physics Communications, 145, 280-297.
Gramacy, R.B., Taddy, M.A., and S. Wild (2011). “Variable Selection and Sensitivity Analysis via Dynamic Trees with an Application to Computer Code Performance Tuning” arXiv:1108.4739
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, predict.dynaTree,
relevance.dynaTree,
varpropuse, varproptotal
## friedman data if(require("tgp")) { f <- friedman.1.data(1000) X <- f[,1:6] Z <- f$Y ## fit the model and do the sensitivity analysis N <- 100 ## use N >= 1000 for better results ## small N is for fast CRAN checks out <- dynaTree(X=X, y=Z, N=N, ab=c(0.01,2)) ## also try with model="linear" ## gather relevance statistics out <- relevance(out) boxplot(out$relevance) abline(h=0, col=2, lty=2) ## relevance stats are not as useful when model="linear" ## since it will appear that x4 and x5 not helpful; these ## interact linearly with the response ## full simulation-based sensitivity analysis, the dynaTree:: ## part is only needed if the tgp package is loaded out <- dynaTree::sens(out, verb=100) ## plot the main effects r <- range(rbind(c(out$MEmean, out$MEq1, out$MEq2))) par(mfrow=c(1,ncol(out$X)), mar=c(5,3,2,1)) plot(out$MEgrid[,1], out$MEmean[,1], type="l", ylim=r, lwd=2, ylab="", xlab=colnames(out$MEmean)[1]) lines(out$MEgrid[,1], out$MEq1[,1], lty=2, lwd=2) lines(out$MEgrid[,1], out$MEq2[,1], lty=2, lwd=2) if(ncol(out$X) > 1) { for(d in 2:ncol(out$X)) { plot(out$MEgrid[,d], out$MEmean[,d], col=d, type="l", ylim=r, lwd=2, xlab=colnames(out$MEmean)[d], ylab="") lines(out$MEgrid[,d], out$MEq1[,d], col=d, lty=2) lines(out$MEgrid[,d], out$MEq2[,d], col=d, lty=2) } } ## Sobol indices par(mfrow=c(1,2), mar=c(5,4,4,2)) boxplot(out$S, main="first order indices", xlab="inputs") boxplot(out$T, main="total indices", xlab="inputs") ## these look better when model="linear" ## clean up deletecloud(out) ## for a classification example using the sensitivity hooks ## in the dynaTrees function, see the class2d demo ## i.e., demo("class2d") }## friedman data if(require("tgp")) { f <- friedman.1.data(1000) X <- f[,1:6] Z <- f$Y ## fit the model and do the sensitivity analysis N <- 100 ## use N >= 1000 for better results ## small N is for fast CRAN checks out <- dynaTree(X=X, y=Z, N=N, ab=c(0.01,2)) ## also try with model="linear" ## gather relevance statistics out <- relevance(out) boxplot(out$relevance) abline(h=0, col=2, lty=2) ## relevance stats are not as useful when model="linear" ## since it will appear that x4 and x5 not helpful; these ## interact linearly with the response ## full simulation-based sensitivity analysis, the dynaTree:: ## part is only needed if the tgp package is loaded out <- dynaTree::sens(out, verb=100) ## plot the main effects r <- range(rbind(c(out$MEmean, out$MEq1, out$MEq2))) par(mfrow=c(1,ncol(out$X)), mar=c(5,3,2,1)) plot(out$MEgrid[,1], out$MEmean[,1], type="l", ylim=r, lwd=2, ylab="", xlab=colnames(out$MEmean)[1]) lines(out$MEgrid[,1], out$MEq1[,1], lty=2, lwd=2) lines(out$MEgrid[,1], out$MEq2[,1], lty=2, lwd=2) if(ncol(out$X) > 1) { for(d in 2:ncol(out$X)) { plot(out$MEgrid[,d], out$MEmean[,d], col=d, type="l", ylim=r, lwd=2, xlab=colnames(out$MEmean)[d], ylab="") lines(out$MEgrid[,d], out$MEq1[,d], col=d, lty=2) lines(out$MEgrid[,d], out$MEq2[,d], col=d, lty=2) } } ## Sobol indices par(mfrow=c(1,2), mar=c(5,4,4,2)) boxplot(out$S, main="first order indices", xlab="inputs") boxplot(out$T, main="total indices", xlab="inputs") ## these look better when model="linear" ## clean up deletecloud(out) ## for a classification example using the sensitivity hooks ## in the dynaTrees function, see the class2d demo ## i.e., demo("class2d") }
Updating an already-initialized dynamic tree model with new input/output pairs, primarily to facilitate sequential design and optimization applications
## S3 method for class 'dynaTree' update(object, X, y, verb = round(length(y)/10), ...)## S3 method for class 'dynaTree' update(object, X, y, verb = round(length(y)/10), ...)
object |
a |
X |
an augmenting design |
y |
an augmenting vector of real-valued responses or integer
categories with |
verb |
a positive scalar integer indicating how many time steps
(iterations) after which a progress statement should be
printed to the console; a value of |
... |
to comply with the generic |
This function updates the dynaTree fit with
new (X,y) pairs by the Particle Learning (PL)
algorithm. The updated fit will be for data combined
as rbind(object$X, X) and c(object$y, y).
The primary use of this function is to facilitate sequential
design by optimization and active learning. Typically one
would use predict.dynaTree to estimate active
learning statistics at candidate location.
These are used to pick new (X,y)
locations to add to the design – the new fit being facilitated
by this function; see the examples below
The returned list is the same as dynaTree –
i.e., a "dynaTree"-class object
The object (object) must contain a pointer to a particle
cloud (object$num) which has not been deleted by
deletecloud. In particular, it cannot be
an object returned from dynaTrees
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Taddy, M.A., Gramacy, R.B., and Polson, N. (2011). “Dynamic trees for learning and design” Journal of the American Statistical Association, 106(493), pp. 109-123; arXiv:0912.1586
Anagnostopoulos, C., Gramacy, R.B. (2013) “Information-Theoretic Data Discarding for Dynamic Trees on Data Streams.” Entropy, 15(12), 5510-5535; arXiv:1201.5568
Carvalho, C., Johannes, M., Lopes, H., and Polson, N. (2008). “Particle Learning and Smoothing”. Discussion Paper 2008-32, Duke University Dept. of Statistical Science.
https://bobby.gramacy.com/r_packages/dynaTree/
predict.dynaTree, dynaTree,
plot.dynaTree, deletecloud,
getBF
## simple function describing (x,y) data f1d <- function(x, sd=0.1){ return( sin(x) - dcauchy(x,1.6,0.15) + rnorm(1,0,sd)) } ## initial (x,y) data X <- seq(0, 7, length=30) y <- f1d(X) ## PL fit to initial data obj <- dynaTree(X=X, y=y, N=1000, model="linear") ## a predictive grid XX <- seq(0,7, length=100) obj <- predict(obj, XX, quants=FALSE) ## follow the ALM algorithm and choose the next ## point with the highest predictive variance m <- which.max(obj$var) xstar <- drop(obj$XX[m,]) ystar <- f1d(xstar) ## plot the next chosen point par(mfrow=c(2,1)) plot(obj, ylab="y", xlab="x", main="fitted surface") points(xstar, ystar, col=3, pch=20) plot(obj$XX, sqrt(obj$var), type="l", xlab="x", ylab="predictive sd", main="active learning") ## update the fit with (xstar, ystar) obj <- update(obj, xstar, ystar) ## new predictive surface obj <- predict(obj, XX, quants=FALSE) ## plotted plot(obj, ylab="y", xlab="x", main="updated fitted surface") plot(obj$XX, sqrt(obj$var), type="l", xlab="x", ylab="predictive sd", main="active learning") ## delete the cloud to prevent a memory leak deletecloud(obj); obj$num <- NULL ## see demo("design") for more iterations and ## design under other active learning heuristics ## like ALC, and EI for optimization; also see ## demo("online") for an online learning example## simple function describing (x,y) data f1d <- function(x, sd=0.1){ return( sin(x) - dcauchy(x,1.6,0.15) + rnorm(1,0,sd)) } ## initial (x,y) data X <- seq(0, 7, length=30) y <- f1d(X) ## PL fit to initial data obj <- dynaTree(X=X, y=y, N=1000, model="linear") ## a predictive grid XX <- seq(0,7, length=100) obj <- predict(obj, XX, quants=FALSE) ## follow the ALM algorithm and choose the next ## point with the highest predictive variance m <- which.max(obj$var) xstar <- drop(obj$XX[m,]) ystar <- f1d(xstar) ## plot the next chosen point par(mfrow=c(2,1)) plot(obj, ylab="y", xlab="x", main="fitted surface") points(xstar, ystar, col=3, pch=20) plot(obj$XX, sqrt(obj$var), type="l", xlab="x", ylab="predictive sd", main="active learning") ## update the fit with (xstar, ystar) obj <- update(obj, xstar, ystar) ## new predictive surface obj <- predict(obj, XX, quants=FALSE) ## plotted plot(obj, ylab="y", xlab="x", main="updated fitted surface") plot(obj$XX, sqrt(obj$var), type="l", xlab="x", ylab="predictive sd", main="active learning") ## delete the cloud to prevent a memory leak deletecloud(obj); obj$num <- NULL ## see demo("design") for more iterations and ## design under other active learning heuristics ## like ALC, and EI for optimization; also see ## demo("online") for an online learning example
Calculates the proportion of particles which use each input to make a tree split and the proportion of all splits in trees of each particle that correspond to each input variable; also provides tree height and leaf size summary information
## S3 method for class 'dynaTree' varpropuse(object) ## S3 method for class 'dynaTree' varproptotal(object) ## S3 method for class 'dynaTree' treestats(object)## S3 method for class 'dynaTree' varpropuse(object) ## S3 method for class 'dynaTree' varproptotal(object) ## S3 method for class 'dynaTree' treestats(object)
object |
a |
varpropuse gives the proportion of times a particle
uses each input variable in a tree split; varproptotal gives
the proportion of total uses by the tree in each particle (i.e.,
averaged over the total number of splits used in the tree).
Usually, varpropuse returns a vector of (nearly) all ones
unless there are variables which are not useful in predicting
the response. Using model = "linear" is not recommended
for this sort of variable selection.
treestats returns the average tree height, and the average
leaf size, both active and retired
For varprop*, a
vector of proportions of length ncol(object$X)) is returned;
for treestats a 1-row, 4-column data.frame is
returned
Robert B. Gramacy [email protected],
Matt Taddy and Christoforos Anagnostopoulos
Gramacy, R.B., Taddy, M.A., and S. Wild (2011). “Variable Selection and Sensitivity Analysis via Dynamic Trees with an Application to Computer Code Performance Tuning” arXiv:1108.4739
https://bobby.gramacy.com/r_packages/dynaTree/
dynaTree, sens.dynaTree,
relevance.dynaTree
## ffit a dynaTree model to the Ozone data X <- airquality[,2:4] y <- airquality$Ozone na <- apply(is.na(X), 1, any) | is.na(y) out <- dynaTree(X=X[!na,], y=y[!na]) ## obtain variable usage proportions varpropuse(out) varproptotal(out) ## gather relevance statistics which are more meaningful out <- relevance(out) boxplot(out$relevance) abline(h=0, col=2, lty=2) ## obtain tree statistics treestats(out) ## clean up deletecloud(out)## ffit a dynaTree model to the Ozone data X <- airquality[,2:4] y <- airquality$Ozone na <- apply(is.na(X), 1, any) | is.na(y) out <- dynaTree(X=X[!na,], y=y[!na]) ## obtain variable usage proportions varpropuse(out) varproptotal(out) ## gather relevance statistics which are more meaningful out <- relevance(out) boxplot(out$relevance) abline(h=0, col=2, lty=2) ## obtain tree statistics treestats(out) ## clean up deletecloud(out)