Install any packages you might be missing:
library(spatstat)
library(splancs)
# tree data
data(finpines)
#get info
print(finpines)
## Marked planar point pattern: 126 points
## Mark variables: diameter, height
## window: rectangle = [-5, 5] x [-8, 2] metres
intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern: 126 points
## Average intensity 1.26 points per square metre
##
## Coordinates are given to 13 decimal places
##
## Mark variables: diameter, height
## Summary:
## diameter height
## Min. :0.000 Min. :0.800
## 1st Qu.:1.000 1st Qu.:1.825
## Median :2.000 Median :2.850
## Mean :2.532 Mean :2.828
## 3rd Qu.:3.000 3rd Qu.:3.600
## Max. :7.000 Max. :5.400
##
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
#creates point pattern data format from the x and y coordinates
finpines.pts<-as.points(finpines$x,finpines$y)
plot(finpines.pts, pch=19, xlab="X", ylab="Y",main="Finish Pines Locations")
#divide into quadrants and count number of points in the quadrant
Q <- quadratcount(finpines, nx = 4, ny = 3)
Q
## x
## y [-5,-2.5) [-2.5,0) [0,2.5) [2.5,5]
## [-1.33,2] 10 8 11 21
## [-4.67,-1.33) 8 7 9 3
## [-8,-4.67) 12 11 22 4
summary(Q)
## Number of cases in table: 126
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 20.215, df = 6, p-value = 0.002535
Q_dat<-data.frame(Q)
var(Q_dat$Freq)
## [1] 33.72727
mean(Q_dat$Freq)
## [1] 10.5
plot(finpines.pts,pch=19, main="Finpines", xlab="X",ylab="Y")
plot(Q, add = TRUE, cex = 2)
# Finpines data
sp_point <- matrix(NA, nrow=length(finpines$x),ncol=2)
sp_point[,1] <- finpines$x
sp_point[,2] <- finpines$y
colnames(sp_point) <- c("x","y")
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="X",ylab="Y",cex=.5)
intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern: 126 points
## Average intensity 1.26 points per square metre
##
## Coordinates are given to 13 decimal places
##
## Mark variables: diameter, height
## Summary:
## diameter height
## Min. :0.000 Min. :0.800
## 1st Qu.:1.000 1st Qu.:1.825
## Median :2.000 Median :2.850
## Mean :2.532 Mean :2.828
## 3rd Qu.:3.000 3rd Qu.:3.600
## Max. :7.000 Max. :5.400
##
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
# Random points
u.x <- runif(n=nrow(sp_point), min=bbox(sp_point)[1,1], max=bbox(sp_point)[1,2])
u.y <- runif(n=nrow(sp_point), min=bbox(sp_point)[2,1], max=bbox(sp_point)[2,2])
# Regular points
r.x <- seq(from=min(sp_point[,1]),to=max(sp_point[,1]),length=sqrt(nrow(sp_point)))
r.y <- seq(from=min(sp_point[,2]),to=max(sp_point[,2]),length=sqrt(nrow(sp_point)))
r.x <- jitter(rep(r.x,length(r.x)),.001)
r.y <- jitter(rep(r.y,each=length(r.y)),.001)
# Plot Finpines, uniform, and regular points
par(mfrow=c(1,3),mar=c(4,4,1.5,0.5))
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="X",ylab="Y",cex=.5,pch=19)
plot(x=u.x,y=u.y,main="Random Points", xlab="X",ylab="Y",cex=.5,pch=19)
plot(x=r.x,y=r.y,main="Regular Points", xlab="X",ylab="Y",cex=.5,pch=19)
# CSR from random uniform point process
## Generates a CSR point process conditionally on a fixed number N of points
## Let N = 100
par(mfrow=c(2,2),mar=c(0.4,0.4,1,0.4))
pp.100.unif<-runifpoint(100,win=owin(c(0,1),c(0,1)))
pp.100.unif
## Planar point pattern: 100 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unif)
## [1] 100
plot(pp.100.unif, pch=19, cex=0.6,main="")
title('100 uniform points')
# repeat a few times to see different realizations of the process
# HPP with intensity 100 in the unit square
pp.100.unit <- rpoispp(100,win=owin(c(0,1),c(0,1)))
pp.100.unit
## Planar point pattern: 109 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unit)
## [1] 109
plot(pp.100.unit,main="")
title('HPP intensity = 100, unit square')
# repeat a few times to see different realizations of the process (including different numbers of points)
# HPP with intensity 1 in a 10 x 10 square
pp.1.10 <- rpoispp(1, win=owin(c(0,10),c(0,10)))
par(mfrow=c(1,2))
plot(pp.100.unit,main="")
title('intensity = 100, unit square')
plot(pp.1.10,main="")
title('intensity = 1, 10 x 10 square')
# HPP with intensity 200 in the unit square
par(mfrow=c(1,1))
pp.200.unit <- rpoispp(200)
plot(pp.200.unit)
# HPP with intensity 100 in the unit square
side.square <-1
pp.100.unit <- rpoispp(10000,win=owin(c(0,side.square),c(0,side.square)))
points.tot <- npoints(pp.100.unit)
plot(pp.100.unit,main="")
# inscribe a circle within the square
radius.circle <- side.square/2
win.circle <- disc(radius.circle, c(radius.circle,radius.circle))
pp.incircle <- pp.100.unit[win.circle]
points.incircle <- npoints(pp.incircle)
summary(pp.incircle)
## Planar point pattern: 7943 points
## Average intensity 10117.4 points per square unit
##
## Coordinates are given to 8 decimal places
##
## Window: polygonal boundary
## single connected closed polygon with 128 vertices
## enclosing rectangle: [0, 1] x [0, 1] units
## Window area = 0.785083 square units
## Fraction of frame area: 0.785
plot(pp.incircle)
# calculate pi using ratio of areas and number of points
piEst <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
piEst
## [1] 3.148236
Function to simulate pi over many iterations
pi.simulation2 = function(iter, lambda, side.square){
piEst.list = rep(0,iter)
points.incircle.list = rep(0,iter)
# inscribe a circle within the square
radius.circle <- side.square/2
win.circle <- disc(radius.circle, c(radius.circle,radius.circle))
for(i in 1:iter){
# HPP with intensity 100 in the unit square
pp.square <- rpoispp(lambda,win=owin(c(0,side.square),c(0,side.square)))
points.tot <- npoints(pp.square)
# find points within the circle
pp.incircle <- pp.square[win.circle]
points.incircle <- npoints(pp.incircle)
# calculate pi using ratio of areas and number of points
pihat <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
piEst.list[i] = pihat
}
return( piEst.list)
}
iter = 1000
lambda = 1000
side.square = .5
est = pi.simulation2(iter,lambda, side.square)
mean(est)
## [1] 3.146699
var(est)
## [1] 0.01027661
# we plot the histogram and see it appears to be binomial. why?
hist(est,nclass = 15)
# create polygon (square) area for finpines data
minfx <- min(finpines$x)
maxfx <- max(finpines$x)
minfy <- min(finpines$y)
maxfy <- max(finpines$y)
polyfin <- as.points(c(minfx,maxfx,maxfx,minfx),c(minfy,minfy,maxfy,maxfy))
UL.khat.fin <- Kenv.csr(length(finpines$x), polyfin, nsim=99, seq(0.05,3.5,0.15))
## Doing simulation 1
## Doing simulation 2
## Doing simulation 3
## Doing simulation 4
## Doing simulation 5
## Doing simulation 6
## Doing simulation 7
## Doing simulation 8
## Doing simulation 9
## Doing simulation 10
## Doing simulation 11
## Doing simulation 12
## Doing simulation 13
## Doing simulation 14
## Doing simulation 15
## Doing simulation 16
## Doing simulation 17
## Doing simulation 18
## Doing simulation 19
## Doing simulation 20
## Doing simulation 21
## Doing simulation 22
## Doing simulation 23
## Doing simulation 24
## Doing simulation 25
## Doing simulation 26
## Doing simulation 27
## Doing simulation 28
## Doing simulation 29
## Doing simulation 30
## Doing simulation 31
## Doing simulation 32
## Doing simulation 33
## Doing simulation 34
## Doing simulation 35
## Doing simulation 36
## Doing simulation 37
## Doing simulation 38
## Doing simulation 39
## Doing simulation 40
## Doing simulation 41
## Doing simulation 42
## Doing simulation 43
## Doing simulation 44
## Doing simulation 45
## Doing simulation 46
## Doing simulation 47
## Doing simulation 48
## Doing simulation 49
## Doing simulation 50
## Doing simulation 51
## Doing simulation 52
## Doing simulation 53
## Doing simulation 54
## Doing simulation 55
## Doing simulation 56
## Doing simulation 57
## Doing simulation 58
## Doing simulation 59
## Doing simulation 60
## Doing simulation 61
## Doing simulation 62
## Doing simulation 63
## Doing simulation 64
## Doing simulation 65
## Doing simulation 66
## Doing simulation 67
## Doing simulation 68
## Doing simulation 69
## Doing simulation 70
## Doing simulation 71
## Doing simulation 72
## Doing simulation 73
## Doing simulation 74
## Doing simulation 75
## Doing simulation 76
## Doing simulation 77
## Doing simulation 78
## Doing simulation 79
## Doing simulation 80
## Doing simulation 81
## Doing simulation 82
## Doing simulation 83
## Doing simulation 84
## Doing simulation 85
## Doing simulation 86
## Doing simulation 87
## Doing simulation 88
## Doing simulation 89
## Doing simulation 90
## Doing simulation 91
## Doing simulation 92
## Doing simulation 93
## Doing simulation 94
## Doing simulation 95
## Doing simulation 96
## Doing simulation 97
## Doing simulation 98
## Doing simulation 99
khat.fin<-khat(finpines.pts, polyfin, seq(0.05,3.5,0.15))
# plot of Khat-pi t^2 from data
plot(seq(0.05,3.5,0.15), khat.fin-pi*seq(0.05,3.5,0.15)^2, type="l", xlab="Distance", ylab="Estimated K-pi h^2")
# plot upper bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$upper-pi*seq(0.05,3.5,0.15)^2, lty=2)
# plot lower bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$lower-pi*seq(0.05,3.5,0.15)^2, lty=2)
# L functions
l<-function(k,h)
{sqrt(k/pi)-h}
# plot Lhat from data
plot(seq(0.05,3.5,0.15), l(khat.fin, seq(0.05,3.5,0.15)), type="l", xlab="Distance", ylab="Estimated L",ylim=c(-.1,.2))
# plot upper bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$upper,seq(0.05,3.5,0.15)), lty=2)
# plot lower bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$lower,seq(0.05,3.5,0.15)), lty=2)
K<-Kest(finpines)
plot(K)
plot(K, cbind(r, sqrt(iso/pi)) ~ r)
plot(K, cbind(trans,iso,border) - theo ~ r)
plot(envelope(finpines,Kest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
# G test: Nearest neighbor distribution test
Gtest<-Gest(finpines)
plot(envelope(finpines,Gest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
# F test: Point to nearest event distribution test
Ftest<-Fest(finpines)
plot(envelope(finpines,Fest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
# J test: J(r) = (1-G(r)/(1-F(r)))
Jtest<-Jest(finpines)
plot(envelope(finpines,Jest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.