Revision 7669be7daa9719001ccbd992cacd4f0dbb4f8f91 authored by Martin Schlather on 19 May 2005, 00:00:00 UTC, committed by Gabor Csardi on 19 May 2005, 00:00:00 UTC
1 parent dcb0ce9
RFtest.Covariance.R
# source("RFtest.Covariance.R")
if (file.exists("source.R")) source("source.R")
epsilon <- 0.000000001
tf <- function(x, l, i=1) {
x <- as.matrix(x)
if (!is.null(l[[i]]$s)) {
return(sqrt(apply(x^2, 1, sum))/l[[i]]$s)
} else {
return(sqrt(apply((x %*% matrix(l[[i]]$a, nrow=ncol(x)))^2, 1, sum)))
}
}
tfst <- function(x, m) {
if (is.matrix(x)) {
if (is.null(m$a)) {
stop("tfst -- no $aniso"); ERR
}
else {
x <- x %*% matrix(m$a, nrow=ncol(x))
return(cbind(sqrt(apply(x[,-ncol(x),drop=FALSE]^2, 1, sum)),
abs(x[,ncol(x)])))
}
} else {
stop("tfst -- iso"); ERR
}
}
st.fct <- function(x, m) {
eval(parse(text=paste(letters[1:6], "<- m$k[", 1:6,"];")))
y <- abs(tfst(x, m))
switch(e,
psi <- sqrt(y[,ncol(y)]^c +1)^d,
psi <- sqrt((d^(-1) * y[,ncol(y)]^c + 1) / (y[,ncol(y)]^c + 1)),
psi <- sqrt(-log(y[,ncol(y)]^c + 1/d)/log(d))
)
if (ncol(y)==2) {
cov <- c("stable", "whittle", "cauchy")[b]
phi <- CovarianceFct(y[,1] / psi, model=cov, param=c(0, m$v, 0, 1, a))
} else {
stopifnot(ncol(y)==1); ERR
phi <- 1
}
return(psi^(-f) * phi)
}
aniso <- function(dim) matrix(runif(dim^2), ncol=dim)
id <- function(dim) diag(dim)
x0 <- list(c(1:4), cbind(1:4, c(3,2,5,6)))
x1 <- list(cbind(0:4, c(0,3,2,5,6)), cbind(0:4, c(0,3,2,5,6), c(0,2:5)))
models <-
list(
list(list(list(model="exp", var=3, scale =5)),
function(x, l) { exp(-tf(x,l,1)) * l[[1]]$v},
x0),
list(list(list(model="exp", var=7, scale =5),"+",
list(model="gencau", var=3, scale =5, k=c(2,3))),
function(x, l) {
exp(-tf(x,l,1)) * l[[1]]$v +
(1+tf(x,l,3)^l[[3]]$k[1])^(-l[[3]]$k[2]/l[[3]]$k[1]) * l[[3]]$v
},
x0),
list(list(list(model="exp", var=7, scale =5),"*",
list(model="gencau", var=3, scale =5, k=c(2,3))),
function(x, l) {
exp(-tf(x,l,1)) * l[[1]]$v *
(1+tf(x,l,3)^l[[3]]$k[1])^(-l[[3]]$k[2]/l[[3]]$k[1]) * l[[3]]$v
},
x0),
list(list(list(model="exp", var=7, aniso=aniso),"*",
list(model="gencau", var=3, aniso=aniso, k=c(2,3))),
function(x, l) {
exp(-tf(x,l,1)) * l[[1]]$v *
(1+tf(x,l,3)^l[[3]]$k[1])^(-l[[3]]$k[2]/l[[3]]$k[1]) * l[[3]]$v
},
x0),
list(list(list(model="nsst", var=7, aniso=aniso,
kappa=c(1.1, 1, 1.5,0.7, 1, 3))),
function(x,l) st.fct(x, l[[1]]),
x1),
list(list(list(model="nsst", var=7, aniso=aniso,
kappa=c(1.1, 2, 1.5,0.7, 2, 3))),
function(x,l) st.fct(x, l[[1]]),
x1),
list(list(list(model="nsst", var=7, aniso=aniso,
kappa=c(1.1, 3, 1.5,0.7, 3, 3))),
function(x,l) st.fct(x, l[[1]]),
x1),
list(list(list(model="nsst", var=7, aniso=aniso,
kappa=c(1.3, 3, 1.2,0.5, 2, 4)),"+",
list(model="nsst", var=5, aniso=aniso,
kappa=c(1.2, 1, 1.5,0.7, 2, 3))),
function(x,l) st.fct(x, l[[1]]) + st.fct(x, l[[3]]),
x1),
list(list(list(model="nsst", var=7, aniso=aniso,
kappa=c(1.3, 3, 1.2,0.5, 2, 4)),"*",
list(model="nsst", var=5, aniso=aniso,
kappa=c(1.2, 1, 1.5,0.7, 2, 3)), "+",
list(model="nsst", var=5, aniso=aniso,
kappa=c(1.2, 1, 1.5,0.7, 1, 3))),
function(x,l)st.fct(x,l[[1]]) * st.fct(x,l[[3]]) + st.fct(x,l[[5]]),
x1),
)
{
for (model in models) {
cat("\n",paste(c(lapply(model[[1]], function(l) l[[1]]), recursive=TRUE)),"..")
for (x in model[[3]]) {
f <- model[[2]]
M <- model[[1]]
for (i in 1:length(model[[1]]))
if (is.function(M[[i]]$aniso))
M[[i]]$aniso <- M[[i]]$aniso(ncol(as.matrix(x)))
#print(M)
if ( any(abs(CovarianceFct(as.matrix(x), M) - f(x, M)) > f(x, M) * epsilon)){
#print(model)
print("X")
print(x)
print("C code")
print(CovarianceFct(x, M))
print("f eval")
print(f(x, M))
print("-----------")
print(abs(CovarianceFct(x, M) - f(x, M)))
print("X-----------")
print(f(x, M) * epsilon)
stop("differences"); ERR
} else cat(". OK")
}
}
cat("\n")
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...