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)