CompareSurfs: Quantify the difference between two surfaces
In sppmix: Modeling Spatial Poisson and Related Point Processes

Description Usage Arguments Details Value Author(s) See Also Examples

This function can be used to compare two intensity surfaces. In particular, if we have the true surface and an estimator of the truth, then we can assess how well the estimate fits the surface, i.e., how close are the two surfaces. Now if we have two estimates of the true surface then the estimate with the smallest measure fits the truth better. We can also compare two estimating surfaces this way.

For examples see

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

1	CompareSurfs(surf1, surf2, LL = 100, truncate = FALSE)

`surf1, surf2`	Either IPPP intensity surfaces (objects of type `intensity_surface`) or images (objects of type `im`) that represent intensity surfaces.
`LL`	Length of the side of the square grid. The intensities are calculated on an L * L grid. The larger this value is, the slower the calculation, but the better the approximation to the difference between the two surfaces.
`truncate`	Requests to truncate the components of the mixture intensities to have all their mass within the observation window.

Since the two surfaces passed to the function can be represented as a 2d intensity surface, any measure between two images can be used for comparison purposes, provided that the window is the same.

If the two windows are different the function will choose the largest one and compare the two surfaces in there.

Returns a list containing all distances computed and the window of comparison, an object of class owin.

Sakis Micheas

rmixsurf, owin, rsppmix, est_mix_bdmcmc, drop_realization, plotmix_3d, plot_avgsurf, GetBMA,

#compare two surfaces first
mixsurf1 = rmixsurf(m = 5, sig0 = .1, df = 5,xlim= c(-1,4), ylim = c(-2,1), rand_m = FALSE)
mixsurf2 = rmixsurf(m = 8, sig0 = .1, df = 5,xlim= c(-4,3), ylim = c(-1,2), rand_m = FALSE)
comp=CompareSurfs(mixsurf1,mixsurf2)
plot(mixsurf1,main = "First IPPP",win=comp$win)
plot(mixsurf2,main = "Second IPPP",win=comp$win)
# now we compare intensity surfaces and image objects that represent intensity surfaces
truemixsurf = rmixsurf(m = 5,xlim= c(-2,2), ylim = c(-2,2))
plot(truemixsurf,main="True IPPP surface")
#get a point pattern
genpp = rsppmix(truemixsurf,truncate=FALSE)
# Run BDMCMC and get posterior realizations
postfit=est_mix_bdmcmc(genpp,m=10,L=30000)
postfit=drop_realization(postfit,.1*postfit$L) #apply burnin
BMA=GetBMA(postfit,burnin=0)
title1 = paste("Bayesian model average of",postfit$L,"posterior realizations")
plotmix_3d(BMA,title1=title1)
comp=CompareSurfs(truemixsurf,BMA,LL=100)
#We compare the average surface and the truth for many cases below. If we pass images, we
#make sure it has the same dimensions or we force it to the same value by setting LL=100.
#We retrieve the average surfaces corresponding to MAP-1, MAP and MAP+1 components and
#compare them against the truth.
#First retrieve the frequency table and MAP estimate for number of components
BDtab=GetBDTable(postfit)
BDtab
MAPm=BDtab$MAPcomp
BDfitMAPcomp_minus1=GetBDCompfit(postfit,MAPm-1,burnin=0)
avgsurfMAPcomp_minus1=plot_avgsurf(BDfitMAPcomp_minus1$BDgens, LL = 100,burnin=0)
comp=CompareSurfs(truemixsurf,avgsurfMAPcomp_minus1,LL=100)
BDfitMAPcomp=GetBDCompfit(postfit,MAPm,burnin=0)
avgsurfMAPcomp=plot_avgsurf(BDfitMAPcomp$BDgens, LL = 100,burnin=0)
comp=CompareSurfs(truemixsurf,avgsurfMAPcomp,LL=100)
BDfitMAPcomp_plus1=GetBDCompfit(postfit,MAPm+1,burnin=0)
avgsurfMAPcomp_plus1=plot_avgsurf(BDfitMAPcomp_plus1$BDgens, LL = 100,burnin=0)
comp=CompareSurfs(truemixsurf,avgsurfMAPcomp_plus1,LL=100)