functions on kernel density estimation in pointpattern analysis

Dear Roger Bivand,Tim Keitt, and Dan Putler,
  Thanks for your answers. I have tried density.ppp(spatstat) and
kernel2d(splancs), but the results are not very satisfied. I think there
should be a higher density in the blue part of the map in the attachment.
  My dataset has been put in the attachment and programs have been pasted in
the following, so that u can use and check it.
   I want to do the work like "non-parametric estimation of a spatially
varying intensity" in Diggle's book(2003.P.116-121).
  BTW, i'm not familiar with locfit,would u please also check it using
locfit? Thanks very much.
###############################################################################
                           Kernel density estimation--spatstat
################################################################################
library(sp)
library(foreign)
library(maptools)
library(mgcv)
library(spatstat)

guichi<-readShapePoly("d:/deleting/kernel/kernel/guichi2.shp")
W <-as(as(guichi, "SpatialPolygons"), "owin")  #boundary polygons

cases<-coordinates(readShapePoints("d:/deleting/kernel/kernel/cases.shp"))
#points
colnames(cases)<-c("x","y")
cases[1:2,]

#plot(W);points(cases)

pointcase <- ppp(cases[,1], cases[,2], window=W)  #generate the ppp object

kdensity<-density.ppp(pointcase, 0.05)
plot(kdensity)
*Q:*there are almost the same density in the whole area,but in fact it may
have a higher density in the blue part of the attached map? I think the
problem may the inappropriate value of sigma, how to determine its value?

################################################################################
                                 Kernel density estimation--splancs
################################################################################
library(sp)
library(splancs)
library(foreign)
library(maptools)

case<-readShapePoints("d:/deleting/kernel/kernel/cases.shp")
guichi<-readShapePoly("d:/deleting/kernel/kernel/guichi2.shp")

#Conversion
case_pts <- coordinates(case)
case <- as.points(case_pts)
splancs_poly <-
getPolygonCoordsSlot(getPolygonsPolygonsSlot(getSpPpolygonsSlot(guichi)[[1]])[[1]])

#to unpack the coordinates of the points and the single ring boundary
polymap(splancs_poly,xlab="x(?)",ylab="y(?)")
pointmap(case_pts, add=TRUE)

m<-mse2d(case,splancs_poly,nsmse=1000,range=5)  #plots the estimated mean
square error as a function of h0
plot(m$h[290:1000],m$mse[290:1000],type="l")

n<-which(m$mse==min(m$mse))
h0<-m$h[n]

#smooth variation
smooth<-kernel2d(case, splancs_poly, h0=h0, nx=100, ny=100)
polymap(splancs_poly)  #sets the axes correctly and draws the polygon
image(smooth,add=T)   #the smoothed image is superimposed
polymap(splancs_poly,add=T) #redrawn the polymap in order not to be obsured
by smooth image
*Q*:The result is still not satisfied,there must be something wrong with my
programs.

functions on kernel density estimation in pointpattern analysis

Thread (8 messages)