R/approximations.R In sppmix: Modeling Spatial Poisson and Related Point Processes

Documented in approx_normmixdnormmixGetBMA

```#' Approximate the masses of bivariate normal mixture components
#'
#' @description
#' Calculates the mass of the density of each component of a normal mixture over a given window.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #approx_normmix}
#'
#' @param mix An object of class \code{normmix}
#' @param xlim,ylim Vectors defining the x-y integration limits. A mixture component mass is estimated within
#' this window.
#'
#' @return A numerical vector with elements corresponding to the mass of each component
#' within the window. These values are required when we
#' use truncation over an observation window in order to handle edge effects.
#' @author Jiaxun Chen, Yuchen Wang
#' @examples
#' \donttest{
#' truemix = normmix(ps = c(.3, .7), mus = list(c(0.2, 0.2),
#'  c(.8, .8)), sigmas =list(.01*diag(2),.01*diag(2)))
#' approx_normmix(truemix,xlim= c(-2, 2), ylim = c(-2, 2))}
#'
#' @export
approx_normmix <- function(mix, xlim = c(0, 1), ylim = c(0, 1)) {

if (!is.normmix(mix)) {
stop("mix must be of class normmix or intensity surface")
}

approx <- numeric(mix\$m)

for(k in 1:mix\$m)
{
approx[k] <- mvtnorm::pmvnorm(lower = c(xlim[1], ylim[1]),  upper = c(xlim[2], ylim[2]),
mean = mix\$mus[[k]],sigma = mix\$sigmas[[k]])
#    approx[k]=ApproxCompMass_sppmix(xlim,ylim,mix\$mus[[k]],mix\$sigmas[[k]])
}
return(approx)
}

#' Calculate the density or intensity of a
#' normal mixture over a fine grid
#'
#' @description
#' When a \code{normmix} object is given, this
#' function calculates the mixture density over
#' a fine grid for the given window. When an \code{intensity_surface} object
#' is given, the function multiplies the density with
#'  the intensity surface parameter \code{lambda}, and returns the Poisson
#'  mixture intensity function over the grid. Used for plotting.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #dnormmix}
#'
#' @param mix An object of class \code{normmix} or \code{intensity_surface}
#' @param xlim,ylim Vectors defining the x-y integration limits. A mixture component mass is estimated within
#' this window.
#' @param L 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 truncate Requests to truncate the components
#' of the mixture intensity to have all their mass
#' within the given x-y limits.
#'
#' @return An object of class \code{\link[spatstat]{im}}. This is a pixel image
#'  on a grid with values corresponding to the density (or intensity surface) at that location.
#'
#'
#' @author Jiaxun Chen, Sakis Micheas, Yuchen Wang
#' @examples
#' \donttest{
#' truemix <- rnormmix(m = 3, sig0 = .1, df = 5,xlim= c(0, 5),
#'  ylim = c(0, 5))
#' normdens=dnormmix(truemix,xlim= c(0, 5), ylim = c(0, 5))
#' #2d plots
#' plot_density(as.data.frame(normdens))+ ggplot2::ggtitle(
#'  "2d mixture density plot\nWindow=[0,5]x[0,5]")
#' plot_density(as.data.frame(normdens),TRUE)+ ggplot2::ggtitle(
#'  "2d mixture contour plot\nWindow=[0,5]x[0,5]")
#' #3d plot
#' plotmix_3d(normdens)
#' #Now build an intensity surface based on the normal mixture
#' intsurf=to_int_surf(truemix,lambda = 100, win =
#'  spatstat::owin( c(0, 5),c(0, 5)))
#' intsurfdens=dnormmix(intsurf,xlim= c(0, 5), ylim = c(0, 5))
#' plot_density(as.data.frame(intsurfdens))+ ggplot2::ggtitle(
#'  "2d mixture intensity plot\nWindow=[0,5]x[0,5]")
#' plot_density(as.data.frame(intsurfdens),TRUE)+ ggplot2::ggtitle(
#'  "2d mixture intensity contour plot\nWindow=[0,5]x[0,5]")
#' plotmix_3d(normdens)#3d plot
#' #For an intensity surface object we use these functions instead
#' plotmix_2d(intsurf)
#' plot(intsurf)}
#'
#' @export
dnormmix <- function(mix, xlim = c(0, 1), ylim = c(0, 1), L = 128,
truncate = TRUE) {

if (!is.normmix(mix)) {
stop("mix must be of class normmix or intensity surface")
}

if (is.intensity_surface(mix)) {
xlim <- mix\$window\$xrange
ylim <- mix\$window\$yrange
}

x <- seq(xlim[1], xlim[2], length.out = L)
y <- seq(ylim[1], ylim[2], length.out = L)

approx=rep(1,mix\$m)
if (truncate) approx <- approx_normmix(mix, xlim, ylim)

mixlist=MakeMixtureListFromNormMix(mix)
est=matrix(0,L,L)
for(i in 1:L)
for(j in 1:L)
est[j,i]=densNormMixatx_sppmix(c(x[i],y[j]),mixlist,approx)

if (is.intensity_surface(mix)) {
spatstat::im(est * mix\$lambda, x, y)
} else {
spatstat::im(est, x, y)
}
}

#' Compute the Bayesian Model average
#'
#' @description
#' This function uses the posterior realizations
#' from a \code{\link{est_mix_bdmcmc}} call, to compute the
#' Bayesian Model Average across different number of components
#' and returns the fitted Poisson point process with mixture of normals intensity surface.
#'
#' For examples see
#'
#' \url{http://faculty.missouri.edu/~micheasa/sppmix/sppmix_all_examples.html
#' #GetBMA}
#'
#' @param fit Object of class \code{bdmcmc_res}.
#' @inheritParams plot.bdmcmc_res
#' @return An image as an object of class \code{\link[spatstat]{im.object}}.
#' @author Sakis Micheas
#' @examples
#' \donttest{
#' fit=est_mix_bdmcmc(pp = spatstat::redwood, m = 5)
#' BMA=GetBMA(fit)
#' burnin=.1*fit\$L
#' title1 = paste("Bayesian model average of",fit\$L-burnin,"posterior realizations")
#' plotmix_3d(BMA,title1=title1)
#' plot_density(as.data.frame(BMA))+ggplot2::ggtitle("Bayesian model average intensity surface")
#' plot_density(as.data.frame(BMA),TRUE)+ggplot2::ggtitle(
#'  "Contours of the Bayesian model average intensity surface")}
#'
#' @export
GetBMA<- function(fit, win = fit\$data\$window,
burnin = fit\$L/10,
LL = 100, zlims = c(0, 0))
{
# Start the clock!
ptm1 <- proc.time()
L <- fit\$L
xlims <- win\$xrange
ylims <- win\$yrange
#  cat("\nTo save all open rgl graphs use Save_AllOpenRglGraphs.\n")
distr_numcomp <- GetCompDistr_sppmix(fit\$numcomp[(burnin+1):L],fit\$maxnumcomp)
zcoord <- ApproxBayesianModelAvgIntensity_sppmix(
fit\$allgens_List[(burnin+1):L],
fit\$genlamdas[(burnin+1):L],
fit\$numcomp[(burnin+1):L],
distr_numcomp,1,fit\$maxnumcomp,LL, xlims, ylims,
fit\$ApproxCompMass[(burnin+1):L,])
# Stop the clock
ptm<-proc.time() - ptm1
cat(paste("\nComputation time of the Bayesian Model Average in seconds:",ptm[[1]],"\n"))
#find the highest z
maxz_height <- max(zcoord)
if (zlims[1] == 0 && zlims[2] == 0) {
zlims <- c(0, 1.1*maxz_height)
}
gridvals=GetGrid_sppmix(LL,xlims,ylims);
xcoord=as.vector(gridvals[[1]]);
ycoord=as.vector(gridvals[[2]]);
dens_image=spatstat::as.im(list(x=xcoord,y=ycoord,z=zcoord))
return(dens_image)
}
```

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.