https://github.com/cran/RandomFields
Tip revision: 73d7f2b428b288bbf57445dd5e6f51fd3a84847e authored by Martin Schlather on 30 June 2009, 00:00:00 UTC
version 1.3.40
version 1.3.40
Tip revision: 73d7f2b
test.mle.R
# source("test.mle.R")
if (EXTENDED.TESTING <- file.exists("source.R")) {
print("A")
source("source.R")
print("B")
source("mle.old.R")
}
print("C")
if (EXTENDED.TESTING) {
try(load("rs"))
options(warn=1) ## 4
## 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.
n <- 10
points <- 100
points <- 15
##
model <-"gencauchy"
x <- runif(points,0,3)
y <- runif(points,0,3) ## 100 random points in square [0, 3]^2
param <- c(0, 1, 0, 1, 1, 2)
estparam <- c(0, NA, 0, NA, NA, 2) ## NA means: "to be estimated"
estparam <- param ## NA means: "to be estimated"
estparam[runif(length(estparam)) < 0.5] <- NA
print(estparam)
save(file="rs", .Random.seed)
d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=n)
if (FALSE)
repeat{
save(file="rs", .Random.seed)
d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=n)
str(fitvario(Print=3,
x=cbind(x,y), data=d, model=model, param=estparam,
lower=c(0.1, 0.1), upper=c(1.9, 5),
lsq.methods=c("self", "plain", "sqrt.nr", "sd.inv"),
#cross="cross.ign"
)[c("variogram", "values")], vec.len=6, digits.d=5)
str(try(Ofitvario(x=cbind(x,y), data=d, model=model, param=estparam,
lower=c(0.1, 0.1), upper=c(1.9, 5))[c("mle", "mle.value")]),
vec.len=6, digits.d=5)
print(estparam)
}
## 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)
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam, sill=1.5))
str(Ofitvario(x=cbind(x,y), data=d, model=model, param=estparam, sill=1.5))
estparam <- c(0, NA, NaN, NA, NA, NA)
parampositions(model=model, param=estparam)
f <- function(param) {
param[5] <- max(0, 1.5 - param[1])
return(param)
}
str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam,
sill=1, transform=f))
str(Ofitvario(x=cbind(x,y), data=d, model=model, param=estparam,
sill=1, transform=f))
## the next call gives a warning, since the user may programme
## strange things in this setup, and the program cannot check it.
estparam <- c(0, NA, NA, NA, NA, NaN)
parampositions(model=model, param=estparam)
f <- function(param) {param[3] <- param[2]; param}
unix.time(str(fitvario(x=cbind(x,y), data=d, model=model,
param=estparam, transform=f, standard.style=TRUE), vec.len=6))
unix.time(str(Ofitvario(x=cbind(x,y), data=d, model=model,
param=estparam, transform=f, standard.style=TRUE), vec.len=6))
## much better programmed, but also much slower:
estmodel <- list(list(model="gencauchy", var=NA, scale=NA,
kappa=list(NA, function(m) m[[1]]$kappa[1])))
unix.time(str(fitvario(x=cbind(x,y), data=d, model=estmodel),
vec.len=6))
unix.time(str(Ofitvario(x=cbind(x,y), data=d, model=estmodel),
vec.len=6))
} # source exists