| weather {RandomFields} | R Documentation |
Meteorological dataset, which consists of difference between forecasts and observations (forcasts minus observations) of temperature and pressure at 157 locations in the North American Pacific Northwest.
data(weather)
This data frame contains the following columns:
in units of Pascal
in units of degree Celcius
longitudinal coordinates of the locations
latitude coordinates of the locations
The forecasts are from the GFS member of the University of Washington regional numerical weather prediction ensemble (UWME; Grimit and Mass 2002; Eckel and Mass 2005); they were valid on December 18, 2003 at 4 pm local time, at a forecast horizon of 48 hours.
The data were obtained from Cliff Mass and Jeff Baars in the University of Washington Department of Atmospheric Sciences.
Eckel, A. F. and Mass, C. F. (2005) Aspects of effective mesoscale, short-range ensemble forecasting Wea. Forecasting 20, 328-350.
Gneiting, T., Kleiber, W. and Schlather, M. (2011) Matern cross-covariance functions for multivariate random fields J. Amer. Statist. Assoc. 105, 1167-1177.
Grimit, E. P. and Mass, C. F. (2002) Initial results of a mesoscale short-range forecasting system over the Pacific Northwest Wea. Forecasting 17, 192-205.
## Not run:
#
############################################################
## ##
## T. Gneiting, W. Kleiber, M. Schlather, JASA 2011 ##
## ##
############################################################
## 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.30
## End(Not run)