x<-runif(100)
y<-runif(100)

par(pty="s")
plot(x,y,pch=16)

dist<-function(x1, y1, x2, y2)  sqrt((x2  x1)^2 + (y2  y1)^2)

r<-numeric(100)
nn<-numeric(100)
d<-numeric(100)
for (i in 1:100) {
d<-0
for (k in 1:100) d[k]<-dist(x[i],y[i],x[k],y[k])
r[i]<-min(d[-i])
nn[i]<-which(d==min(d[-i]))
}


for (i in 1:100) lines(c(x[i],x[nn[i]]),c(y[i],y[nn[i]]))

mean(r)

distances<-numeric(100)
for(i in 1:100) distances[i]<-dist(x[1],y[1],x[i],y[i])

sum(distances<0.1)  1

dd<-numeric(10000)
dd<-matrix(dd,nrow=100)

for (j in 1:100) {for(i in 1:100) dd[j,i]<-dist(x[j],y[j],x[i],y[i])}

sum(dd[1,]<0.1)-1

sum(dd<0.1)-100

d<-seq(0.01,1,0.01)

count<-numeric(100)

for (i in 1:100) count[i]<-sum(dd<d[i])-100

K<-count/10000

plot(d,K,type="l")
lines(d,pi*d^2)

# Check the location of the spatial library on your system

library(spatial)
pines<-ppinit("pines.dat")

par(mfrow=c(1,2),pty="s")

plot(pines,pch=16)
plot(Kfn(pines,5),type="s",xlab="distance",ylab="L(t)")

lims<-Kenvl(5,100,Psim(71))
lines(lims$x,lims$lower,lty=2)
lines(lims$x,lims$upper,lty=2)

attach(trees)
names(trees)

aa<-data.frame(xl=0,xu=1,yl=0,yu=1) 

map<-NULL
map$y<-y
map$x<-x
map$area<-aa

par(mfrow=c(1,2),pty="s")
plot(x,y,pch=16)
plot(Kfn(map,1),type="s",xlab="distance",ylab="L(t)")

lims<-Kenvl(1,100,Psim(257))
lines(lims$x,lims$upper,lty=2)
lines(lims$x,lims$lower,lty=2)

yt<-cut(y,seq(0,1,0.1))
barplot(as.vector(table(yt)),
          names=as.character(seq(.1,1,.1)),xlab="y-distance",ylab="intensity")

xt<-cut(x,seq(0,1,0.1))
yt<-cut(y,seq(0,1,0.1))

count<-as.vector(table(xt,yt))
table(count)

expected<-100*exp(-1)/factorial(0:5)
expected
d<-numeric(200)
map<-matrix(0,200,200)
ds<-matrix(0,200,200)
for (i in 1:200) {
d<-0
for (k in 1:200) d[k]<-dist(x[i],y[i],x[k],y[k])
map[i,]<-order(d)
ds[i,]<-d
}

attach(sizes)
names(sizes)

plan<-interp(x,y,z)

par(mfrow=c(2,2))
plot(x,y,pch=16)
persp(plan)
contour(plan)
image(plan)

neighbours<-ds[,2]
sizes<-z[map[,2]]

model1<-lm(z~neighbours+sizes)
summary(model1)

model2<-lm(z~neighbours+sizes+x+y)
summary(model2)

model3<-lm(z~neighbours+sizes+x+I(x^2)+y+I(y^2))
summary(model3)

attach(spatialdata)
names(spatialdata)

par(mfrow=c(1,2))
plot(latitude,yield)
plot(longitude,yield)

par(mfrow=c(1,1))
barplot(sort(tapply(yield,variety,mean)))

tapply(yield,Block,mean)

model1<-aov(yield~variety)
summary(model1)

Block<-factor(Block)
model2<-aov(yield~Block+variety+Error(Block))
summary(model2)

model3<-aov(yield~Block+variety+latitude+longitude)
summary(model3)

space<-groupedData(yield~variety|Block)

model4<-gls(yield~variety-1,space)
summary(model4)

plot(Variogram(model4,form=~latitude+longitude))

model5<-update(model4, corr=corSpher(c(28,0.2),form=~latitude+longitude,nugget=T))
summary(model5)

model6<-update(model4, corr=corRatio(c(12.5,0.2),form=~latitude+longitude,nugget=T))

anova(model5,model6)

anova(model4,model6)

plot(Variogram(model6,resType="n"))

plot(model6,resid( ., type="n")~fitted(.),abline=0)

qqnorm(model6,~resid(.,type="n"))

model7<-update(model6,model=yield~variety)
anova(model7)

anova(model6,L=c(-1,0,1))

