plot.bdmcmc_res: Plot results from a BDMCMC fit

In sppmix: Modeling Spatial Poisson and Related Point Processes

Description

This function uses the posterior realizations from a est_mix_bdmcmc call, to produce a plethora of plots about the fitted Poisson point process with mixture intensity surface.

For examples see

http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html#plot.bdmcmc_res

Usage

## S3 method for class 'bdmcmc_res'
plot(x, win = fit$data$window,
  burnin = floor(fit$L/10), LL = 100, zlims = c(0, 0), ...)

Arguments

x	Object of class bdmcmc_res.
win	An object of class owin.
burnin	Number of initial realizations to discard. By default, it is 1/10 of the total number of iterations.
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.
zlims	The limits of the z axis. Defaults to [0,max(z)].
...	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.

Value

A list with the following objects

FreqTab	Frequency table for the number of components.
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.
BMA	The Bayesian Model Average is returned only if we answer "Y" to request it. Alternatively, use function GetBMA to compute the BMA. This is an image, i.e., an object of class im.object.

Author(s)

Jiaxun Chen, Sakis Micheas, Yuchen Wang

See Also

est_mix_bdmcmc, PlotUSAStates, plotmix_3d, plot_density, check_labels, FixLS_da, plot_chains, GetBMA GetBDTable,

Examples

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)

sppmix documentation built on Jan. 13, 2021, 10:04 p.m.