Nothing
R/mcmc_plots.R
In sppmix: Modeling Spatial Poisson and Related Point Processes
Defines functions
Save_AllOpenRglGraphs
Plots_off
plot_convdiags
plot_autocorr
plot_runmean
plot_CompDist
plot.bdmcmc_res
plot.damcmc_res
plot_avgsurf
plot_ind
plot_chains
Documented in
plot_autocorr
plot_avgsurf
plot.bdmcmc_res
plot_chains
plot_CompDist
plot_convdiags
plot.damcmc_res
plot_ind
plot_runmean
Plots_off
Save_AllOpenRglGraphs
#' Plot MCMC chains
#'
#' @description
#' Plot the MCMC chains for all component means and probabilities, generated by \code{\link{est_mix_damcmc}} or
#' \code{\link{est_mix_bdmcmc}}.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_chains}
#'
#' @param fit Object of class \code{damcmc_res} or \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param chain Character vector choosing from c("p", "x", "y"). Multiple
#' choices are supported. This will plot the MCMC chain for the requested variables.
#' @param ncol Number of columns in each plot.
#' @param separate Logical flag to request the the chains should be shown in
#' separate plots or shown in one plot with different colors per component. The
#' latter (\code{separate=FALSE}) is useful for spotting label switching visually.
#' @param open_new_window Open a new window for the plot.
#'
#' @seealso \code{\link{PlotUSAStates}},
#' \code{\link{normmix}},
#' \code{\link{rsppmix}},
#' \code{\link{est_mix_damcmc}},
#' \code{\link{FixLS_da}}
#' @author Jiaxun Chen, Sakis Micheas
#' @examples
#' \donttest{
#' fit <- est_mix_damcmc(pp = spatstat::redwood, m = 10)
#' plot(fit)
#' plot_chains(fit)
#' #plot the chains in the same plot with different colors
#' plot_chains(fit, separate = FALSE)
#' # Only plot the realizations for the component means
#' plot_chains(fit, chain = c("x", "y"))
#' #check labels
#' check_labels(fit)
#' #fix labels and plot the chains again
#' post_fixed = FixLS_da(fit, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed, separate = FALSE)
#' #We work with the California Earthquake data. We fit an IPPP with intensity surface
#' #modeled by a mixture with 8 normal components.
#' CAfit=est_mix_damcmc(CAQuakes2014.RichterOver3.0, m=5, L = 20000)
#' #Now retrieve the surface of Maximum a Posteriori (MAP) estimates of the mixture parameter.
#' #Note that the resulting surface is not affected by label switching.
#' MAPsurf=GetMAPEst(CAfit)
#' #Plot the states and the earthquake locations along with the fitted MAP IPPP intensity surface
#' ret=PlotUSAStates(states=c('California','Nevada','Arizona'), showcentroids=FALSE,
#' shownames=TRUE, main= "Earthquakes in CA, 2014", pp=CAQuakes2014.RichterOver3.0, surf=MAPsurf,
#' boundarycolor="white", namescolor="white")
#' CAfit=est_mix_damcmc(pp = CAQuakes2014.RichterOver3.0, m = 5)
#' plot(CAfit)
#' check_labels(CAfit)
#' plot_chains(CAfit, separate = FALSE)
#' #fix labels and plot the chains again
#' post_fixedCA = FixLS_da(CAfit, plot_result = TRUE)
#' plot_chains(post_fixedCA, separate = FALSE)}
#'
#' @export
plot_chains <- function(fit, burnin = floor(fit$L / 10), chain = c("p", "x", "y"),
ncol = fit$m %% 3 + 1, separate = TRUE,
open_new_window = FALSE) {
iter=ps=xs=ys=Component=NULL
if(class(fit)=="bdmcmc_res")
{
tab=GetBDTable(fit,FALSE)
fit=GetBDCompfit(fit,tab$MAPcomp)$BDgens
burnin = floor(fit$L / 10)
ncol = fit$m %% 3 + 1
}
fit <- drop_realization(fit, burnin)
chain <- match.arg(chain, c("p", "x", "y"), several.ok = TRUE)
plot_df <- data.frame(Component = gl(fit$m, k = fit$L,
labels = paste("Component", 1:fit$m)),
iter = 1:fit$L,
ps = as.vector(fit$genps),
xs = as.vector(t(fit$genmus[, 1, ])),
ys = as.vector(t(fit$genmus[, 2, ])))
# rangeps=range(plot_df$ps)
# rangexs=range(plot_df$xs)
# rangeys=range(plot_df$ys)
if (separate == TRUE) {
plot_p <- ggplot(plot_df, aes(iter, ps)) + geom_path() +
facet_wrap(~ Component, ncol = ncol, scales = "free_y") +
labs(x = "Iteration", y = "Probability") +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeps,
expand=FALSE)+
add_title("Generated mixture probabilities", m = fit$m, L = fit$L)
plot_x <- ggplot(plot_df, aes(iter, xs)) + geom_path() +
facet_wrap(~ Component, ncol = ncol, scales = "free_y") +
labs(x = "Iteration", y = expression(mu[x])) +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangexs,
expand=FALSE)+
add_title("Generated mixture mean, X-coordinate", m = fit$m, L = fit$L)
plot_y <- ggplot(plot_df, aes(iter, ys)) + geom_path() +
facet_wrap(~ Component, ncol = ncol, scales = "free_y") +
labs(x = "Iteration", y = expression(mu[y])) +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeys,
expand=FALSE)+
add_title("Generated mixture mean, Y-coordinate", m = fit$m, L = fit$L)
} else {
plot_p <- ggplot(plot_df, aes(iter, ps, color = Component)) + geom_path() +
labs(x = "Iteration", y = "Probability") +
theme_classic() + theme(panel.border = ggplot2::element_rect(fill = NA,
size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeps,
expand=FALSE)+
add_title("Generated mixture probabilities", m = fit$m, L = fit$L)
plot_x <- ggplot(plot_df, aes(iter, xs, color = Component)) + geom_path() +
labs(x = "Iteration", y = expression(mu[x])) +
theme_classic() + theme(panel.border = ggplot2::element_rect(fill = NA,
size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangexs,
expand=FALSE)+
add_title("Generated mixture mean, X-coordinate", m = fit$m, L = fit$L)
plot_y <- ggplot(plot_df, aes(iter, ys, color = Component)) + geom_path() +
labs(x = "Iteration", y = expression(mu[y])) +
theme_classic() + theme(panel.border = ggplot2::element_rect(fill = NA,
size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeys,
expand=FALSE)+
add_title("Generated mixture mean, Y-coordinate", m = fit$m, L = fit$L)
}
if(open_new_window) {
lapply(paste0("plot_", chain),
function(x) {openwin_sppmix(TRUE); print(eval(parse(text = x)))})
} else {
lapply(paste0("plot_", chain),
function(x) print(eval(parse(text = x))))
}
invisible(NULL)
}
#' Plot membership indicators
#'
#' @description
#' The function plots the posterior means of the
#' membership indicators (or allocation variables)
#' of each point to one of the mixture components,
#' based on a DAMCMC fit. These are the posterior probabilities
#' of a point belonging to a component.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_ind}
#'
#' @param fit Object of class \code{damcmc_res} or \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param open_new_window Open a new window for the plot.
#'
#' @seealso \code{\link{est_mix_damcmc}}
#' @author Sakis Micheas, Yuchen Wang
#' @examples
#' \donttest{
#' fit <- est_mix_damcmc(pp = spatstat::redwood, m = 10)
#' plot_ind(fit)}
#'
#' @export
plot_ind <- function(fit, burnin = floor(fit$L / 10),
open_new_window=FALSE) {
point=component=probability=NULL
if(class(fit)=="bdmcmc_res")
{
tab=GetBDTable(fit,FALSE)
fit=GetBDCompfit(fit,tab$MAPcomp)$BDgens
burnin =floor(fit$L/10)
}
openwin_sppmix(check2open=open_new_window)
fit <- drop_realization(fit, burnin)
labsy=1:fit$m
n=ncol(fit$genzs)
labsx=as.integer(seq(1,n,length=10))
probs <- GetAvgLabelsDiscrete2Multinomial_sppmix(fit$genzs, fit$m)
plot_df <- data.frame(probability = as.vector(probs),
point = 1:n,
component = rep(1:fit$m, each =n))
ggplot(plot_df, aes(point, component, xend = point,
yend = component + 1,
col = probability,
xmin=1,xmax=n,
ymin=1,ymax=fit$m)) +
geom_segment(size = I(5)) +
ggplot2::scale_color_gradient(low = "white", high = "grey18",
guide = guide_colorbar(nbin = 100, barheight = 15)) +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::labs(x = "Point", y = "Component", colour = "Probability") +
ggplot2::scale_x_discrete(limits = labsx)+
ggplot2::scale_y_discrete(limits = labsy+.5,labels=labsy)+
add_title("Membership indicators", m = fit$m, n = n)
}
#' Plot the average intensity surface
#'
#' @description
#' This function calculates the intensity surface
#' at each posterior realization and then computes
#' the average for the intensity surface over a fine grid.
#' The result is a much smoother posterior estimator
#' of the intensity surface, which is not necessarily
#' the same as the surface of posterior means,
#' which is obtained by \code{\link{GetPMEst}}.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_avgsurf}
#'
#' @param fit An object the contains all posterior realizations, e.g., the return
#' value from \code{\link{est_mix_damcmc}} or \code{\link{est_mix_bdmcmc}}.
#' @param win An object of class \code{\link[spatstat]{owin}}.
#' @param LL Length of the side of the square grid.
#' The density or intensity is calculated on an L * L grid.
#' The larger this value is, the slower the calculation,
#' but the better the approximation.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param zlims The limits of the z axis. Defaults to [0,1.1*max(intensity)].
#' @param grayscale Logical flag to request a gray scale plot.
#' @param showplot Logical flag to request that the plot will be shown. Set to FALSE if you want to return the \code{\link[spatstat]{im.object}}, but do not want to produce the 3d plot.
#' @return An image as an object of class \code{\link[spatstat]{im.object}}.
#' @author Jiaxun Chen, Sakis Micheas
#' @seealso \code{\link{rnormmix}},
#' \code{\link{to_int_surf}},
#' \code{\link{rsppmix}},
#' \code{\link{est_mix_damcmc}},
#' \code{\link{plot_density}},
#' \code{\link[ggplot2]{ggtitle}},
#' \code{\link[ggplot2]{geom_point}},
#' \code{\link{plotmix_3d}},
#' \code{\link{GetPMEst}}
#' @examples
#' \donttest{
#' truemix <- rnormmix(m = 5, sig0 = .1, df = 5, xlim= c(-1, 1), ylim =c(0, 3))
#' trueintsurf=to_int_surf(truemix, lambda = 200, win =spatstat::owin( c(-1, 1),c(0, 3)))
#' plot(trueintsurf, main = "True Poisson intensity surface (mixture of normal components)")
#' pp1 <- rsppmix(trueintsurf)
#' # Run Data augmentation MCMC and get posterior realizations
#' postfit=est_mix_damcmc(pp1,m=5)
#' # Plot the average of the surfaces of the posterior realizations
#' avgsurf=plot_avgsurf(postfit, LL = 50)
#' p<-plot_density(as.data.frame(avgsurf))+ggplot2::ggtitle(
#' "Average surface of the posterior realization surfaces\n x denotes a true component mean")
#' #show the point pattern points
#' pp_df <- data.frame(pp1$x,pp1$y)
#' names(pp_df) <- c("x", "y")
#' p<-p + ggplot2::geom_point(data = pp_df,size=0.8)
#' #show the true means
#' mean_df <- data.frame(do.call(rbind, trueintsurf$mus))
#' names(mean_df) <- c("x", "y")
#' p + ggplot2::geom_point(data = mean_df, color = "red", shape = "x", size = 5)
#' #repeat for the contour plot
#' p<-plot_density(as.data.frame(avgsurf),contour = TRUE)+ggplot2::ggtitle(
#' "Average surface of the posterior realization surfaces\n x denotes a true component mean")
#' #show the point pattern points
#' pp_df <- data.frame(pp1$x,pp1$y)
#' names(pp_df) <- c("x", "y")
#' p<-p + ggplot2::geom_point(data = pp_df,size=0.8)
#' #show the true means
#' mean_df <- data.frame(do.call(rbind, trueintsurf$mus))
#' names(mean_df) <- c("x", "y")
#' p + ggplot2::geom_point(data = mean_df, color = "red", shape = "x", size = 5)
#' #plot the 3d surface again based on the returned object
#' plotmix_3d(avgsurf,title1 = paste("Average of", .9*postfit$L,
#' "posterior realizations of the intensity surface"))}
#'
#' @export
plot_avgsurf <- function(fit,
win = fit$data$window, LL = 100,
burnin = floor(fit$L/10), zlims = c(0, 0),
grayscale = FALSE,showplot =TRUE)
{
if(class(fit)=="bdmcmc_res")
{
tab=GetBDTable(fit,FALSE)
fit=GetBDCompfit(fit,tab$MAPcomp)$BDgens
burnin = floor(fit$L/ 10)
}
# cat("\nTo save all open rgl graphs use Save_AllOpenRglGraphs.\n")
# get limits
xlims <- c(win$xrange)
ylims <- c(win$yrange)
L <- fit$L
gridvals <- GetGrid_sppmix(LL,xlims,ylims)
xcoord <- as.vector(gridvals[[1]])
ycoord <- as.vector(gridvals[[2]])
#############ATTENTION, DO NOT TRANSPOSE
zcoord <- ApproxAvgPostIntensity(
fit$allgens_List, fit$genlamdas, LL, burnin,
xlims, ylims, fit$ApproxCompMass)
if(showplot)
{
title1 = paste("Average of",L-burnin,
"posterior realizations of the intensity surface")
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red",
"#7F0000"))
if (grayscale == TRUE) {
col <- gray.colors(100,start = 1, end = 0)[findInterval(zcoord, seq(min(zcoord),
max(zcoord), length = 100))]
} else {
col <- jet.colors(100)[findInterval(zcoord, seq(min(zcoord),
max(zcoord), length = 100))]
}
# rgl::layout3d(matrix(1:2, 1, 2), widths = c(5, 1))
rgl::open3d(windowRect=c(50,50,1000,800),
zoom=1.2)
U=rgl::par3d("userMatrix")
rgl::par3d(userMatrix=
rgl::rotate3d(U,pi/4,0,0,1))
#find the highest z
maxz_height <- max(zcoord)
if (zlims[1] == 0 && zlims[2] == 0) {
zlims <- c(0, 1.1*maxz_height)
}
# Rangez=zlims[2]-zlims[1];
#############ATTENTION, DO NOT TRANSPOSE
rgl::persp3d(x = xcoord, y = ycoord, z = zcoord,
color = col, xlab="x",ylab="y",zlab="",
zlim=c(zlims[1]-0.01,zlims[2]),
box = FALSE, axes = FALSE)
rgl::axis3d('x')
rgl::axis3d('y')
rgl::axis3d('z-+', pos = c(xlims[1], ylims[2], 0))
rgl::rgl.lines(c(xlims[1], xlims[1]),
c(ylims[2], ylims[2]),
c(0,maxz_height),
color = 'black')
rgl::title3d(main=NULL)
rgl::text3d(xlims[2],ylims[2],
zlims[2],texts= title1)
rgl::bgplot3d(suppressWarnings(
fields::image.plot(legend.only = TRUE,
smallplot= c(.8,.82,0.05,.7),
zlim = zlims,
col = jet.colors(100))))
}
return( as.im(list(x=xcoord,
y=ycoord,
z=zcoord)))
}
#' Plot results from a DAMCMC fit
#'
#' @description
#' This function uses the posterior realizations
#' from a \code{\link{est_mix_damcmc}} call, to produce a plethora
#' of plots about the fitted Poisson point process with mixture intensity surface.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot.damcmc_res}
#'
#' @param x Object of class \code{damcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param ... Additional arguments for the S3 method.
#' @author Jiaxun Chen, Sakis Micheas, Yuchen Wang
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{PlotUSAStates}},
#' \code{\link{GetPMEst}},
#' \code{\link{check_labels}},
#' \code{\link{FixLS_da}},
#' \code{\link{plot_chains}}
#' @examples
#' \donttest{
#' fit <- est_mix_damcmc(pp = spatstat::redwood, m = 10)
#' plot(fit)
#' #Tornadoes
#' Tornfit=est_mix_damcmc(Tornadoes2011MO, m=5, L = 20000)
#' MAPsurf=GetMAPEst(Tornfit)
#' ret=PlotUSAStates(states=c('Iowa','Arkansas', 'Missouri','Illinois','Indiana',
#' 'Kentucky','Tennessee', 'Kansas','Nebraska','Texas','Oklahoma','Mississippi',
#' 'Alabama','Louisiana'),showcentroids=FALSE, shownames=TRUE, plotlevels = FALSE,
#' main="Tornadoes about MO, 2011", pp=Tornadoes2011MO, surf=MAPsurf,
#' boundarycolor="white", namescolor="white")
#' print(ret)
#' plot(Tornfit)
#' # get the intensity of posterior means
#' post_mean = GetPMEst(Tornfit)
#' # plot the estimated intensity surface
#' plot(post_mean)
#' #check labels
#' check_labels(Tornfit)
#' # If present then fix label switching, start with approx=TRUE
#' post_fixed = FixLS_da(Tornfit, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed,separate = FALSE)}
#'
#' @export
#' @method plot damcmc_res
plot.damcmc_res <- function(
x,burnin = floor(length(fit$allgens) / 10)
,...) {
fit<-x
post_mix <- GetPMEst(fit, burnin = burnin)
plot(post_mix)
print(plotmix_2d(post_mix, pattern = fit$data))
print(plot_ind(fit, burnin)+add_title("Posterior Means of the membership indicators", m = fit$m, n = fit$data$n))
print(plot_chains(fit, burnin,separate=FALSE))
invisible(NULL)
}
#' Plot results from a BDMCMC fit
#'
#' @description
#' This function uses the posterior realizations
#' from a \code{\link{est_mix_bdmcmc}} call, to produce a plethora
#' of plots about the fitted Poisson point process with mixture intensity surface.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot.bdmcmc_res}
#'
#' @param x Object of class \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param win An object of class \code{\link[spatstat]{owin}}.
#' @param LL Length of the side of the square grid.
#' The density or intensity is calculated on an L * L grid.
#' The larger this value is, the slower the calculation,
#' but the better the approximation as well as the smoother the resulting plots.
#' @param zlims The limits of the z axis. Defaults to [0,max(z)].
#' @param ... Additional arguments for the S3 method.
#' @details Unlike the corresponding output from DAMCMC (fixed number of components), the BDMCMC algorithm
#' allows us to obtain a distribution for the number of components
#' which can be thought of as a statistical inference procedure for model
#' selection. In particular, the Bayesian model average of all the realizations
#' is perhaps the best possible estimator of the Poisson intensity surface, however,
#' it can be very slow to compute for moderate number of iterations and maximum
#' number of components allowed.
#' @return A list with the following objects
#' \item{FreqTab}{Frequency table for the number of components.}
#' \item{Mapsurf}{The Maximum A Posteriori (MAP) Poisson intensity surface based on the corresponding posterior means (label switching might be present). The MAP is the mode of the distribution of the number of components.}
#' \item{BMA}{The Bayesian Model Average is returned only if we answer "Y" to request it. Alternatively, use function \code{\link{GetBMA}} to compute the BMA. This is an image, i.e., an object of class \code{\link[spatstat]{im.object}}.}
#' @author Jiaxun Chen, Sakis Micheas, Yuchen Wang
#' @seealso \code{\link{est_mix_bdmcmc}},
#' \code{\link{PlotUSAStates}},
#' \code{\link{plotmix_3d}},
#' \code{\link{plot_density}},
#' \code{\link{check_labels}},
#' \code{\link{FixLS_da}},
#' \code{\link{plot_chains}},
#' \code{\link{GetBMA}}
#' \code{\link{GetBDTable}},
#' @examples
#' \donttest{
#' fit <- est_mix_bdmcmc(pp = spatstat::redwood, m = 10)
#' plot(fit)
#' #Tornadoes
#' ret=PlotUSAStates(states=c('Iowa','Arkansas', 'Missouri','Illinois','Indiana','Kentucky',
#' 'Tennessee', 'Kansas','Nebraska','Texas','Oklahoma','Mississippi', 'Alabama','Louisiana'),
#' showcentroids=FALSE,shownames=TRUE, plotlevels = FALSE, main="Tornadoes about MO, 2011",
#' pp=Tornadoes2011MO)
#' print(ret)
#' Tornfit=est_mix_bdmcmc(pp = Tornadoes2011MO, m = 7)
#' TornResults=plot(Tornfit)#if we plot the Bayesian model average return it
#' TornResults
#' if(!is.null(TornResults$BMA)){
#' BMA_image=TornResults$BMA#must answer yes above or compute it using GetBMA
#' burnin=.1*Tornfit$L
#' title1 = paste("Bayesian model average of", Tornfit$L-burnin,"posterior realizations")
#' plotmix_3d(BMA_image,title1=title1)
#' plot_density(as.data.frame(BMA_image))+ggplot2::ggtitle(
#' "Bayesian Model Average Intensity")
#' plot_density(as.data.frame(BMA_image),TRUE)+ggplot2::ggtitle(
#' "Contours of the Bayesian Model Average Intensity")}
#' # Work with the MAP intensity
#' Mapsurf=TornResults$Mapsurf
#' plot(Mapsurf)
#' #retrieve realizations for the MAP number of components only
#' tab=GetBDTable(Tornfit)
#' MAPm=tab$MAPcomp
#' BDfitMAPcomp=GetBDCompfit(Tornfit,MAPm)
#' summary(BDfitMAPcomp)
#' summary(BDfitMAPcomp$BDgens)
#' plot(BDfitMAPcomp$BDsurf,main=paste(
#' "Poisson Mixture intensity surface, MAP # of components=",MAPm))
#' #check labels
#' check_labels(BDfitMAPcomp$BDgens)
#' # If present then fix label switching, start with approx=TRUE
#' post_fixed = FixLS_da(BDfitMAPcomp$BDgens, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed,separate = FALSE)
#' #this one works better
#' post_fixed = FixLS_da(BDfitMAPcomp$BDgens,approx=FALSE, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed,separate = FALSE)}
#'
#' @export
#' @method plot bdmcmc_res
plot.bdmcmc_res <- function(x, win = fit$data$window,
burnin = floor(fit$L/10),
LL = 100, zlims = c(0, 0),...) {
fit<-x
L <- fit$L
xlims <- win$xrange
ylims <- win$yrange
BDtab=GetBDTable(fit,FALSE)
tab=BDtab$FreqTab
MAPcomp=BDtab$MAPcomp
plot_CompDist(fit)
mapcompgens=GetBDCompfit(fit,MAPcomp)
mapgens=mapcompgens$BDgens
BDnormmix=mapcompgens$BDnormmix
mean_lambda=mean(mapgens$genlamdas);
PP=fit$data
n=ncol(fit$genzs)
plot2dPP(PP,BDnormmix$mus,title1="Point pattern along with the mixture means for MAP m")
intsurfMAP=mapcompgens$BDsurf
print(summary(intsurfMAP))
print(plotmix_2d(intsurfMAP, PP)+add_title("MAP Poisson intensity surface along with the point pattern",lambda =intsurfMAP$lambda,m=intsurfMAP$m,n=PP$n,L=mapcompgens$BDgens$L))
print(plotmix_2d(intsurfMAP, PP,contour = TRUE)+add_title("MAP Poisson intensity surface contours along with the point pattern",lambda =intsurfMAP$lambda,m=intsurfMAP$m,n=PP$n,L=mapcompgens$BDgens$L))
print(plot_chains(fit, burnin,separate=FALSE))
dens_image=NULL
if(Get_User_Input_sppmix("Compute Bayesian model average (slow operation)?"))
{
dens_image=GetBMA(fit)
title1 = paste("Bayesian model average of",L-burnin,"posterior realizations")
plotmix_3d(dens_image,title1=title1)
}
if(is.null(dens_image))
return(list(FreqTab=tab,Mapsurf=intsurfMAP))
else
return(list(FreqTab=tab,Mapsurf=intsurfMAP,BMA=dens_image))
}
#' Plots for the number of components
#'
#' @description
#' The function produces two plots: the trace plot for the number of components based on
#' all realizations of a BDMCMC fit, and a barplot describing the distribution
#' of the number of components.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_CompDist}
#'
#' @param fit Object of class \code{bdmcmc_res}.
#' @param open_new_window Open new windows for the two plots.
#' @seealso \code{\link{est_mix_bdmcmc}}
#' @author Sakis Micheas
#' @examples
#' \donttest{
#' fitBD <- est_mix_bdmcmc(spatstat::redwood, m = 10)
#' plot_CompDist(fitBD)
#' CAfitBD=est_mix_bdmcmc(pp = CAQuakes2014.RichterOver3.0, m = 10)
#' plot_CompDist(CAfitBD)}
#'
#' @export
plot_CompDist=function(fit,open_new_window=FALSE)
{
openwin_sppmix(open_new_window)
BDtab=GetBDTable(fit)
tab=BDtab$FreqTab
oldpar=par(mai=c(1.2,1.2,1,1))
mp=barplot(tab,names.arg=1:fit$maxnumcomp,
xlab="Number of Components",ylab=paste("Iterations,",fit$L,"total"),
main="Distribution of the number of components"
,ylim=c(0,1.2*max(tab)))
for(i in 1:fit$maxnumcomp)
text(x=mp,y=tab+0.1*max(tab),labels=tab)
openwin_sppmix(open_new_window)
plot(fit$numcomp,ylim=c(1,fit$maxnumcomp),
xlab="Iteration",
ylab="Number of components",
type="l",main="Generated chain for the number of components")
suppressWarnings( par(oldpar))
}
#' Checking convergence: running means plot
#'
#' @description
#' This function produces a running mean plot, that is,
#' a plot of the iterations against the mean of
#' the draws up to each iteration. If the plot is not a near
#' constant line then convergence has not been achieved (e.g., label switching is present).
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_runmean}
#'
#' @param chain An \code{Lx1} vector containing the \code{L} posterior realizations.
#' @param open_new_window Open a new window for the plot.
#'
#' @author Sakis Micheas
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{rmixsurf}},
#' \code{\link{rsppmix}}
#' @examples
#' \donttest{
#' truemix_surf <- rmixsurf(m = 3, lambda=100, xlim = c(-3,3), ylim = c(-3,3))
#' plot(truemix_surf)
#' genPPP=rsppmix(intsurf = truemix_surf, truncate = FALSE)
#' fit <- est_mix_damcmc(pp = genPPP, m = 3)
#' plot_runmean(fit$genps[,1])
#' plot_runmean(fit$genps[,2])
#' plot_runmean(fit$genps[,3])}
#'
#' @export
plot_runmean<- function(chain,open_new_window = FALSE)
{
Iteration=PostMean=NULL
openwin_sppmix(check2open=open_new_window)
L=length(chain)
runmean=rep(0,L)
for(i in 1:L)
{
runmean[i]=mean(chain[1:i])
}
# plot(1:L,runmean,type="l")
plot_df <- data.frame(Iteration = 1:L,
PostMean =runmean)
plot_p <- ggplot(plot_df,
aes(Iteration, PostMean)) + geom_path() +
labs(x = "Iteration", y = "Posterior Mean") +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:L,expand=FALSE)+
ggplot2::theme(plot.margin=ggplot2::unit(c(1,1,1,0), "lines"))+
add_title("Running mean plot",L = L)
plot_p
}
#' Checking convergence: autocorrelation plot
#'
#' @description
#' This function can be used to assess
#' convergence by visualizing the
#' autocorrelations between the draws of the
#' Markov chain \code{chain}. The \code{lag k}
#' autocorrelation \code{rho_k} is the correlation between every draw
#' and its \code{kth} lag. We would expect the
#' \code{kth} lag autocorrelation to be
#' smaller as \code{k} increases (that is, the
#' 100th and 1000th draws should be less
#' correlated than the 100th and 105th draws).
#' For higher values of k we anticipate small
#' autocorrelation values, otherwise the chain
#' is not \code{mixing} well (in other words we do not
#' explore the parameter space adequately).
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_autocorr}
#'
#' @param chain An \code{Lx1} vector containing the \code{L} posterior realizations.
#' @param open_new_window Open a new window for the plot.
#' @param maxlag The maximum lag value to consider. Default is 100.
#' @author Sakis Micheas
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{rmixsurf}},
#' \code{\link{rsppmix}}
#' @examples
#' \donttest{
#' truemix_surf <- rmixsurf(m = 3, lambda=100, xlim = c(-3,3), ylim = c(-3,3))
#' plot(truemix_surf)
#' genPPP=rsppmix(intsurf = truemix_surf, truncate = FALSE)
#' fit <- est_mix_damcmc(pp = genPPP, m = 3)
#' plot_autocorr(fit$genps[,1])
#' plot_autocorr(fit$genps[,2])
#' plot_autocorr(fit$genps[,3])}
#'
#' @export
plot_autocorr<- function(chain,
open_new_window = FALSE,maxlag=100)
{
Autocorr=NULL
openwin_sppmix(check2open=open_new_window)
L=length(chain)
if(L<=maxlag)
stop("Reduce lag or increase number of iterations.")
autocorr=rep(0,maxlag)
xbar=mean(chain)
xvar=(L-1)*var(chain)
for(k in 1:maxlag)
{
autocorr[k]=sum((chain[1:(L-k)]-xbar)*
(chain[(k+1):L]-xbar))/xvar
}
plot_df <- data.frame(lag = 1:L,
Autocorr =autocorr)
plot_p <- ggplot(plot_df,
aes(lag, Autocorr)) + geom_path() +
labs(x = "Lag", y = "Autocorrelation") +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:maxlag,expand=FALSE)+
ggplot2::theme(plot.margin=ggplot2::unit(c(1,1,1,0), "lines"))+
add_title("Autocorrelation plot")
plot_p
}
#' Checking convergence visually
#'
#' @description
#' Based on a `damcmc_res` object, this function will produce
#' many graphs to help assess convergence visually,
#' including running mean plots and
#' autocorrelation plots for all the parameters. This function calls
#' \code{\link{plot_runmean}} and \code{\link{plot_autocorr}} for all parameters
#' so we do not have to it individually.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_convdiags}
#'
#' @param fit Object of class \code{damcmc_res} or \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param open_new_window Open a new window for the plot.
#' @param maxlag The maximum lag value to consider. Default is 100.
#' @author Sakis Micheas
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{rmixsurf}},
#' \code{\link{plot_runmean}},
#' \code{\link{plot_autocorr}},
#' \code{\link{rsppmix}}
#' @examples
#' \donttest{
#' truemix_surf <- rmixsurf(m = 3, lambda=100, xlim = c(-3,3), ylim = c(-3,3))
#' plot(truemix_surf)
#' genPPP=rsppmix(intsurf = truemix_surf, truncate = FALSE)
#' fit = est_mix_damcmc(pp = genPPP, m = 3)
#' plot_convdiags(fit)}
#'
#' @export
plot_convdiags<- function(fit,burnin = floor(fit$L / 10),
open_new_window = FALSE,maxlag=100)
{
fit <- drop_realization(fit, burnin)
open_new_window1=open_new_window
maxlag1=maxlag
if(fit$L<=maxlag1)
stop("Reduce lag or increase number of iterations.")
#probs first
for(j in 1:fit$m)
{
print(plot_runmean(fit$genps[,j],open_new_window=open_new_window1)+ggtitle(paste("Component probability",j)))
print(plot_autocorr(fit$genps[,j],
open_new_window=open_new_window1, maxlag =maxlag1)+ggtitle(paste("Component probability",j)))
}
#mu x-coords
for(j in 1:fit$m)
{
print(plot_runmean(fit$genmus[j,1,],open_new_window=open_new_window1)+ggtitle(paste("x-coord for mu, component",j)))
print(plot_autocorr(fit$genmus[j,1,],
open_new_window=open_new_window1, maxlag =maxlag1)+ggtitle(paste("x-coord for mu, component",j)))
}
#mu y-coords
for(j in 1:fit$m)
{
print(plot_runmean(fit$genmus[j,2,],open_new_window=open_new_window1)+ggtitle(paste("y-coord for mu, component",j)))
print(plot_autocorr(fit$genmus[j,2,],
open_new_window=open_new_window1, maxlag =maxlag1)+ggtitle(paste("y-coord for mu, component",j)))
}
}
#' Closes all open plots
#'
#' @description
#' The function closes all Rgl plots, as well as, any graphics devices.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #Plots_off}
#'
#' @seealso \code{\link{Save_AllOpenRglGraphs}}
#'
#' @author Sakis Micheas
#' @examples
#' \donttest{
#' Plots_off()}
#'
#' @export
Plots_off<- function()
{
graphics.off()
while (rgl::rgl.cur()>0)
{
rgl::rgl.close()
}
}
#' Saves RGL plots
#'
#' @description
#' The function saves all open RGL plots (3d plots) to the specified directory
#' and using a template name.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #Save_AllOpenRglGraphs}
#'
#' @param dir1 Directory to save the plots.
#' @param filename1 Template filename. Each open plot is saved as filename1.png, filename2.png, and so forth.
#' @author Sakis Micheas
#' @seealso \code{\link{Plots_off}}
#' @examples
#' \donttest{
#' # use a temporary directory to save the plots
#' Save_AllOpenRglGraphs(dir1=tempdir())}
#'
#' @export
Save_AllOpenRglGraphs<- function(
dir1,filename1="RglGraph")
{
if(missing(dir1))
stop("Please specify a directory wherein to save the RGL plots.")
# graphics.off()
count=1;
while (rgl::rgl.cur()>0)
{
rgl::rgl.snapshot(paste(filename1,count,".png",sep = ""))
rgl::rgl.close()
count=count+1;
}
if (count==2)
cat(paste("Saved 1 plot in",dir1))
else
cat(paste("Saved",count-1,"plots in",dir1))
}
Try the sppmix package in your browser
Any scripts or data that you put into this service are public.
sppmix documentation built on Jan. 13, 2021, 10:04 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Save_AllOpenRglGraphs Plots_off plot_convdiags plot_autocorr plot_runmean plot_CompDist plot.bdmcmc_res plot.damcmc_res plot_avgsurf plot_ind plot_chains
Documented in plot_autocorr plot_avgsurf plot.bdmcmc_res plot_chains plot_CompDist plot_convdiags plot.damcmc_res plot_ind plot_runmean Plots_off Save_AllOpenRglGraphs
#' Plot MCMC chains
#'
#' @description
#' Plot the MCMC chains for all component means and probabilities, generated by \code{\link{est_mix_damcmc}} or
#' \code{\link{est_mix_bdmcmc}}.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_chains}
#'
#' @param fit Object of class \code{damcmc_res} or \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param chain Character vector choosing from c("p", "x", "y"). Multiple
#' choices are supported. This will plot the MCMC chain for the requested variables.
#' @param ncol Number of columns in each plot.
#' @param separate Logical flag to request the the chains should be shown in
#' separate plots or shown in one plot with different colors per component. The
#' latter (\code{separate=FALSE}) is useful for spotting label switching visually.
#' @param open_new_window Open a new window for the plot.
#'
#' @seealso \code{\link{PlotUSAStates}},
#' \code{\link{normmix}},
#' \code{\link{rsppmix}},
#' \code{\link{est_mix_damcmc}},
#' \code{\link{FixLS_da}}
#' @author Jiaxun Chen, Sakis Micheas
#' @examples
#' \donttest{
#' fit <- est_mix_damcmc(pp = spatstat::redwood, m = 10)
#' plot(fit)
#' plot_chains(fit)
#' #plot the chains in the same plot with different colors
#' plot_chains(fit, separate = FALSE)
#' # Only plot the realizations for the component means
#' plot_chains(fit, chain = c("x", "y"))
#' #check labels
#' check_labels(fit)
#' #fix labels and plot the chains again
#' post_fixed = FixLS_da(fit, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed, separate = FALSE)
#' #We work with the California Earthquake data. We fit an IPPP with intensity surface
#' #modeled by a mixture with 8 normal components.
#' CAfit=est_mix_damcmc(CAQuakes2014.RichterOver3.0, m=5, L = 20000)
#' #Now retrieve the surface of Maximum a Posteriori (MAP) estimates of the mixture parameter.
#' #Note that the resulting surface is not affected by label switching.
#' MAPsurf=GetMAPEst(CAfit)
#' #Plot the states and the earthquake locations along with the fitted MAP IPPP intensity surface
#' ret=PlotUSAStates(states=c('California','Nevada','Arizona'), showcentroids=FALSE,
#' shownames=TRUE, main= "Earthquakes in CA, 2014", pp=CAQuakes2014.RichterOver3.0, surf=MAPsurf,
#' boundarycolor="white", namescolor="white")
#' CAfit=est_mix_damcmc(pp = CAQuakes2014.RichterOver3.0, m = 5)
#' plot(CAfit)
#' check_labels(CAfit)
#' plot_chains(CAfit, separate = FALSE)
#' #fix labels and plot the chains again
#' post_fixedCA = FixLS_da(CAfit, plot_result = TRUE)
#' plot_chains(post_fixedCA, separate = FALSE)}
#'
#' @export
plot_chains <- function(fit, burnin = floor(fit$L / 10), chain = c("p", "x", "y"),
ncol = fit$m %% 3 + 1, separate = TRUE,
open_new_window = FALSE) {
iter=ps=xs=ys=Component=NULL
if(class(fit)=="bdmcmc_res")
{
tab=GetBDTable(fit,FALSE)
fit=GetBDCompfit(fit,tab$MAPcomp)$BDgens
burnin = floor(fit$L / 10)
ncol = fit$m %% 3 + 1
}
fit <- drop_realization(fit, burnin)
chain <- match.arg(chain, c("p", "x", "y"), several.ok = TRUE)
plot_df <- data.frame(Component = gl(fit$m, k = fit$L,
labels = paste("Component", 1:fit$m)),
iter = 1:fit$L,
ps = as.vector(fit$genps),
xs = as.vector(t(fit$genmus[, 1, ])),
ys = as.vector(t(fit$genmus[, 2, ])))
# rangeps=range(plot_df$ps)
# rangexs=range(plot_df$xs)
# rangeys=range(plot_df$ys)
if (separate == TRUE) {
plot_p <- ggplot(plot_df, aes(iter, ps)) + geom_path() +
facet_wrap(~ Component, ncol = ncol, scales = "free_y") +
labs(x = "Iteration", y = "Probability") +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeps,
expand=FALSE)+
add_title("Generated mixture probabilities", m = fit$m, L = fit$L)
plot_x <- ggplot(plot_df, aes(iter, xs)) + geom_path() +
facet_wrap(~ Component, ncol = ncol, scales = "free_y") +
labs(x = "Iteration", y = expression(mu[x])) +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangexs,
expand=FALSE)+
add_title("Generated mixture mean, X-coordinate", m = fit$m, L = fit$L)
plot_y <- ggplot(plot_df, aes(iter, ys)) + geom_path() +
facet_wrap(~ Component, ncol = ncol, scales = "free_y") +
labs(x = "Iteration", y = expression(mu[y])) +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeys,
expand=FALSE)+
add_title("Generated mixture mean, Y-coordinate", m = fit$m, L = fit$L)
} else {
plot_p <- ggplot(plot_df, aes(iter, ps, color = Component)) + geom_path() +
labs(x = "Iteration", y = "Probability") +
theme_classic() + theme(panel.border = ggplot2::element_rect(fill = NA,
size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeps,
expand=FALSE)+
add_title("Generated mixture probabilities", m = fit$m, L = fit$L)
plot_x <- ggplot(plot_df, aes(iter, xs, color = Component)) + geom_path() +
labs(x = "Iteration", y = expression(mu[x])) +
theme_classic() + theme(panel.border = ggplot2::element_rect(fill = NA,
size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangexs,
expand=FALSE)+
add_title("Generated mixture mean, X-coordinate", m = fit$m, L = fit$L)
plot_y <- ggplot(plot_df, aes(iter, ys, color = Component)) + geom_path() +
labs(x = "Iteration", y = expression(mu[y])) +
theme_classic() + theme(panel.border = ggplot2::element_rect(fill = NA,
size = 1)) +
ggplot2::coord_cartesian(xlim =1:fit$L,
#ylim =rangeys,
expand=FALSE)+
add_title("Generated mixture mean, Y-coordinate", m = fit$m, L = fit$L)
}
if(open_new_window) {
lapply(paste0("plot_", chain),
function(x) {openwin_sppmix(TRUE); print(eval(parse(text = x)))})
} else {
lapply(paste0("plot_", chain),
function(x) print(eval(parse(text = x))))
}
invisible(NULL)
}
#' Plot membership indicators
#'
#' @description
#' The function plots the posterior means of the
#' membership indicators (or allocation variables)
#' of each point to one of the mixture components,
#' based on a DAMCMC fit. These are the posterior probabilities
#' of a point belonging to a component.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_ind}
#'
#' @param fit Object of class \code{damcmc_res} or \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param open_new_window Open a new window for the plot.
#'
#' @seealso \code{\link{est_mix_damcmc}}
#' @author Sakis Micheas, Yuchen Wang
#' @examples
#' \donttest{
#' fit <- est_mix_damcmc(pp = spatstat::redwood, m = 10)
#' plot_ind(fit)}
#'
#' @export
plot_ind <- function(fit, burnin = floor(fit$L / 10),
open_new_window=FALSE) {
point=component=probability=NULL
if(class(fit)=="bdmcmc_res")
{
tab=GetBDTable(fit,FALSE)
fit=GetBDCompfit(fit,tab$MAPcomp)$BDgens
burnin =floor(fit$L/10)
}
openwin_sppmix(check2open=open_new_window)
fit <- drop_realization(fit, burnin)
labsy=1:fit$m
n=ncol(fit$genzs)
labsx=as.integer(seq(1,n,length=10))
probs <- GetAvgLabelsDiscrete2Multinomial_sppmix(fit$genzs, fit$m)
plot_df <- data.frame(probability = as.vector(probs),
point = 1:n,
component = rep(1:fit$m, each =n))
ggplot(plot_df, aes(point, component, xend = point,
yend = component + 1,
col = probability,
xmin=1,xmax=n,
ymin=1,ymax=fit$m)) +
geom_segment(size = I(5)) +
ggplot2::scale_color_gradient(low = "white", high = "grey18",
guide = guide_colorbar(nbin = 100, barheight = 15)) +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::labs(x = "Point", y = "Component", colour = "Probability") +
ggplot2::scale_x_discrete(limits = labsx)+
ggplot2::scale_y_discrete(limits = labsy+.5,labels=labsy)+
add_title("Membership indicators", m = fit$m, n = n)
}
#' Plot the average intensity surface
#'
#' @description
#' This function calculates the intensity surface
#' at each posterior realization and then computes
#' the average for the intensity surface over a fine grid.
#' The result is a much smoother posterior estimator
#' of the intensity surface, which is not necessarily
#' the same as the surface of posterior means,
#' which is obtained by \code{\link{GetPMEst}}.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_avgsurf}
#'
#' @param fit An object the contains all posterior realizations, e.g., the return
#' value from \code{\link{est_mix_damcmc}} or \code{\link{est_mix_bdmcmc}}.
#' @param win An object of class \code{\link[spatstat]{owin}}.
#' @param LL Length of the side of the square grid.
#' The density or intensity is calculated on an L * L grid.
#' The larger this value is, the slower the calculation,
#' but the better the approximation.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param zlims The limits of the z axis. Defaults to [0,1.1*max(intensity)].
#' @param grayscale Logical flag to request a gray scale plot.
#' @param showplot Logical flag to request that the plot will be shown. Set to FALSE if you want to return the \code{\link[spatstat]{im.object}}, but do not want to produce the 3d plot.
#' @return An image as an object of class \code{\link[spatstat]{im.object}}.
#' @author Jiaxun Chen, Sakis Micheas
#' @seealso \code{\link{rnormmix}},
#' \code{\link{to_int_surf}},
#' \code{\link{rsppmix}},
#' \code{\link{est_mix_damcmc}},
#' \code{\link{plot_density}},
#' \code{\link[ggplot2]{ggtitle}},
#' \code{\link[ggplot2]{geom_point}},
#' \code{\link{plotmix_3d}},
#' \code{\link{GetPMEst}}
#' @examples
#' \donttest{
#' truemix <- rnormmix(m = 5, sig0 = .1, df = 5, xlim= c(-1, 1), ylim =c(0, 3))
#' trueintsurf=to_int_surf(truemix, lambda = 200, win =spatstat::owin( c(-1, 1),c(0, 3)))
#' plot(trueintsurf, main = "True Poisson intensity surface (mixture of normal components)")
#' pp1 <- rsppmix(trueintsurf)
#' # Run Data augmentation MCMC and get posterior realizations
#' postfit=est_mix_damcmc(pp1,m=5)
#' # Plot the average of the surfaces of the posterior realizations
#' avgsurf=plot_avgsurf(postfit, LL = 50)
#' p<-plot_density(as.data.frame(avgsurf))+ggplot2::ggtitle(
#' "Average surface of the posterior realization surfaces\n x denotes a true component mean")
#' #show the point pattern points
#' pp_df <- data.frame(pp1$x,pp1$y)
#' names(pp_df) <- c("x", "y")
#' p<-p + ggplot2::geom_point(data = pp_df,size=0.8)
#' #show the true means
#' mean_df <- data.frame(do.call(rbind, trueintsurf$mus))
#' names(mean_df) <- c("x", "y")
#' p + ggplot2::geom_point(data = mean_df, color = "red", shape = "x", size = 5)
#' #repeat for the contour plot
#' p<-plot_density(as.data.frame(avgsurf),contour = TRUE)+ggplot2::ggtitle(
#' "Average surface of the posterior realization surfaces\n x denotes a true component mean")
#' #show the point pattern points
#' pp_df <- data.frame(pp1$x,pp1$y)
#' names(pp_df) <- c("x", "y")
#' p<-p + ggplot2::geom_point(data = pp_df,size=0.8)
#' #show the true means
#' mean_df <- data.frame(do.call(rbind, trueintsurf$mus))
#' names(mean_df) <- c("x", "y")
#' p + ggplot2::geom_point(data = mean_df, color = "red", shape = "x", size = 5)
#' #plot the 3d surface again based on the returned object
#' plotmix_3d(avgsurf,title1 = paste("Average of", .9*postfit$L,
#' "posterior realizations of the intensity surface"))}
#'
#' @export
plot_avgsurf <- function(fit,
win = fit$data$window, LL = 100,
burnin = floor(fit$L/10), zlims = c(0, 0),
grayscale = FALSE,showplot =TRUE)
{
if(class(fit)=="bdmcmc_res")
{
tab=GetBDTable(fit,FALSE)
fit=GetBDCompfit(fit,tab$MAPcomp)$BDgens
burnin = floor(fit$L/ 10)
}
# cat("\nTo save all open rgl graphs use Save_AllOpenRglGraphs.\n")
# get limits
xlims <- c(win$xrange)
ylims <- c(win$yrange)
L <- fit$L
gridvals <- GetGrid_sppmix(LL,xlims,ylims)
xcoord <- as.vector(gridvals[[1]])
ycoord <- as.vector(gridvals[[2]])
#############ATTENTION, DO NOT TRANSPOSE
zcoord <- ApproxAvgPostIntensity(
fit$allgens_List, fit$genlamdas, LL, burnin,
xlims, ylims, fit$ApproxCompMass)
if(showplot)
{
title1 = paste("Average of",L-burnin,
"posterior realizations of the intensity surface")
jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red",
"#7F0000"))
if (grayscale == TRUE) {
col <- gray.colors(100,start = 1, end = 0)[findInterval(zcoord, seq(min(zcoord),
max(zcoord), length = 100))]
} else {
col <- jet.colors(100)[findInterval(zcoord, seq(min(zcoord),
max(zcoord), length = 100))]
}
# rgl::layout3d(matrix(1:2, 1, 2), widths = c(5, 1))
rgl::open3d(windowRect=c(50,50,1000,800),
zoom=1.2)
U=rgl::par3d("userMatrix")
rgl::par3d(userMatrix=
rgl::rotate3d(U,pi/4,0,0,1))
#find the highest z
maxz_height <- max(zcoord)
if (zlims[1] == 0 && zlims[2] == 0) {
zlims <- c(0, 1.1*maxz_height)
}
# Rangez=zlims[2]-zlims[1];
#############ATTENTION, DO NOT TRANSPOSE
rgl::persp3d(x = xcoord, y = ycoord, z = zcoord,
color = col, xlab="x",ylab="y",zlab="",
zlim=c(zlims[1]-0.01,zlims[2]),
box = FALSE, axes = FALSE)
rgl::axis3d('x')
rgl::axis3d('y')
rgl::axis3d('z-+', pos = c(xlims[1], ylims[2], 0))
rgl::rgl.lines(c(xlims[1], xlims[1]),
c(ylims[2], ylims[2]),
c(0,maxz_height),
color = 'black')
rgl::title3d(main=NULL)
rgl::text3d(xlims[2],ylims[2],
zlims[2],texts= title1)
rgl::bgplot3d(suppressWarnings(
fields::image.plot(legend.only = TRUE,
smallplot= c(.8,.82,0.05,.7),
zlim = zlims,
col = jet.colors(100))))
}
return( as.im(list(x=xcoord,
y=ycoord,
z=zcoord)))
}
#' Plot results from a DAMCMC fit
#'
#' @description
#' This function uses the posterior realizations
#' from a \code{\link{est_mix_damcmc}} call, to produce a plethora
#' of plots about the fitted Poisson point process with mixture intensity surface.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot.damcmc_res}
#'
#' @param x Object of class \code{damcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param ... Additional arguments for the S3 method.
#' @author Jiaxun Chen, Sakis Micheas, Yuchen Wang
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{PlotUSAStates}},
#' \code{\link{GetPMEst}},
#' \code{\link{check_labels}},
#' \code{\link{FixLS_da}},
#' \code{\link{plot_chains}}
#' @examples
#' \donttest{
#' fit <- est_mix_damcmc(pp = spatstat::redwood, m = 10)
#' plot(fit)
#' #Tornadoes
#' Tornfit=est_mix_damcmc(Tornadoes2011MO, m=5, L = 20000)
#' MAPsurf=GetMAPEst(Tornfit)
#' ret=PlotUSAStates(states=c('Iowa','Arkansas', 'Missouri','Illinois','Indiana',
#' 'Kentucky','Tennessee', 'Kansas','Nebraska','Texas','Oklahoma','Mississippi',
#' 'Alabama','Louisiana'),showcentroids=FALSE, shownames=TRUE, plotlevels = FALSE,
#' main="Tornadoes about MO, 2011", pp=Tornadoes2011MO, surf=MAPsurf,
#' boundarycolor="white", namescolor="white")
#' print(ret)
#' plot(Tornfit)
#' # get the intensity of posterior means
#' post_mean = GetPMEst(Tornfit)
#' # plot the estimated intensity surface
#' plot(post_mean)
#' #check labels
#' check_labels(Tornfit)
#' # If present then fix label switching, start with approx=TRUE
#' post_fixed = FixLS_da(Tornfit, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed,separate = FALSE)}
#'
#' @export
#' @method plot damcmc_res
plot.damcmc_res <- function(
x,burnin = floor(length(fit$allgens) / 10)
,...) {
fit<-x
post_mix <- GetPMEst(fit, burnin = burnin)
plot(post_mix)
print(plotmix_2d(post_mix, pattern = fit$data))
print(plot_ind(fit, burnin)+add_title("Posterior Means of the membership indicators", m = fit$m, n = fit$data$n))
print(plot_chains(fit, burnin,separate=FALSE))
invisible(NULL)
}
#' Plot results from a BDMCMC fit
#'
#' @description
#' This function uses the posterior realizations
#' from a \code{\link{est_mix_bdmcmc}} call, to produce a plethora
#' of plots about the fitted Poisson point process with mixture intensity surface.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot.bdmcmc_res}
#'
#' @param x Object of class \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param win An object of class \code{\link[spatstat]{owin}}.
#' @param LL Length of the side of the square grid.
#' The density or intensity is calculated on an L * L grid.
#' The larger this value is, the slower the calculation,
#' but the better the approximation as well as the smoother the resulting plots.
#' @param zlims The limits of the z axis. Defaults to [0,max(z)].
#' @param ... Additional arguments for the S3 method.
#' @details Unlike the corresponding output from DAMCMC (fixed number of components), the BDMCMC algorithm
#' allows us to obtain a distribution for the number of components
#' which can be thought of as a statistical inference procedure for model
#' selection. In particular, the Bayesian model average of all the realizations
#' is perhaps the best possible estimator of the Poisson intensity surface, however,
#' it can be very slow to compute for moderate number of iterations and maximum
#' number of components allowed.
#' @return A list with the following objects
#' \item{FreqTab}{Frequency table for the number of components.}
#' \item{Mapsurf}{The Maximum A Posteriori (MAP) Poisson intensity surface based on the corresponding posterior means (label switching might be present). The MAP is the mode of the distribution of the number of components.}
#' \item{BMA}{The Bayesian Model Average is returned only if we answer "Y" to request it. Alternatively, use function \code{\link{GetBMA}} to compute the BMA. This is an image, i.e., an object of class \code{\link[spatstat]{im.object}}.}
#' @author Jiaxun Chen, Sakis Micheas, Yuchen Wang
#' @seealso \code{\link{est_mix_bdmcmc}},
#' \code{\link{PlotUSAStates}},
#' \code{\link{plotmix_3d}},
#' \code{\link{plot_density}},
#' \code{\link{check_labels}},
#' \code{\link{FixLS_da}},
#' \code{\link{plot_chains}},
#' \code{\link{GetBMA}}
#' \code{\link{GetBDTable}},
#' @examples
#' \donttest{
#' fit <- est_mix_bdmcmc(pp = spatstat::redwood, m = 10)
#' plot(fit)
#' #Tornadoes
#' ret=PlotUSAStates(states=c('Iowa','Arkansas', 'Missouri','Illinois','Indiana','Kentucky',
#' 'Tennessee', 'Kansas','Nebraska','Texas','Oklahoma','Mississippi', 'Alabama','Louisiana'),
#' showcentroids=FALSE,shownames=TRUE, plotlevels = FALSE, main="Tornadoes about MO, 2011",
#' pp=Tornadoes2011MO)
#' print(ret)
#' Tornfit=est_mix_bdmcmc(pp = Tornadoes2011MO, m = 7)
#' TornResults=plot(Tornfit)#if we plot the Bayesian model average return it
#' TornResults
#' if(!is.null(TornResults$BMA)){
#' BMA_image=TornResults$BMA#must answer yes above or compute it using GetBMA
#' burnin=.1*Tornfit$L
#' title1 = paste("Bayesian model average of", Tornfit$L-burnin,"posterior realizations")
#' plotmix_3d(BMA_image,title1=title1)
#' plot_density(as.data.frame(BMA_image))+ggplot2::ggtitle(
#' "Bayesian Model Average Intensity")
#' plot_density(as.data.frame(BMA_image),TRUE)+ggplot2::ggtitle(
#' "Contours of the Bayesian Model Average Intensity")}
#' # Work with the MAP intensity
#' Mapsurf=TornResults$Mapsurf
#' plot(Mapsurf)
#' #retrieve realizations for the MAP number of components only
#' tab=GetBDTable(Tornfit)
#' MAPm=tab$MAPcomp
#' BDfitMAPcomp=GetBDCompfit(Tornfit,MAPm)
#' summary(BDfitMAPcomp)
#' summary(BDfitMAPcomp$BDgens)
#' plot(BDfitMAPcomp$BDsurf,main=paste(
#' "Poisson Mixture intensity surface, MAP # of components=",MAPm))
#' #check labels
#' check_labels(BDfitMAPcomp$BDgens)
#' # If present then fix label switching, start with approx=TRUE
#' post_fixed = FixLS_da(BDfitMAPcomp$BDgens, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed,separate = FALSE)
#' #this one works better
#' post_fixed = FixLS_da(BDfitMAPcomp$BDgens,approx=FALSE, plot_result = TRUE)
#' plot_chains(post_fixed)
#' plot_chains(post_fixed,separate = FALSE)}
#'
#' @export
#' @method plot bdmcmc_res
plot.bdmcmc_res <- function(x, win = fit$data$window,
burnin = floor(fit$L/10),
LL = 100, zlims = c(0, 0),...) {
fit<-x
L <- fit$L
xlims <- win$xrange
ylims <- win$yrange
BDtab=GetBDTable(fit,FALSE)
tab=BDtab$FreqTab
MAPcomp=BDtab$MAPcomp
plot_CompDist(fit)
mapcompgens=GetBDCompfit(fit,MAPcomp)
mapgens=mapcompgens$BDgens
BDnormmix=mapcompgens$BDnormmix
mean_lambda=mean(mapgens$genlamdas);
PP=fit$data
n=ncol(fit$genzs)
plot2dPP(PP,BDnormmix$mus,title1="Point pattern along with the mixture means for MAP m")
intsurfMAP=mapcompgens$BDsurf
print(summary(intsurfMAP))
print(plotmix_2d(intsurfMAP, PP)+add_title("MAP Poisson intensity surface along with the point pattern",lambda =intsurfMAP$lambda,m=intsurfMAP$m,n=PP$n,L=mapcompgens$BDgens$L))
print(plotmix_2d(intsurfMAP, PP,contour = TRUE)+add_title("MAP Poisson intensity surface contours along with the point pattern",lambda =intsurfMAP$lambda,m=intsurfMAP$m,n=PP$n,L=mapcompgens$BDgens$L))
print(plot_chains(fit, burnin,separate=FALSE))
dens_image=NULL
if(Get_User_Input_sppmix("Compute Bayesian model average (slow operation)?"))
{
dens_image=GetBMA(fit)
title1 = paste("Bayesian model average of",L-burnin,"posterior realizations")
plotmix_3d(dens_image,title1=title1)
}
if(is.null(dens_image))
return(list(FreqTab=tab,Mapsurf=intsurfMAP))
else
return(list(FreqTab=tab,Mapsurf=intsurfMAP,BMA=dens_image))
}
#' Plots for the number of components
#'
#' @description
#' The function produces two plots: the trace plot for the number of components based on
#' all realizations of a BDMCMC fit, and a barplot describing the distribution
#' of the number of components.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_CompDist}
#'
#' @param fit Object of class \code{bdmcmc_res}.
#' @param open_new_window Open new windows for the two plots.
#' @seealso \code{\link{est_mix_bdmcmc}}
#' @author Sakis Micheas
#' @examples
#' \donttest{
#' fitBD <- est_mix_bdmcmc(spatstat::redwood, m = 10)
#' plot_CompDist(fitBD)
#' CAfitBD=est_mix_bdmcmc(pp = CAQuakes2014.RichterOver3.0, m = 10)
#' plot_CompDist(CAfitBD)}
#'
#' @export
plot_CompDist=function(fit,open_new_window=FALSE)
{
openwin_sppmix(open_new_window)
BDtab=GetBDTable(fit)
tab=BDtab$FreqTab
oldpar=par(mai=c(1.2,1.2,1,1))
mp=barplot(tab,names.arg=1:fit$maxnumcomp,
xlab="Number of Components",ylab=paste("Iterations,",fit$L,"total"),
main="Distribution of the number of components"
,ylim=c(0,1.2*max(tab)))
for(i in 1:fit$maxnumcomp)
text(x=mp,y=tab+0.1*max(tab),labels=tab)
openwin_sppmix(open_new_window)
plot(fit$numcomp,ylim=c(1,fit$maxnumcomp),
xlab="Iteration",
ylab="Number of components",
type="l",main="Generated chain for the number of components")
suppressWarnings( par(oldpar))
}
#' Checking convergence: running means plot
#'
#' @description
#' This function produces a running mean plot, that is,
#' a plot of the iterations against the mean of
#' the draws up to each iteration. If the plot is not a near
#' constant line then convergence has not been achieved (e.g., label switching is present).
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_runmean}
#'
#' @param chain An \code{Lx1} vector containing the \code{L} posterior realizations.
#' @param open_new_window Open a new window for the plot.
#'
#' @author Sakis Micheas
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{rmixsurf}},
#' \code{\link{rsppmix}}
#' @examples
#' \donttest{
#' truemix_surf <- rmixsurf(m = 3, lambda=100, xlim = c(-3,3), ylim = c(-3,3))
#' plot(truemix_surf)
#' genPPP=rsppmix(intsurf = truemix_surf, truncate = FALSE)
#' fit <- est_mix_damcmc(pp = genPPP, m = 3)
#' plot_runmean(fit$genps[,1])
#' plot_runmean(fit$genps[,2])
#' plot_runmean(fit$genps[,3])}
#'
#' @export
plot_runmean<- function(chain,open_new_window = FALSE)
{
Iteration=PostMean=NULL
openwin_sppmix(check2open=open_new_window)
L=length(chain)
runmean=rep(0,L)
for(i in 1:L)
{
runmean[i]=mean(chain[1:i])
}
# plot(1:L,runmean,type="l")
plot_df <- data.frame(Iteration = 1:L,
PostMean =runmean)
plot_p <- ggplot(plot_df,
aes(Iteration, PostMean)) + geom_path() +
labs(x = "Iteration", y = "Posterior Mean") +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:L,expand=FALSE)+
ggplot2::theme(plot.margin=ggplot2::unit(c(1,1,1,0), "lines"))+
add_title("Running mean plot",L = L)
plot_p
}
#' Checking convergence: autocorrelation plot
#'
#' @description
#' This function can be used to assess
#' convergence by visualizing the
#' autocorrelations between the draws of the
#' Markov chain \code{chain}. The \code{lag k}
#' autocorrelation \code{rho_k} is the correlation between every draw
#' and its \code{kth} lag. We would expect the
#' \code{kth} lag autocorrelation to be
#' smaller as \code{k} increases (that is, the
#' 100th and 1000th draws should be less
#' correlated than the 100th and 105th draws).
#' For higher values of k we anticipate small
#' autocorrelation values, otherwise the chain
#' is not \code{mixing} well (in other words we do not
#' explore the parameter space adequately).
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_autocorr}
#'
#' @param chain An \code{Lx1} vector containing the \code{L} posterior realizations.
#' @param open_new_window Open a new window for the plot.
#' @param maxlag The maximum lag value to consider. Default is 100.
#' @author Sakis Micheas
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{rmixsurf}},
#' \code{\link{rsppmix}}
#' @examples
#' \donttest{
#' truemix_surf <- rmixsurf(m = 3, lambda=100, xlim = c(-3,3), ylim = c(-3,3))
#' plot(truemix_surf)
#' genPPP=rsppmix(intsurf = truemix_surf, truncate = FALSE)
#' fit <- est_mix_damcmc(pp = genPPP, m = 3)
#' plot_autocorr(fit$genps[,1])
#' plot_autocorr(fit$genps[,2])
#' plot_autocorr(fit$genps[,3])}
#'
#' @export
plot_autocorr<- function(chain,
open_new_window = FALSE,maxlag=100)
{
Autocorr=NULL
openwin_sppmix(check2open=open_new_window)
L=length(chain)
if(L<=maxlag)
stop("Reduce lag or increase number of iterations.")
autocorr=rep(0,maxlag)
xbar=mean(chain)
xvar=(L-1)*var(chain)
for(k in 1:maxlag)
{
autocorr[k]=sum((chain[1:(L-k)]-xbar)*
(chain[(k+1):L]-xbar))/xvar
}
plot_df <- data.frame(lag = 1:L,
Autocorr =autocorr)
plot_p <- ggplot(plot_df,
aes(lag, Autocorr)) + geom_path() +
labs(x = "Lag", y = "Autocorrelation") +
ggplot2::theme_classic() +
ggplot2::theme(panel.border = ggplot2::element_rect(fill = NA, size = 1)) +
ggplot2::coord_cartesian(xlim =1:maxlag,expand=FALSE)+
ggplot2::theme(plot.margin=ggplot2::unit(c(1,1,1,0), "lines"))+
add_title("Autocorrelation plot")
plot_p
}
#' Checking convergence visually
#'
#' @description
#' Based on a `damcmc_res` object, this function will produce
#' many graphs to help assess convergence visually,
#' including running mean plots and
#' autocorrelation plots for all the parameters. This function calls
#' \code{\link{plot_runmean}} and \code{\link{plot_autocorr}} for all parameters
#' so we do not have to it individually.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #plot_convdiags}
#'
#' @param fit Object of class \code{damcmc_res} or \code{bdmcmc_res}.
#' @param burnin Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
#' @param open_new_window Open a new window for the plot.
#' @param maxlag The maximum lag value to consider. Default is 100.
#' @author Sakis Micheas
#' @seealso \code{\link{est_mix_damcmc}},
#' \code{\link{rmixsurf}},
#' \code{\link{plot_runmean}},
#' \code{\link{plot_autocorr}},
#' \code{\link{rsppmix}}
#' @examples
#' \donttest{
#' truemix_surf <- rmixsurf(m = 3, lambda=100, xlim = c(-3,3), ylim = c(-3,3))
#' plot(truemix_surf)
#' genPPP=rsppmix(intsurf = truemix_surf, truncate = FALSE)
#' fit = est_mix_damcmc(pp = genPPP, m = 3)
#' plot_convdiags(fit)}
#'
#' @export
plot_convdiags<- function(fit,burnin = floor(fit$L / 10),
open_new_window = FALSE,maxlag=100)
{
fit <- drop_realization(fit, burnin)
open_new_window1=open_new_window
maxlag1=maxlag
if(fit$L<=maxlag1)
stop("Reduce lag or increase number of iterations.")
#probs first
for(j in 1:fit$m)
{
print(plot_runmean(fit$genps[,j],open_new_window=open_new_window1)+ggtitle(paste("Component probability",j)))
print(plot_autocorr(fit$genps[,j],
open_new_window=open_new_window1, maxlag =maxlag1)+ggtitle(paste("Component probability",j)))
}
#mu x-coords
for(j in 1:fit$m)
{
print(plot_runmean(fit$genmus[j,1,],open_new_window=open_new_window1)+ggtitle(paste("x-coord for mu, component",j)))
print(plot_autocorr(fit$genmus[j,1,],
open_new_window=open_new_window1, maxlag =maxlag1)+ggtitle(paste("x-coord for mu, component",j)))
}
#mu y-coords
for(j in 1:fit$m)
{
print(plot_runmean(fit$genmus[j,2,],open_new_window=open_new_window1)+ggtitle(paste("y-coord for mu, component",j)))
print(plot_autocorr(fit$genmus[j,2,],
open_new_window=open_new_window1, maxlag =maxlag1)+ggtitle(paste("y-coord for mu, component",j)))
}
}
#' Closes all open plots
#'
#' @description
#' The function closes all Rgl plots, as well as, any graphics devices.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #Plots_off}
#'
#' @seealso \code{\link{Save_AllOpenRglGraphs}}
#'
#' @author Sakis Micheas
#' @examples
#' \donttest{
#' Plots_off()}
#'
#' @export
Plots_off<- function()
{
graphics.off()
while (rgl::rgl.cur()>0)
{
rgl::rgl.close()
}
}
#' Saves RGL plots
#'
#' @description
#' The function saves all open RGL plots (3d plots) to the specified directory
#' and using a template name.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #Save_AllOpenRglGraphs}
#'
#' @param dir1 Directory to save the plots.
#' @param filename1 Template filename. Each open plot is saved as filename1.png, filename2.png, and so forth.
#' @author Sakis Micheas
#' @seealso \code{\link{Plots_off}}
#' @examples
#' \donttest{
#' # use a temporary directory to save the plots
#' Save_AllOpenRglGraphs(dir1=tempdir())}
#'
#' @export
Save_AllOpenRglGraphs<- function(
dir1,filename1="RglGraph")
{
if(missing(dir1))
stop("Please specify a directory wherein to save the RGL plots.")
# graphics.off()
count=1;
while (rgl::rgl.cur()>0)
{
rgl::rgl.snapshot(paste(filename1,count,".png",sep = ""))
rgl::rgl.close()
count=count+1;
}
if (count==2)
cat(paste("Saved 1 plot in",dir1))
else
cat(paste("Saved",count-1,"plots in",dir1))
}
Try the sppmix package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.