# some examples of estimating variograms

 # locations over a 20x20 grid

locs_expand.grid(1:20,1:20)
locs_cbind(locs[,1],locs[,2])
r_seq(0,13,length=50)

 # get a random sample from the 20x20 points
samp1_sample(1:400,20)
samp2_sample(1:400,80)

 # tell Splus to access the spatial module
module(spatial)

par(mfrow=c(1,2),oma=c(25,0,0,0),mar=c(1,3,0,0))
 # make the bm surface with p = 1.5
cbm_covbm(locs,pow=1.5)
z_rmultnorm(1,rep(0,20*20),cbm)
 # now sample from it
df1_data.frame(z=z[1,samp1],x=locs[samp1,1],y=locs[samp1,2])
df2_data.frame(z=z[1,samp2],x=locs[samp2,1],y=locs[samp2,2])
z_matrix(z,ncol=20)
pout_persp(z)
points(perspp(df1$x,df1$y,df1$z,pout),pch=1)
points(perspp(df2$x,df2$y,df2$z,pout))

 # now show the empirical variogram
par(mar=c(6,6,2,2))
plot(variogram(z~x+y,data=df1))
v_variogram(z~x+y,data=df2)
points(v$distance,v$gamma,pch="*")
lines(r,r^(1.5))



par(mfrow=c(1,2),oma=c(25,0,0,0),mar=c(1,3,0,0))
# gaussian surface
cga_covpow(locs,pow=2,scale=4)
z_rmultnorm(1,rep(0,20*20),cga)
 # now sample from it
df1_data.frame(z=z[1,samp1],x=locs[samp1,1],y=locs[samp1,2])
df2_data.frame(z=z[1,samp2],x=locs[samp2,1],y=locs[samp2,2])
z_matrix(z,ncol=20)
pout_persp(z)
points(perspp(df1$x,df1$y,df1$z,pout),pch=1)
points(perspp(df2$x,df2$y,df2$z,pout))
 # now show the empirical variogram
par(mar=c(6,6,2,2))
plot(variogram(z~x+y,data=df1))
v_variogram(z~x+y,data=df2)
points(v$distance,v$gamma,pch="*")
lines(r,(1-exp(-(r/4)^(2))))

par(mfrow=c(1,2),oma=c(25,0,0,0),mar=c(1,3,0,0))
# exponential surface
cga_covpow(locs,pow=1,scale=12)
z_rmultnorm(1,rep(0,20*20),cga)
 # now sample from it
df1_data.frame(z=z[1,samp1],x=locs[samp1,1],y=locs[samp1,2])
df2_data.frame(z=z[1,samp2],x=locs[samp2,1],y=locs[samp2,2])
z_matrix(z,ncol=20)
pout_persp(z)
points(perspp(df1$x,df1$y,df1$z,pout),pch=1)
points(perspp(df2$x,df2$y,df2$z,pout))
 # now show the empirical variogram
par(mar=c(6,6,2,2))
plot(variogram(z~x+y,data=df1))
v_variogram(z~x+y,data=df2)
points(v$distance,v$gamma,pch="*")
lines(r,(1-exp(-(r/12)^(1))))


 # a bunch of realizations - just to see
par(mfrow=c(3,4),oma=c(0,0,3,0),mar=c(4,4,1,1))
for(i in 1:12){
 cga_covpow(locs,pow=2,scale=4)
 z_rmultnorm(1,rep(0,20*20),cga)
 # now sample from it
 df1_data.frame(z=z[1,samp1],x=locs[samp1,1],y=locs[samp1,2])
 df2_data.frame(z=z[1,samp2],x=locs[samp2,1],y=locs[samp2,2])
  # now show the empirical variogram
 plot(variogram(z~x+y,data=df1),xlim=c(0,12),ylim=c(0,1.8))
 v_variogram(z~x+y,data=df2)
 points(v$distance,v$gamma,pch="*")
 lines(r,(1-exp(-(r/4)^(2))))
}
mtext("Realizations from a Gaussian process with Gaussian variogram",
      side=3,line=.5,cex=1,outer=T)

par(mfrow=c(3,4),oma=c(0,0,3,0),mar=c(4,4,1,1))
for(i in 1:12){
 cga_covpow(locs,pow=1,scale=12)
 z_rmultnorm(1,rep(0,20*20),cga)
 # now sample from it
 df1_data.frame(z=z[1,samp1],x=locs[samp1,1],y=locs[samp1,2])
 df2_data.frame(z=z[1,samp2],x=locs[samp2,1],y=locs[samp2,2])
  # now show the empirical variogram
 plot(variogram(z~x+y,data=df1),xlim=c(0,12),ylim=c(0,1.8))
 v_variogram(z~x+y,data=df2)
 points(v$distance,v$gamma,pch="*")
 lines(r,(1-exp(-(r/12)^(1))))
}
mtext("Realizations from a Gaussian process with exponential variogram",
      side=3,line=.5,cex=1,outer=T)

par(mfrow=c(3,4),oma=c(0,0,3,0),mar=c(4,4,1,1))
for(i in 1:12){
 cga_covbm(locs,pow=1.5,scale=10)
 z_rmultnorm(1,rep(0,20*20),cga)
 # now sample from it
 df1_data.frame(z=z[1,samp1],x=locs[samp1,1],y=locs[samp1,2])
 df2_data.frame(z=z[1,samp2],x=locs[samp2,1],y=locs[samp2,2])
  # now show the empirical variogram
 plot(variogram(z~x+y,data=df1),xlim=c(0,12),ylim=c(0,1.8))
 v_variogram(z~x+y,data=df2)
 points(v$distance,v$gamma,pch="*")
 lines(r,(r/10)^(1.5))
}
mtext("Realizations from a Gaussian process with power (p=1.5) variogram",
      side=3,line=.5,cex=1,outer=T)


 # now look at anisotropy plot
anisotropy.plot(z~x+y,data=df2)