# Code to generate unconditional realizations source("dace.sim.r") # # Generate a 2-D realization # set.seed(99) #remove for a "real" realization library(rgl) #for nice plots n=10 x1=x2=seq(0,1,length=n) design=as.matrix(expand.grid(x1,x2)) l1=list(m1=abs(outer(design[,1],design[,1],"-"))) l2=list(m2=abs(outer(design[,2],design[,2],"-"))) l.dez=list(l1=l1,l2=l2) rho=c(0.5,0.5) alpha=2 s2=1 se2=0 z=sim.field(l.dez,rho,s2,se2=se2,alpha=alpha, conditioning=1e-14) persp3d(x1,x2,z,col="grey",xlab="x1",ylab="x2",zlab="GP draw",type='h') points3d(design[,1],design[,2],z) # Example library(geoR) gsim=cbind(design,z) gsim=as.geodata(gsim) vgram=variog(gsim) eyefit(vgram) # Explore different models, etc. # Construct envelope using permutation test. vgram.env=variog.mc.env(gsim, obj.var = vgram) plot(vgram, envelope = vgram.env) # Construct envelope by simulating from the model vgram.ml=likfit(gsim, ini.cov.pars = c(2.5, 1), cov.model="gaussian", fix.nugget = TRUE, nugget=1e-10) vgram.env=variog.model.env(gsim, obj.v = vgram, model.pars = vgram.ml,nsim=100,save.sim=TRUE) plot(vgram, envelope = vgram.env) # Instead of plotting max/min envelope, plot the 5% # and 95% quantiles from the simulations gtmp=gsim nbins=length(vgram$u) vgsims=matrix(0,nrow=100,ncol=nbins) for(i in 1:100) { gtmp$data=vgram.env$simulated.data[,i] vg=variog(gtmp) vgsims[i,]=vg$v plot(vgram, envelope = vgram.env) for(i in 1:nbins) points(vg$u[i],quantile(vgsims[,i],0.05),pch=20) for(i in 1:nbins) points(vg$u[i],quantile(vgsims[,i],0.95),pch=20) # Add in a line at the MLE estimate of the model lines.variomodel(vg$u,cov.model="gaussian", cov.pars=vgram.ml$cov.pars,nugget=1e-10) # Continuity c=1 eps=4 h=seq(0,1,length=100) # Bound from Adler, 1981 bound=c/( abs(log(h))^(1+eps) ) # Gaussian correlation with theta=1 Rh=exp(-h^2) # Exponential correlation with theta=1 Rhe=exp(-abs(h)) plot(h,bound,type='l',lwd=3,col="grey",ylim=c(0,1),xlab="h",ylab="1-R(H)") lines(h,1-Rh,col="blue") lines(h,1-Rhe,col="orange") # Example - single path inference/ergodicity source("dace.sim.r") n=100 design=matrix(seq(0,1,length=n),ncol=1) l1=list(m1=abs(outer(design[,1],design[,1],"-"))) l.dez=list(l1=l1) rho=0.2 alpha=1 s2s=10 rs=1/s2s seed=sample(1:1e5,1) set.seed(seed) s2=1 rho=0.2 # -log(.2)*s2=1.6ish se2=0 z0=sim.field(l.dez,rho,s2,se2=se2,alpha=alpha) z0=z0-mean(z0) plot(design,z0,type='l',xlab="x",ylab="GP draw", lwd=2,ylim=c(-5,5)) set.seed(seed) s2=1*s2s rho=0.2 z1=sim.field(l.dez,rho,s2,se2=se2,alpha=alpha) z1=z1-mean(z1) plot(design,z0,type='l',xlab="x",ylab="GP draw", lwd=2,ylim=c(-5,5)) lines(design,z1,col="red",lwd=2,lty=3) set.seed(seed) s2=1 rho=0.2^rs z2=sim.field(l.dez,rho,s2,se2=se2,alpha=alpha) z2=z2-mean(z2) plot(design,z0,type='l',xlab="x",ylab="GP draw", lwd=2,ylim=c(-5,5)) lines(design,z1,col="red",lwd=2,lty=3) lines(design,z2,col="blue",lwd=2,lty=2) set.seed(seed) s2=1*s2s rho=0.2^rs # -log(rho)*s2=1.6ish z3=sim.field(l.dez,rho,s2,se2=se2,alpha=alpha) z3=z3-mean(z3) plot(design,z0,type='l',xlab="x",ylab="GP draw", lwd=2,ylim=c(-2,2)) lines(design,z1,col="red",lwd=2,lty=3) lines(design,z2,col="blue",lwd=2,lty=2) lines(design,z3,col="green",lwd=2)