Nothing
#'
#' @rdname dsample
#' @title Generating Random Samples via Wang-Lee algorithm
#' @description \code{dsample} generates a sample of specified size \code{n} from the target density function (up to a normalizing constant) based on the Wang-Lee algorithm.
#' @param expr expression of a target density function
#' @param rpmat matrix containing random points for discretization
#' @param nk positive integer, the number of contours. See \sQuote{Details}.
#' @param n non-negative integer, the desired sample size.
#' @param wconst real number between 0 and 1. See \sQuote{Details}.
#' @details \code{X} has the number of rows equals to the number of discrete base points. In each row, the first element contains the functional value of the target density and the rest elements are the coordinates at which the density is evaluated.
#' \code{wconst} is a constant for adjusting the volume of the last contour.
#' @return \code{dsample} gives the samples in \code{data.frame} with number of rows \code{n} and number of columns \code{ncol(rpmat)}.
#' @references
#' Wang, L. and Lee, C.H. (2014). Discretization-based direct random sample generation. Computational Statistics and Data Analysis, 71, 1001-1010.
#' Lee, C.H. (2009). Efficient Monte Carlo Random Sample Generation through Discretization, MSc thesis, Department of Satistics, University of Manitoba, Canada
#' @keywords sampling discretization
#' @examples
#' ## Example on page 414 in West (1993)
#' expr <- expression((x1*(1-x2))^5 * (x2*(1-x1))^3 * (1-x1*(1-x2)-x2*(1-x1))^37)
#' sets <- list(x1=runif(1e3), x2=runif(1e3))
#' smp <- dsample(expr=expr, rpmat=sets, nk=1e2, n=1e3)
#' @export
dsample <- function(expr, rpmat, n=1e3, nk=1e4, wconst){
variables <- all.vars(expr)
if(missing(wconst)) wconst <- 1
cl <- match.call(expand.dots=FALSE)
stopifnot(is.expression(expr), is.list(rpmat), is.numeric(nk), is.numeric(n), is.numeric(wconst))
y <- eval(expr=expr, envir=rpmat)
X <- as.data.frame(rpmat)
yX <- as.data.frame(cbind(y, X))
yX <- yX[which(y>0),] # data.frame
yX <- yX[order(yX$y, decreasing=TRUE), ]
cnt <- graphics::hist(yX$y, breaks=seq(from=min(yX$y), to=max(yX$y), length.out=nk+1), plot=FALSE)
cnames <- paste("e", seq_len(nk), sep="") # contour name (in order)
yX$cid <- cut(yX$y, cnt$breaks, include.lowest=TRUE)
yX$cnt.name <- cnames[match(yX$cid, names(table(yX$cid)))]
cnt.counts <- cnt$counts
cnt.mids <- cnt$mids
names(cnt.mids) <- names(cnt.counts) <- cnames
gpdf <- rev(cnt.counts * cnt.mids)
gpdf[nk] <- wconst*gpdf[nk]
rev.cntc <- rev(cnt.counts)
cdf <- cumsum(gpdf)/sum(gpdf)
cumcdf <- cumsum(cdf)
pptns <- graphics::hist(stats::runif(n), breaks=c(0,cdf), plot=FALSE)$counts
names(pptns) <- rev(paste("e", seq_len(nk), sep=""))
scnt <- mapply(FUN=sample, MoreArgs=list(replace=TRUE), rev.cntc, pptns)
idx <- unlist( mapply("+", as.list( c(0, cumsum(rev.cntc))[-(nk+1)] ), scnt) )
yX <- yX[order(yX$y, decreasing=TRUE)[idx], ]
robj <- list(expr=expr, yX=yX, X=yX[names(rpmat)], cnt.counts=cnt.counts, cnt.mids=cnt.mids, gpdf=gpdf, cdf=cdf, cumcdf=cumcdf, pptns=pptns, scnt=scnt, idx=idx)
class(robj) <- "dsample"
return(robj)
}
#' @rdname summary.dsample
#' @title Summary Statistics of Marginal Distributions
#' @description Producing basic summary statistics (mean, standard deviation and the first five modes) from the sample drawn for all marginal distributions.
#' @param object \code{data.frame} containing the samples drawn
#' @param n the first n samples
#' @param k number of clusters
#' @param ... arguments passing to the functions used internally
#' @return \code{summary.dsample} gives a list of summary statistics.
#' \item{means}{Means}
#' \item{stdevs}{Standard deviations}
#' \item{modes}{Modes}
#' \item{hc}{object produced by \code{hclust}}
#' \item{grp}{cluster members produced by \code{hclust}}
#' \item{X}{samples generated by \code{dsample}}
#' \item{cdf}{cumulative distributions}
#' @export
summary.dsample <- function(object, n=5, k=1, ...) {
stopifnot(inherits(object, "dsample"))
X <- object$X
hc <- stats::hclust(stats::dist(X))
grp <- stats::cutree(hc, k=k, ...)
cdf <- object$cdf
means <- colMeans(X)
stdevs <- do.call(c, lapply(X, stats::sd, na.rm=TRUE))
modes <- cbind(X)[1:n,]
robj <- list(means=means, stdevs=stdevs, modes=modes, hc=hc, grp=grp, X=X, cdf=cdf)
class(robj) <- "dsample"
return(robj)
}
#' @rdname plot.dsample
#' @title Visualizing Wang-Lee Samples
#' @description The samples generated by the Wang-Lee algorithm are plotted for visual examination. The plot is useful when multiple modes exist.
#' @param x an object produced by \code{dsample}.
#' @param which plot type, 1: CDF, 2: Contours, and 3: Histogram.
#' @param ... arguments passing functions inside
#' @return \code{plot.dsample} has no return value.
#' @export
plot.dsample <- function(x, which, ...){
X <- x$X
cdf <- x$cdf
grp <- x$grp
# graphics::par(mfrow=c(2,2))
if(which==1) graphics::plot(cdf, main="CDF", xlab="E", ylab="F(E)", cex=0.5)
if(which==2) {
graphics::plot(X, cex=0.5, main="Scatter and Contour Plots", xlab="x1", ylab="x2", col=grp)
density <- MASS::kde2d(X[,1], X[,2], n=1e3)
graphics::contour(density, nlevels=5, add=TRUE)
}
if(which==3) graphics::hist(X[,1], main="Histogram", ylab="density", xlab=expression(x[1]), prob=TRUE, breaks=20)
}
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.