model <-"gencauchy"
param <- c(0, 1, 0, 1, 1, 2)
estparam <- c(0, NA, 0, NA, NA, 2) ## NA means: "to be estimated"
## sequence in `estparam' is
## mean, variance, nugget, scale, (+ further model parameters)
## So, mean, variance, and scale will be estimated here.
## Nugget is fixed and equals zero.
points <- 100
x <- runif(points,0,3)
y <- runif(points,0,3) ## 300 random points in square [0, 3]^2
## simulate data according to the model:
d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=1000) #1000
## fit the data:
Print(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
lower=c(0.1, 0.1, 0.1), upper=c(1.9, 5, 2)))
#########################################################
## The next two estimations give about the same result.
## For the first the sill is fixed to 1.5. For the second the sill
## is reached if the estimated variance is smaller than 1.5
estparam <- c(0, NA, NA, NA, NA, NA)
Print(v <- fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
sill=1, use.nat=FALSE)) ## gencauchy works better with use.nat=FALSE
estmodel <- list("+",
list("$", var=NA, scale=NA,
list("gencauchy", alpha=NA, beta=NA)
),
list("$", var=NA, list("nugget"))
)
parampositions(model=estmodel, dim=2)
f <- function(variab) c(variab, max(0, 1.0 - variab[1]))
Print(v2 <- fitvario(x=cbind(x,y), data=d, model=estmodel,
lower = c(TRUE, TRUE, TRUE, TRUE, FALSE),
transform=f, use.nat=FALSE))
#########################################################
## estimation of coupled parameters (alpha = beta, here)
# source("RandomFields/tests/source.R")
f <- function(param) param[c(1:3,3,4)]
Print(fitvario(x=cbind(x,y), data=d, model=estmodel,
lower=c(TRUE, TRUE, TRUE, FALSE, TRUE),
transform=f))
#########################################################
## estimation in a anisotropic framework
x <- y <- (1:6)/4
model <- list("$", aniso=matrix(nc=2, c(4,2,-2,1)), var=1.5,
list("exp"))
z <- GaussRF(x=x, y=y, grid=TRUE, model=model, n=10)
estmodel <-list("$", aniso=matrix(nc=2, c(NA,NA,-2,1)), var=NA,
list("exp"))
Print(fitvario(as.matrix(expand.grid(x, y)), data=z,
model=estmodel, nphi=20))
#########################################################
## estimation with trend (formula)
model <- list("$", var=1, scale=2, list("gauss"))
estmodel <- list("$", var=NA, scale=NA, list("gauss"))
x <- seq(-pi,pi,pi/2)
n <- 5
data <- GaussRF(x, x, gridtri=FALSE, model=model,
trend=function(X1,X2) sin(X1) + 2*cos(X2),n=n)
Print(v <- fitvario(x, x, data=data, gridtrip=FALSE,
model=estmodel,
trend=~sin(X1)+cos(X1)+sin(X2)+cos(X2)))
########################################################
## estimation of anisotropy matrix with two identical ##
## diagonal elements ##
x <- c(0, 5, 0.4)
model <- list("$", var=1, scale=1, list("exponential"))
z <- GaussRF(x, x, x, model=model, gridtriple=TRUE, n=10, Print=2)
est.model <- list("+",
list("$", var=NA, aniso=diag(c(NA, NA, NA)), list("exponen")),
list("$", var=NA, list("nugget")))
parampositions(est.model, dim=3)
trafo <- function(variab) {variab[c(1:2, 2:4)]}
lower <- c(TRUE, TRUE, FALSE, TRUE, TRUE) # which parameter to be estimated
fitlog <- fitvario(x, x, x, gridtriple=TRUE, data=z, model=est.model,
transform=trafo, lower=lower)
str(fitlog$ml)Run the code above in your browser using DataLab