| vario.pgs {RGeostats} | R Documentation |
Calculate the variogram of the underlying GRFs from the categories (e.g lithofacies)
vario.pgs(db, dbprop = NA, rule = rule.input(), props = NA,
model.type = "residuals", flag.rho = FALSE,
dirvect = NA, lag = NA, toldis = 0.5, tolang = NA,
bench=0, cylrad=0, nlag = NA, breaks = NA, calcul = "cov",
opt.code = 0, tolcode = 0, dates=NA,
verbose=FALSE, vario.add=NA)
db |
The |
dbprop |
The |
rule |
The |
props |
Array giving the (constant) proportions of the different facies involved in the Lithotype rule |
model.type |
Character string specifying the model structure to use. One of "full", "symetrical" or "residuals" (default). See details. |
flag.rho |
When TRUE, the correlation coefficient linking the two underlying
Gaussian Random Functions is inferred. Otherwise, the correlation
coefficient as defined in the input |
dirvect |
The directions in which the variogram must be calculated.
For more information, refer to |
lag |
Array containing the distance lags for each calculation direction. If no corresponding lag is not defined, a default lag is calculated so that the maximum distance is equal to half of the field diagonal. For one distance, if the lag is not defined, it is set to the default lag. |
toldis |
Tolerance on the distance, expressed as a proportion of the lag. This parameter is only used in the case of regular lags |
tolang |
Array containing the angular tolerance for each calculation direction. If no angular tolerance is not defined, a default tolerance is calculated as one half of the angle. For one direction, if the angular tolerance is not defined, it is set to the default tolerance. |
bench |
Array containing the bench separation tolerance for each calculation direction. When bench is defined and positive, a pair is discarded if the separation distance between the two points along the third axis is larger than 'bench'. |
cylrad |
Array containing the cylrad tolerance for each calculation direction. When specified and positive, a pair is discarded if the separation distance in the plane orthogonal to the calculation distance is larger than 'cylrad'. |
nlag |
Array containing the number of lags for each calculation direction. If not defined, the defaulted number of lags is set to 10. For one direction, if the number of lags is not defined, it is set to 10. |
breaks |
Vector giving the boundaries of irregular distance lags. For multidrectional computations, breaks is either: a vector (The breaks are then the same for all directions) a list of vectors of the breaks for each direction. If provided, nlag is adapted to the breaks. Intervals are open left and closed right (e.g ]0;15]). |
calcul |
Character string giving the type of structure to compute:
|
opt.code |
Array containing the option concerning the use of the code for each direction. This option defines the strategy for the code usage when considering a pair of samples:
|
tolcode |
Array giving the maximum distance between the codes of the two samples for each direction. |
dates |
When defined and if the 'time' locator is defined, this value corresponds to the bounds of the Date Intervals. It is organized as follows:
|
verbose |
Verbose argument |
vario.add |
When defined, the current variogram computation is concatenated to the contents of the current variogram. |
The model.type argument indicates the structure of the multivariate covariance matrix to fit. The "residuals" option means that the second GRF is supposed to be of the form:
Y_2 = ρ Y_1 + √{1-ρ^2} R
where Y_1 is the first GRF and R is a second GRF independent from Y_1.
The "symmetrical" option indicates that the resulting covariance function is such as C12(h)=C21(h).
The "full" option indicates that there is no constraint.
The "residual" option has to be used if the two GRF's are supposed to be uncorrelated (it is then specified by the field "rho" of the input rule which is equal to 0 and by the input flag.rho which is set to FALSE).
The output is the experimental Covariance of the GRFs.
The application also calculates the scores during this iterative procedure. The stack of these scores can be displayed through a call to the keypair mechanism after the function 'vario.pgs' has been run successfully. The syntax is:
stack <- get.keypair("vario.pgs_stack")
The resulting object 'stack' is organized as a matrix with the following information described per column:
Rank of the Direction
Rank of the Lag
If the GRFs are uncorrelated: the score for each GRF
If the GRFs are correlated: the number of iterations, the score and the gradients
## Not run:
#This part contains several independent examples.
# Example 1: independent bigaussian model with stationary proportions
# Example 2: correlated bigaussian model ("residuals" type)
with estimation of the correlation coefficient
# Example 3: correlated bigaussian model ("symmetrical" type).
The interest is the particular form of the
cross-correlation.
# Example 4: Shifted bigaussian ("full" type).
# Example 5: Non-stationary proportions
#########################################
# Example 1 #
#########################################
#Simulation of 2 independent GRFs and one categorical variable:
set.seed(234)
data(Exdemo_simpgs_std.rule)
db=db.create(x1=runif(1000),x2=runif(1000),ndim=2)
model=model.create(vartype="Gaussian",range=.2,ndim=2)
model2=model.create(vartype="Exponential",range=.3,ndim=2)
db=simtub(dbout=db,model=model,nbtuba=2000)
vario1=vario.calc(db,nlag=30)
db=simtub(dbout=db,model=model2,nbtuba=2000)
vario2=vario.calc(db,nlag=30)
db=db.locate(db,4:5,"z")
db=db.rule(db,rule=Exdemo_simpgs_std.rule)
#Call vario.pgs
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,nlag=30)
#Graphical representation of the Rule
plot(Exdemo_simpgs_std.rule)
title("Composition rule")
dev.off()
#Graphical representation of the data and the variograms
split.screen(c(2,2))
#The data
screen(1)
plot(db[,2],db[,3],col=db[,6],pch=16,xlab="",ylab="",cex=.6)
title("Data: 3 categories")
#Legend for the variograms
screen(2)
plot(c(0,10),c(0,10),cex=0,xaxt="n",yaxt="n",xlab="",ylab="")
legend(0,5,legend=c("Models","Vario from GRFs","Vario from Categories"),
col=c(6,4,1),lty=c(1,1,NA),pch=c(NA,NA,1),box.col=0)
title("Legend for the variograms")
#GRF1
screen(3)
v1=vario.reduce(vario=vario,varcols=1,flag.vario=TRUE)
plot(vario1,col=4,xlab="h",ylab=expression(gamma),reset=FALSE)
points(v1[1]$hh,v1[1]$gg)
plot(model,vario=v1,add=TRUE,col=6)
title("Variograms of the GRF #1")
#GRF2
screen(4)
v2=vario.reduce(vario=vario,varcols=2,flag.vario=TRUE)
plot(vario2,col=4,xlab="h",ylab=expression(gamma),reset=FALSE)
points(v2[1]$hh,v2[1]$gg)
plot(model2,vario=v2,pch=1,add=TRUE,col=6)
title("Variograms of the GRF #2")
rm(db,Exdemo_simpgs_std.rule,model,model2,v1,v2,vario,vario1,vario2)
par(mfrow=c(1,1))
#########################################
# Example 2 #
#########################################
set.seed(23)
data(Exdemo_simpgs_std.rule)
db=db.create(flag.grid=TRUE,nx=5000,ndim=1)
r=-.7
model=model.create(vartype="Gaussian",range=20,sill=matrix(c(1,r,r,r^2),2),
ndim=1)
model=model.create(vartype="Exponential",range=50,sill=matrix(c(0,0,0,1-r^2),2),
model=model,ndim=1)
Exdemo_simpgs_std.rule$rho=r
db=simtub(dbout=db,model=model,nbtuba=2000)
db=db.rule(db,rule=Exdemo_simpgs_std.rule)
#1st call of vario.pgs to determine rho
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,nlag=5,lag=1,flag.rho=T)
Exdemo_simpgs_std.rule$rho=vario$vars[2]
#Second call with fixed and estimated rho
Exdemo_simpgs_std.rule$rho=vario$vars[2]
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,nlag=50,lag=1,flag.rho=F)
plot(model,as.cov=TRUE,vario=vario)
plot(vario,add=TRUE,col=2)
plot(db[,2],db[,3],col=db[,6],pch=16,xlab="",ylab="",cex=.6)
#########################################
# Example 3 #
#########################################
r1=-1
r2=0
s11=.4
s12=.4
range1=500
s21=.6
s22=.6
range2=10000
p1=.8
p2=p3=.09
p4=1-p1-p2-p3
rule=rule.create(c("S","T","F1","F2","T","F3","F4"))
model=model.create(melem.name(5),
sill=matrix(c(s11,sqrt(s11*s12)*r1,sqrt(s11*s12)*r1,s12),2),range=range1,ndim=2)
model=model.create(melem.name(5),
sill=matrix(c(s21,sqrt(s21*s22)*r2,sqrt(s21*s22)*r2,s22),2),range=range2,
model=model,ndim=2)
r=model.eval(model,0,as.cov=T,jvar=2)
rule$rho=r
curve(model.eval(model,x,as.cov=T,ivar=1,jvar=2,),0,6000,ylab="",
ylim=c(-1,1),n=1000)
db=db.create(nx=c(200,200),dx=c(100,100))
db=simtub(dbout=db,model=model,nbtuba=1000)
db=db.rule(db,rule=rule,props=c(p1,p2,p3,p4))
plot(db[,2],db[,3],pch=16,col=db[,6])
vario=vario.pgs(db,rule=rule,lag=1,nlag=50,flag.rho=F,model.type="symmetrical",props=c(p1,p2,p3,p4))
plot(vario)
plot(model,vario=vario,add=T,col=2)
#########################################
# Example 4 #
#########################################
set.seed(2233)
data(Exdemo_simpgs_std.rule)
model=model.create(melem.name(2),range=50,ndim=1)
data(Exdemo_simpgs_std.rule)
Exdemo_simpgs_std.rule$rho=model.eval(model,10,as.cov=T)
db=db.create(nx=5010)
db=simtub(dbout=db,model=model,nbtuba=2000)
db=db.create(nx=5000,z1=db[1:5000,3],z2=db[11:5010,3])
db=db.rule(db,rule=Exdemo_simpgs_std.rule)
Exdemo_simpgs_std.rule$rho=0
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,lag=1,nlag=5,flag.rho=T,model.type="full")
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,lag=1,nlag=50,flag.rho=F,model.type="full")
plot(vario)
x11()
plot(vario[1]$hh[103:201],vario[1]$gg[103:201],xlab="h",ylab=expression(gamma),type="l")
abline(v=0)
curve(model.eval(model,x+10,as.cov=T),add=T,col=2)
#########################################
# Example 5 #
#########################################
set.seed(23)
data(Exdemo_simpgs_std.rule)
db=db.create(nx=5000)
r=-.7
model=model.create(vartype="Gaussian",range=20,
sill=matrix(c(1,r,r,r^2),2))
model=model.create(vartype="Exponential",range=50,
sill=matrix(c(0,0,0,1-r^2),2),model=model)
Exdemo_simpgs_std.rule$rho=r
db=simtub(dbout=db,model=model,nbtuba=2000)
#Simulation of the proportions
model.prop=model.create(vartype="Gaussian",range=1000)
db.prop=db.create(nx=5000)
db.prop=simtub(dbout=db.prop,model=model.prop,seed=12213,nbtuba=2000,nbsimu=3)
proportions=exp(db.prop[,3:4])/(apply(exp(db.prop[,3:4]),1,sum)+1)
proportions=cbind(proportions,1-apply(proportions[,1:2],1,sum))
db=db.add(db,proportions,locnames=c("p1","p2","p3"))
db=db.rule(db,rule=Exdemo_simpgs_std.rule)
par(mfrow=c(1,5))
plot(c(0,1),c(1,5000),cex=0,xlab="Proportions",ylab="",xaxp=c(0,1,2),yaxp=c(0,5000,5),cex.lab=1.2,cex.axis=1.4)
polygon(c(0,db[,5],0),c(1,1:5000,5000),col=1)
polygon(c(db[5000:1,5],1-db[,7]),c(5000:1,1:5000),col=2)
polygon(c(1-db[,7],1,1),c(1:5000,5000,1),col=3)
db.prop=db.prop.thresh(db.prop,db,rule=Exdemo_simpgs_std.rule,model=model)
plot(db[,3],1:5000,type="l",xlab="1st GRF",ylab="",xaxp=c(3,3,2),yaxt="n",cex.lab=1.2,cex.axis=1.4)
lines(db.prop[,6],1:5000,col=2,lwd=1.5)
plot(c(0,1),c(1,5000),cex=0,xlab="Facies",ylab="",xaxt="n",yaxt="n",cex.axis=1.4)
from=c(1,.5+(1:4999)[db[-1,8]!=db[-5000,8]])
to=c(from[-1],5000)
for(i in 1:length(from))
polygon(c(0,1,1,0),c(from[i],from[i],to[i],to[i]),col=db[1,8],db[-1][db[-1,8]!=db[-5000,8],8][i])
vario=vario.pgs(db,rule=Exdemo_simpgs_std.rule,lag=1,nlag=50,flag.rho=F)
plot(vario)
plot(model,vario=vario,add=T,col=2)
## End(Not run)