#
############################################################
## ##
## T. Gneiting, W. Kleiber, M. Schlather, JASA 2010 ##
## ##
############################################################
## Here the analysis in the above mentioned paper is replicated.
## The results differ slightly from those in the paper, as the algorithm
## has been improved, and the estimation has been nearly fully automated.
## In particular, user supplied lower and upper bound for estimating
## the independent model are no longer required.
## As a result, the corresponding MLE fit deteriorates, whereas
## the other fits improve slightly.
#################################
## get the data ##
#################################
library(fields)
miles2km <- 1.608
data(weather)
## P and T are so different in scale so that they have
## to be normalized first:
sdP <- sd(weather[, 1])
sdT <- sd(weather[, 2])
PT <- cbind( weather[, 1] / sdP, weather[, 2] / sdT )
Dist.mat <- rdist.earth(weather[,3:4])
Dist.mat <- Dist.mat[lower.tri(Dist.mat)] ## in miles
#################################
## first: marginal estimation ##
#################################
uni.est <-
list("+",
list("$", var=NA, list("nugget")),
list("$", var=NA, scale=NA, list("whittle", nu=NA))
)
fvP <- fitvario(Distances=Dist.mat, truedim=2, data=PT[, 1],
model=uni.est, grid=FALSE, ml="ml", Print=1)$ml # -153.9
fvT <- fitvario(Distances=Dist.mat, truedim=2, data=PT[, 2],
model=uni.est, grid=FALSE, ml="ml", Print=1)$ml # -138.45
#
#################################
## second: parsimoninous model ##
#################################
puni2bi <- function(mP, mT, lower) {
nugg1 <- mP[[2]]$var
nugg2 <- mT[[2]]$var
nu1 <- mP[[3]][[4]]$nu
nu2 <- mT[[3]][[4]]$nu
s <- mean(c(mP[[3]]$scale, mT[[3]]$scale))
c1 <- mP[[3]]$var
c2 <- mT[[3]]$var
if (missing(lower)) {
rhored <- 0
f <- 1
} else if (lower) {
rhored <- -1
f <- 0.2
nugg1 <- nugg2 <- 0
} else {
rhored <- 1
f <- 4
nugg1 <- c1 / 2
nugg2 <- c2 / 2
}
return(list("+",
list("M", M=matrix(nc=2, c(sqrt(nugg1), 0, 0, sqrt(nugg2))),
list("nugget")),
list("parsbiWM",
nu = c(nu1 * f, nu2 * f),
s = s * f,
c = c(c1 * f, c2 * f), rhored=rhored
)
) )
}
p.est.model <-
list("+",
list("M", M=matrix(nc=2, c(NA, 0, 0, NA)),
list("nugget")),
list("parsbiWM",
nu = c(NA, NA),
s = NA,
c = c(NA, NA), rhored=NA
))
## takes a rather long time (15 min in 2011)
pars <- fitvario(Distances=Dist.mat, truedim = 2, data=PT,
model=p.est.model, grid=FALSE, trend=list(0,0),
lower= puni2bi(fvP$model, fvT$model, lower=TRUE),
upper= puni2bi(fvP$model, fvT$model, lower=FALSE),
users.guess=puni2bi(fvP$model, fvT$model),
Print=2, ml="ml")$ml ## -280.9791
##
#################################
## third: full biwm model ##
#################################
pars2full <- function(model, lower){
s <- model[[3]]$s
nuP <- model[[3]]$nu[1]
nuT <- model[[3]]$nu[2]
tauP <- model[[2]]$M[1]
tauT <- model[[2]]$M[4]
cP <- model[[3]]$c[1]
cT <- model[[3]]$c[2]
Rhored <- model[[3]]$rhored
nured <- 1
if (missing(lower)) {
f <- 1
} else if (lower) {
nured <- 0.1
f <- 0.5
Rhored <- -1
tauP <- tauT <- 0
} else {
f <- 2
tauP <- max(f * tauP, cP / 10)
tauT <- max(f * tauT, cT / 10)
Rhored <- 0 ## as estimated Rhored is negativ
}
list("+",
list("M", M=matrix(nc=2, c(tauP, 0, 0, tauT)),
list("nugget")),
list("biWM",
nu = c(nuP * f, nuT * f), nured=nured,
s = c(s * f, s * f), s12=s * f,
c = c(cP * f, cT * f), rhored=Rhored
) )
}
est.model <-
list("+",
list("M", M=matrix(nc=2, c(NA, 0, 0, NA)),
list("nugget")),
list("biWM",
nu = c(NA, NA), nured=NA,
s = c(NA, NA), s12=NA,
c = c(NA, NA), rhored=NA
))
fv <- fitvario(Distances=Dist.mat, truedim = 2, data=PT, Print=2,
model=est.model, grid=FALSE, trend=list(0,0),
lower=pars2full(pars$model, lower=TRUE),
upper=pars2full(pars$model, lower=FALSE),
users.guess=pars2full(pars$model), ml="ml")$ml # -280.6910
###
#################################
## output ##
#################################
Factor <- function(nu1, nu2, nu12, a1, a2, a12, d) {
gamma(nu1 + d/2) * gamma(nu2 + d/2) / gamma(nu1) / gamma(nu2) *
(gamma(nu12) / gamma(nu12+d/2))^2 * (a1^nu1 * a2^nu2 / a12^(2*nu12) )^2
}
InfQ <- function(nu1, nu2, nu12, a1, a2, a12, d) {
gamma <- (2*nu12+d) * a1^2 * a2^2 - (nu2 +d/2)*a1^2*a12^2 -
(nu1 +d/2) *a2^2*a12^2
beta <- (2 * nu12 - nu2 + d/2) * a1^2 + (2 * nu12 - nu1 + d/2) * a2^2 -
(nu1 + nu2 + d) *a12^2
alpha <- 2 * nu12 - nu1 -nu2
rednu12 <- 0.5 * (nu1 + nu2) / nu12
if (rednu12 == 1.0) {
t1sq <- t2sq <- max(0, if (beta==0) 0 else -gamma / beta)
inf <- 1
} else {
inf <- Inf
discr <- beta^2 - 4 * alpha * gamma
t1sq <- t2sq <- 0
if (discr >= 0) {
discr <- sqrt(discr)
t1sq = max(0, (-beta + discr) / (2.0 * alpha))
t2sq = max(0, (-beta - discr) / (2.0 * alpha))
}
}
tsq <- c(0, t1sq, t2sq)
return(min(inf, (a12^2 + tsq)^(2.0 * nu12 + d) /
(a1^2 + tsq)^(nu1 + d/2) / (a2^2 + tsq)^(nu2 + d/2) ))
}
uni2full <- function(mP, mT) {
nuggP <- mP[[2]]$var
nuggT <- mT[[2]]$var
nuP <- mP[[3]][[4]]$nu
nuT <- mT[[3]][[4]]$nu
sP <- mP[[3]]$scale
sT <- mT[[3]]$scale
c1 <- mP[[3]]$var
c2 <- mT[[3]]$var
return(list("+",
list("M", M=matrix(nc=2, c(sqrt(nuggP), 0, 0, nuggT)),
list("nugget")),
list("biWM",
nu = c(nuP, nuT), nured=1,
s = c(sP, sT), s12=1,
c = c(c1, c2), rhored=0)
) )
}
PaperOutput <- function(m, sdP, sdT) {
m[[2]]$M <- m[[2]]$M * c(sdP, 0, 0, sdT)
m[[3]]$c <- m[[3]]$c * c(sdP, sdT)^2
sP <- m[[3]]$s[1] * miles2km
sT <- m[[3]]$s[2] * miles2km
sPT <- m[[3]]$s12 * miles2km
m[[3]]$s <- c(sP, sT)
m[[3]]$s12 <- sPT
d <- 2
ml <- fitvario(Distances=Dist.mat * miles2km,
truedim = d, data=t(t(PT) * c(sdP, sdT)),
m=m, grid=FALSE, trend=list(0,0), ml="ml")$ml$ml
nuP <- m[[3]]$nu[1]
nuT <- m[[3]]$nu[2]
nuPT <- 0.5 * (nuP + nuT) / m[[3]]$nured
f <- Factor(nuP, nuT, nuPT, 1/sP, 1/sT, 1/sPT, d)
infQ <- InfQ(nuP, nuT, nuPT, 1/sP, 1/sT, 1/sPT, d)
return(list(
model = m,
sigmaP = sqrt(m[[3]]$c[1]),
sigmaT = sqrt(m[[3]]$c[2]),
nuP = nuP,
nuT = nuT,
nuPT = nuPT,
inv.aP = sP,
inv.aT = sT,
inv.aPT= sPT,
rhoPT = m[[3]]$rhored * sqrt(f * infQ),
tauP = m[[2]]$M[1],
tauT = m[[2]]$M[4],
ml = ml
))
}
Print(PaperOutput(uni2full(fvP$model, fvT$model), sdP, sdT)) # -1277.11
Print(PaperOutput(pars2full(pars$model), sdP, sdT)) # -1265.73
Print(PaperOutput(fv$model, sdP, sdT)) # -1265.30Run the code above in your browser using DataLab