| Title: | Probabilistic Forecasting using Ensembles and Bayesian Model Averaging |
|---|---|
| Description: | Bayesian Model Averaging to create probabilistic forecasts from ensemble forecasts and weather observations <https://stat.uw.edu/sites/default/files/files/reports/2007/tr516.pdf>. |
| Authors: | Chris Fraley, Adrian E. Raftery, J. McLean Sloughter, Tilmann Gneiting, University of Washington. |
| Maintainer: | Chris Fraley <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 5.1.8 |
| Built: | 2024-12-23 06:21:32 UTC |
| Source: | CRAN |
Computes Brier Scores for climatology, raw ensemble, and ensemble forecasting models given observation thresholds.
brierScore( fit, ensembleData, thresholds, dates = NULL, ...)brierScore( fit, ensembleData, thresholds, dates = NULL, ...)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
thresholds |
One or more threshold values for the Brier score computations. |
dates |
The dates for which the Brier score will be computed.
These dates must be consistent with |
... |
Included for generic function compatibility. |
A data frame giving the Brier Scores for climatology
(empirical distribution of the verifying observations),
ensemble (voting), and ensemble foreacsting models
for the specified thresholds.
A logistic Brier score is also given for the BMAgamma0 model.
G. W. Brier, Verification of forecasts expressed in terms of probability, Monthly Weather Review, 78:1–3, 1950.
T. Gneiting and A. E. Raftery, Strictly proper scoring rules, prediction and estimation, Journal of the American Statistical Association 102:359–378, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) ## End(Not run) hist(prcpTestData$obs) brierScore(prcpTestFit, prcpTestData, thresholds = seq(from = 0, to = .5, by = .1))data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) ## End(Not run) hist(prcpTestData$obs) brierScore(prcpTestFit, prcpTestData, thresholds = seq(from = 0, to = .5, by = .1))
Computes the cumulative distribution function (CDF) of an ensemble forecasting model at observation locations.
cdf( fit, ensembleData, values, dates = NULL, ...)cdf( fit, ensembleData, values, dates = NULL, ...)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
values |
The vector of desired values at which the CDF of the ensemble forecasting model is to be evaluated. |
dates |
The dates for which the CDF will be computed.
These dates must be consistent with |
... |
Included for generic function compatibility. |
This method is generic, and can be applied to any ensemble forecasting
model.
Note the model may have been applied to a power transformation of the data,
but that information is included in the input fit, and
the output is transformed appropriately.
A vector of probabilities corresponding to the CDF at the desired values. Useful for determining propability of freezing, precipitation, etc.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
ensembleBMA,
fitBMA,
quantileForecast
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMAnormal( tempTestData[1:20,], trainingDays = 8) tempTestForc <- quantileForecast( tempTestFit, tempTestData) range(tempTestForc) tempTestCDF <- cdf( tempTestFit, tempTestData, values = seq(from=277, to=282, by = 1)) tempTestCDFdata(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMAnormal( tempTestData[1:20,], trainingDays = 8) tempTestForc <- quantileForecast( tempTestFit, tempTestData) range(tempTestForc) tempTestCDF <- cdf( tempTestFit, tempTestData, values = seq(from=277, to=282, by = 1)) tempTestCDF
Combines BMA models having the same characteristics for different dates.
combine( x, y, ...)combine( x, y, ...)
x |
An |
y |
An |
... |
Other |
Input models are checked for compatibility, and entries from different
inputs having the same dates are eliminated. Dates are ordered
chronologically and ensemble members are ordered in the order in which
they occur in inout x.
An ensembleBMA model that merges the models from each input into
a single model for the common dates.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit12 <- ensembleBMAnormal( tempTestData, trainingDays = 30, dates = c("2008010100","2008010200")) tempTestFit34 <- ensembleBMAnormal( tempTestData, trainingDays = 30, dates = c("2008010300","2008010400")) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit12 <- ensembleBMAnormal( tempTestData, trainingDays = 8, dates = c("2008010100","2008010200")) tempTestFit34 <- ensembleBMAnormal( tempTestData, trainingDays = 8., dates = c("2008010300","2008010400")) tempTestFit <- combine( tempTestFit12, tempTestFit34)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit12 <- ensembleBMAnormal( tempTestData, trainingDays = 30, dates = c("2008010100","2008010200")) tempTestFit34 <- ensembleBMAnormal( tempTestData, trainingDays = 30, dates = c("2008010300","2008010400")) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit12 <- ensembleBMAnormal( tempTestData, trainingDays = 8, dates = c("2008010100","2008010200")) tempTestFit34 <- ensembleBMAnormal( tempTestData, trainingDays = 8., dates = c("2008010300","2008010400")) tempTestFit <- combine( tempTestFit12, tempTestFit34)
Specifies a list of values controling the Bayesian Model Averaging fit of a mixture of gammas to ensemble forecasts for wind speed.
controlBMAgamma(maxIter, tol, power = 1, startupSpeed = NULL, init, optim.control = list(ndeps = rep( sqrt(.Machine$double.eps), 2)))controlBMAgamma(maxIter, tol, power = 1, startupSpeed = NULL, init, optim.control = list(ndeps = rep( sqrt(.Machine$double.eps), 2)))
maxIter |
An integer specifying an upper limit on the number of iterations'
for fitting the BMA mixture via EM. The default is
|
tol |
A numeric convergence tolerance. The EM fit for the mixture of
gammas is terminated when the relative error in successive
objective values in the M-step falls below |
power |
A scalar value giving the power by which the data will be transformed to fit the model for mean of the observations. The default is not to transform the data. The untransformed forecast is used to fit the variance model. |
startupSpeed |
A scalar value giving a global value for the anemometer startup speed,
or the threshold below which a value of 0 is recorded. As this can
vary from station to station and network to network, it may be
preferable to include |
init |
An optional list of initial values for variance coefficients and weights. The default is to start with the variance coefficients equal to 1, and with equal weights for each member of the ensemble. |
optim.control |
Control parameters for the optim function used in the M-step of EM.
The default here is list(ndeps = rep( sqrt(.Machine$double.eps), 2)),
which assigns a smaller finite-difference step size than the
|
A list whose components are the input arguments and their assigned values.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Ensemble Forecasting
using Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTestFit1 <- ensembleBMAgamma(winsTestData, trainingDays = 30, control = controlBMAgamma(maxIter = 100, tol = 1.e-6, startupSpeed =1)) ## End(Not run) # for quick run only; use more training days for forecasting winsTestFit1 <- ensembleBMAgamma(winsTestData[1:14,], trainingDays = 5, control = controlBMAgamma(maxIter = 100, tol = 1.e-6, startupSpeed = 1))data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTestFit1 <- ensembleBMAgamma(winsTestData, trainingDays = 30, control = controlBMAgamma(maxIter = 100, tol = 1.e-6, startupSpeed =1)) ## End(Not run) # for quick run only; use more training days for forecasting winsTestFit1 <- ensembleBMAgamma(winsTestData[1:14,], trainingDays = 5, control = controlBMAgamma(maxIter = 100, tol = 1.e-6, startupSpeed = 1))
Specifies a list of values controling the Bayesian Model Averaging fit of a mixture of gammas with a point mass at 0 to ensemble forecasts for precipitation.
controlBMAgamma0(maxIter = Inf, tol = sqrt(.Machine$double.eps), power = (1/3), rainobs = 10, init = list(varCoefs = NULL, weights = NULL), optim.control = list(ndeps = rep( sqrt(.Machine$double.eps), 2)))controlBMAgamma0(maxIter = Inf, tol = sqrt(.Machine$double.eps), power = (1/3), rainobs = 10, init = list(varCoefs = NULL, weights = NULL), optim.control = list(ndeps = rep( sqrt(.Machine$double.eps), 2)))
maxIter |
An integer specifying an upper limit on the number of iterations
for fitting the BMA mixture via EM. The default is
|
tol |
A numeric convergence tolerance. The EM fit for the mixture of
gammas is terminated when the relative error in successive
objective values in the M-step falls below |
power |
A scalar value giving the power by which the data will be transformed to fit the models for the point mass at 0 and mean of nonzero observations. The default is to use the 1/3 power of the data. The untransformed forecast is used to fit the variance model. |
rainobs |
An integer specifying a minimum number of observations with nonzero precipitation in the training data. When necessary and possible, the training period will be extended backward in increments of days to meet the minimum requirement. It is not possible to fit the BMA model for precipitation without sufficient nonzero observations. The default minimum number is 10. It many instances fewer nonzero observations may suffice, but it could also be that more are needed to model precipitation in some datasets. |
init |
An optional list of initial values for variance coefficients and weights. The default is to start with the variance coefficients equal to 1, and with equal weights for each member of the ensemble. |
optim.control |
Control parameters for the optim function used in the M-step of EM.
The default here is list(ndeps = rep( sqrt(.Machine$double.eps), 2)),
which assigns a smaller finite-difference step size than the
|
A list whose components are the input arguments and their assigned values.
J. M. Sloughter, A. E. Raftery, T Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Ensemble Forecasting
using Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
ensembleBMAgamma0,
fitBMAgamma0
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit1 <- ensembleBMAgamma0( prcpTestData, trainingDays = 30, control = controlBMAgamma0(power = (1/4))) ## End(Not run) # for quick run only; use more training days for forecasting prcpTestFit1 <- ensembleBMAgamma0( prcpTestData[1:14,], trainingDays = 6, control = controlBMAgamma0(power = (1/4)))data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit1 <- ensembleBMAgamma0( prcpTestData, trainingDays = 30, control = controlBMAgamma0(power = (1/4))) ## End(Not run) # for quick run only; use more training days for forecasting prcpTestFit1 <- ensembleBMAgamma0( prcpTestData[1:14,], trainingDays = 6, control = controlBMAgamma0(power = (1/4)))
Specifies a list of values controling the Bayesian Model Averaging fit of a mixture of normals to ensemble forecasts.
controlBMAnormal(maxIter, tol, equalVariance, biasCorrection, init)controlBMAnormal(maxIter, tol, equalVariance, biasCorrection, init)
maxIter |
An integer specifying an upper limit on the number of iterations
for fitting the BMA mixture via EM. The default is |
tol |
A numeric convergence tolerance. The EM fit for the mixture
model is terminated when the relative error in successive
objective values in the M-step falls below |
equalVariance |
A logical value indicating whether or not the variances for the mixture components should to be equal. The default is to constrain them to be equal. |
biasCorrection |
A character string describing the type of bias correction to be used.
|
init |
An optional list of initial values for standard deviations and weights. The default is to start with all standard deviations equal to 1, and with equal weights for each member of the ensemble. |
A list whose components are the input arguments and their assigned values.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
ensembleBMAnormal,
fitBMAnormal
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit1 <- ensembleBMAnormal(tempTestData, trainingDays = 30, control = controlBMAnormal(maxIter = 100, biasCorrection = "additive")) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit1 <- ensembleBMAnormal(tempTestData[1:20,], trainingDays = 5, control = controlBMAnormal(maxIter = 100, biasCorrection = "additive"))data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit1 <- ensembleBMAnormal(tempTestData, trainingDays = 30, control = controlBMAnormal(maxIter = 100, biasCorrection = "additive")) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit1 <- ensembleBMAnormal(tempTestData[1:20,], trainingDays = 5, control = controlBMAnormal(maxIter = 100, biasCorrection = "additive"))
Computes the continuous ranked probability score (CRPS) for univariate ensemble forecasting models.
crps( fit, ensembleData, dates=NULL, nSamples=NULL, seed=NULL, ...) CRPS( fit, ensembleData, dates=NULL, nSamples=NULL, seed=NULL, ...)crps( fit, ensembleData, dates=NULL, nSamples=NULL, seed=NULL, ...) CRPS( fit, ensembleData, dates=NULL, nSamples=NULL, seed=NULL, ...)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
nSamples |
The number of simulation samples for CRPS via simulation.
For the normal model, the default is analytic computation of the CRPS.
For the gamma model with a point mass at 0 (precipitation),
the CRPS is always computed by simulation,
with default |
seed |
Argument to |
dates |
The dates for which the CRPS will be computed.
These dates must be consistent with |
... |
Included for generic function compatibility. |
These methods are generic, and can be applied to all ensemble forecasting
models.
For gamma0 model for precipitation and the gamma model
for wind speed the CRPS is only available through simulation.
The default number of simulation samples is 10,000.
Note that the gamma0 model for precipitation and the
gamma model for wind speed may have been applied to a power
transformation of the data.
For normal models for temperature and pressure, analytic computation
of the CRPS is the default. CRPS will be computed via simulation for
normal models only if nSamples is set to a positive value.
For the bivariate normal model for wind speed and direction, the
CRPS is computed for the marginal wind speed distribution.
crps is a matrix giving the CRPS for each instance in the data
for both the raw ensemble and the median probabilistic forecast. CRPS is a vector giving the mean of the CRPS over all
instances for the raw ensemble and the median probabilistic forecast.
T. Gneiting and A. E. Raftery, Strictly proper scoring rules, prediction and estimation, Journal of the American Statistical Association 102:359–378, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMAnormal( tempTestData[1:20,], trainingDays = 8) crpsValues <- crps( tempTestFit, tempTestData) colMeans(crpsValues) CRPS( tempTestFit, tempTestData)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMAnormal( tempTestData[1:20,], trainingDays = 8) crpsValues <- crps( tempTestFit, tempTestData) colMeans(crpsValues) CRPS( tempTestFit, tempTestData)
Checks that the character form of a vector of dates conforms to YYYYMMDDHH or YYYYMMDD.
dateCheck(YYYYMMDDHH)dateCheck(YYYYMMDDHH)
YYYYMMDDHH |
A character vector (or its factor equivalent) of dates which should be in the form YYYYMMDDHH or YYYYMMDD, in which YYYY specifies the year, MM the month, DD the day, and (optionally) HH the hour. |
If both YYYYMMDDHH and YYYYMMDD are present,
the YYYYMMDD dates are assumed to be in error
even if HH == 00 for all of the longer dates.
Requires the chron library.
A logical vector indicating whether or not each element of YYYYMMDDHH has the correct format.
dateCheck(c("2008043000", "20080431", "20080501"))dateCheck(c("2008043000", "20080431", "20080501"))
This data set gives 48-hour forecasts for 2-m temperature,
precipitation accumulated over the last 24 hours, and maximum wind speed
at SeaTac (KSEA) and Portland (PDX) ariports in 2007/2008 initialized at
00 hours UTC using a 12km grid. The forecasts are based on an 8 member
version of the University of Washington mesoscale ensemble
(Grimit and Mass 2002; Eckel and Mass 2005).
A data frame with 66 rows and 34 columns: idate the initialization date of each forecast/observation,
format YYYYMMDDHH (categorical). vdate the validation date of each forecast/observation,
format YYYYMMDDHH (categorical). latitude the latitude of each forecast/observation (numeric). longitude the longitude of each forecast/observation (numeric). longitude the elevation (in meters) above sea level (numeric). station weather station identifier (categorical). network weather network identifier (categorical).
*.gfs,*.cmcg,*.eta,*.gasp,*.jma,*.ngps,*.tcwb
forecasts from the 8 members of the ensemble (numeric).
*.obs observed values for the weather parameters.
The prefix * is one of T2 for temperature,
PCP24 for precipitation, MAXWSP10 for wind speed.
Temperature is given in Kelvin.
Precipitation amounts are quantized to hundredths of an inch.
Maximum wind speed is defined as the maximum of the hourly
'instantaneous' wind speeds over the previous 18 hours, where an
hourly 'instantaneous' wind speed is a 2-minute average from the
period of two minutes before the hour to on the hour.
The wind speed observations are measured at 10-m above the ground and
discretized when recorded by rounding to the
nearest meter per second.
This is a small dataset provided for the purposes of testing.
Typically forecasting would be performed on much larger datasets.
F. A. Eckel and C. F. Mass, Effective mesoscale, short-range ensemble forecasting, Weather and Forecasting 20:328–350, 2005.
E. P. Grimit and C. F. Mass, Initial results of a mesoscale short-range ensemble forecasting system over the Pacific Northwest, Weather and Forecasting 17:192–205, 2002.
## Not run: # R check data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") #---------------------------------------------------------------------------- obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) MAE( tempFit, tempTestData) CRPS( tempFit, tempTestData) #---------------------------------------------------------------------------- obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) MAE( prcpTestFit, prcpTestData) CRPS( prcpTestFit, prcpTestData) #---------------------------------------------------------------------------- obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") winsTestFit <- ensembleBMAgamma(winsTestData, trainingDays = 30) MAE( winsTestFit, winsTestData) CRPS( winsTestFit, winsTestData) ## End(Not run)## Not run: # R check data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") #---------------------------------------------------------------------------- obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) MAE( tempFit, tempTestData) CRPS( tempFit, tempTestData) #---------------------------------------------------------------------------- obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) MAE( prcpTestFit, prcpTestData) CRPS( prcpTestFit, prcpTestData) #---------------------------------------------------------------------------- obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") winsTestFit <- ensembleBMAgamma(winsTestData, trainingDays = 30) MAE( winsTestFit, winsTestData) CRPS( winsTestFit, winsTestData) ## End(Not run)
Fits a BMA mixture model to ensemble forecasts. Allows specification of a model, training rule, and forecasting dates.
ensembleBMA( ensembleData, trainingDays, dates = NULL, control = NULL, model = NULL, exchangeable = NULL, minCRPS = NULL)ensembleBMA( ensembleData, trainingDays, dates = NULL, control = NULL, model = NULL, exchangeable = NULL, minCRPS = NULL)
ensembleData |
An |
dates |
The dates for which BMA forecasting models are desired.
By default, this will be all dates in |
trainingDays |
An integer giving the number of time steps (e.g. days) in the training period. There is no default. |
control |
A list of control values for the fitting functions.
The default is |
model |
A character string describing the BMA model to be fit.
Current choices are |
exchangeable |
A numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The model fit will have equal weights and parameters within each group.
The default determines exchangeability from |
minCRPS |
A logical variable indicating whether or not to add a postprocessing step after a normal BMA fit to choose the standard deviation so as to minimize the CRPS for the training data. This argument is used only for normal models, and the default is to not do the CRPS minimization for those models because it may require consderably more computation time, expecially when there are many ensemble members. |
If dates are specified in dates that cannot be forecast
with the training rule, the corresponding BMA model parameter outputs will
be missing (NA) but not NULL.
The training rule uses the number of days corresponding to its
length regardless of whether or not the dates are consecutive.
The following methods are available for the output of ensembleBMA:
cdf, quantileForecast, modelParameters,
brierScore, crps, CRPS and MAE.
A list with the following output components:
dateTable |
The table of observations corresponding to the dates in |
trainingDays |
The number of days in the training period as specified in input. |
... |
One or more components corresponding to fitted coefficients for the model. |
weights |
The fitted BMA weights for the mixture components for each ensemble member at each date. |
power |
A scalar value giving the power (if any) by which the data was
transformed for modeling.
The untransformed forecast is used to fit the variance model.
This is input as part of |
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian Model Averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
ensembleData,
ensembleBMAnormal,
ensembleBMAgamma0,
ensembleBMAgamma,
cdf,
quantileForecast,
modelParameters,
brierScore,
crps,
MAE,
controlBMAnormal,
controlBMAgamma0,
controlBMAgamma
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMA( tempTestData, trainingDays = 30, model = "normal") ## equivalent to ## tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMA( tempTestData[1:20,], trainingDays = 8, model = "normal") set.seed(0); exch <- sample(1:length(ens),replace=TRUE) tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], exchangeable = exch, dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMA( tempTestData[1:20,], trainingDays = 8, model = "normal")data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMA( tempTestData, trainingDays = 30, model = "normal") ## equivalent to ## tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMA( tempTestData[1:20,], trainingDays = 8, model = "normal") set.seed(0); exch <- sample(1:length(ens),replace=TRUE) tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], exchangeable = exch, dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMA( tempTestData[1:20,], trainingDays = 8, model = "normal")
Fits a Bayesian Model Averaging mixture of gammas to ensemble forecasts. Intended for predicting wind speed. Allows specification of a training period and forecasting dates.
ensembleBMAgamma( ensembleData, trainingDays, dates = NULL, control = controlBMAgamma(), exchangeable = NULL)ensembleBMAgamma( ensembleData, trainingDays, dates = NULL, control = controlBMAgamma(), exchangeable = NULL)
ensembleData |
An |
trainingDays |
An integer giving the number of time steps (e.g. days) in the training period. There is no default. |
dates |
The dates for which forecasting models are desired.
By default, this will be all dates in |
control |
A list of control values for the fitting functions. The defaults are
given by the function |
exchangeable |
A numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The models fit will have equal weights and parameters within each group.
The default determines exchangeability from |
The output is for all of the dates in ensembleBMA, so there
will be missing entries denoted by NA for dates that are too recent
to be forecast with the training rule.
The following methods are available for ensembleBMAgamma0 objects:
cdf, quantileForecast, modelParameters,
brierScore, crps, CRPS and MAE.
A list with the following output components:
training |
A list containing information on the training length and lag and the number of instances used for training for each modeling day. |
prob0coefs |
The fitted coefficients in the model for the point mass at 0 (probability of zero precipitaion) for each member of the ensemble at each date. |
biasCoefs |
The fitted coefficients in the model for the mean of the gamma components for each member of the ensemble at each date (bias correction). |
varCoefs |
The fitted coefficients for the model for the variance of gamma components for each date. The coefficients are the same for all members of the ensemble. |
weights |
The fitted BMA weights for the gamma components for each ensemble member at each date. |
power |
A scalar value giving to the power by which the data was transformed
to fit the models for the point mass at 0 and the bias model.
The untransformed forecast is used to fit the variance model.
This is input as part of |
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
ensembleData,
controlBMAgamma,
fitBMAgamma,
cdf,
quantileForecast,
modelParameters,
brierScore,
crps,
MAE
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTestFit <- ensembleBMAgamma(winsTestData, trainingDays = 30, control = controlBMAgamma(startupSpeed = 1)) ## equivalent to ## winsTestFit <- ensembleBMA(winsTestData, trainingDays = 30, ## model = "gamma") ## End(Not run) # for quick run only; use more training days for forecasting winsTestFit <- ensembleBMAgamma(winsTestData[1:14,], trainingDays = 5, control = controlBMAgamma(startupSpeed = 1))data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTestFit <- ensembleBMAgamma(winsTestData, trainingDays = 30, control = controlBMAgamma(startupSpeed = 1)) ## equivalent to ## winsTestFit <- ensembleBMA(winsTestData, trainingDays = 30, ## model = "gamma") ## End(Not run) # for quick run only; use more training days for forecasting winsTestFit <- ensembleBMAgamma(winsTestData[1:14,], trainingDays = 5, control = controlBMAgamma(startupSpeed = 1))
Fits a Bayesian Model Averaging mixture of gammas with a point mass at 0 to ensemble forecasts. Intended for predicting precipitation. Allows specification of a training rule and forecasting dates.
ensembleBMAgamma0( ensembleData, trainingDays, dates = NULL, control = controlBMAgamma0(), exchangeable = NULL)ensembleBMAgamma0( ensembleData, trainingDays, dates = NULL, control = controlBMAgamma0(), exchangeable = NULL)
ensembleData |
An |
trainingDays |
An integer giving the number of time steps (e.g. days) in the training period. There is no default. |
dates |
The dates for which forecasting models are desired.
By default, this will be all dates in |
control |
A list of control values for the fitting functions. The defaults are
given by the function |
exchangeable |
A numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The models fit will have equal weights and parameters within each group.
The default determines exchangeability from |
The output is for all of the dates in ensembleBMA, so there
will be missing entries denoted by NA for dates that are too recent
to be forecast with the training rule.
The following methods are available for ensembleBMAgamma0 objects:
cdf, quantileForecast, modelParameters,
brierScore, crps, CRPS and MAE.
A list with the following output components:
training |
A list containing information on the training length and lag and the number of instances used for training for each modeling day. |
prob0coefs |
The fitted coefficients in the model for the point mass at 0 (probability of zero precipitaion) for each member of the ensemble at each date. |
biasCoefs |
The fitted coefficients in the model for the mean of the gamma components for each member of the ensemble at each date (bias correction). |
varCoefs |
The fitted coefficients for the model for the variance of gamma components for each date. The coefficients are the same for all members of the ensemble. |
weights |
The fitted BMA weights for the gamma components for each ensemble member at each date. |
power |
A scalar value giving to the power by which the data was transformed
to fit the models for the point mass at 0 and the bias model.
The untransformed forecast is used to fit the variance model.
This is input as part of |
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
ensembleData,
controlBMAgamma0,
fitBMAgamma0,
cdf,
quantileForecast,
modelParameters,
brierScore,
crps,
MAE
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) ## equivalent to ## prcpTestFit <- ensembleBMA( prcpTestData, trainingDays = 30, ## model = "gamma0") ## End(Not run) # for quick run only; use more training days for forecasting prcpTestFit <- ensembleBMAgamma0( prcpTestData[3:16,], trainingDays = 6)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) ## equivalent to ## prcpTestFit <- ensembleBMA( prcpTestData, trainingDays = 30, ## model = "gamma0") ## End(Not run) # for quick run only; use more training days for forecasting prcpTestFit <- ensembleBMAgamma0( prcpTestData[3:16,], trainingDays = 6)
Fits a Bayesian Model Averaging mixture of normals to ensemble forecasts. Allows specification of a training rule and forecasting dates.
ensembleBMAnormal(ensembleData, trainingDays, dates = NULL, control = controlBMAnormal(), exchangeable = NULL, minCRPS = FALSE)ensembleBMAnormal(ensembleData, trainingDays, dates = NULL, control = controlBMAnormal(), exchangeable = NULL, minCRPS = FALSE)
ensembleData |
An |
trainingDays |
An integer giving the number of time steps (e.g. days) in the training period. There is no default. |
dates |
The dates for which BMA forecasting models are desired.
By default, this will be all dates in |
control |
A list of control values for the fitting functions. The defaults are
given by the function |
exchangeable |
A numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The modeling will have equal weights and parameters within each group.
The default determines exchangeability from |
minCRPS |
A logical variable indicating whether or not to add a postprocessing step after the BMA fit to choose the standard deviation so as to minimize the CRPS for the training data. The default is not to do the CRPS minimization, because it can add considerable extra cost to the computation, particularly when there are many ensemble members. |
The output is for all of the dates in ensembleData, so there
will be missing entries denoted by NA for dates that are too recent
to be forecast with the training rule.
The following methods are available for ensembleBMAnormal objects:
cdf, quantileForecast, modelParameters,
brierScore, crps, CRPS and MAE.
A list with the following output components:
training |
A list containing information on the training length and lag and the number of instances used for training for each modeling day. |
biasCoefs |
The fitted bias-correction coefficients for each ensemble member at each date. |
sd |
The fitted standard deviations for the mixture of normals model at each date. |
weights |
The fitted BMA weights for the normal components for each ensemble member at each date. |
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155-1174, 2005.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
ensembleData,
controlBMAnormal,
fitBMAnormal,
cdf,
quantileForecast,
modelParameters,
brierScore,
crps,
MAE
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## equivalent to ## tempTestFit <- ensembleBMA( tempTestData, trainingDays = 30, ## model = "normal") ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMAnormal( tempTestData[1:20,], trainingDays = 8)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## equivalent to ## tempTestFit <- ensembleBMA( tempTestData, trainingDays = 30, ## model = "normal") ## End(Not run) # for quick run only; use more training days for forecasting tempTestFit <- ensembleBMAnormal( tempTestData[1:20,], trainingDays = 8)
Creates an ensembleData object including ensemble forecasts along
with dates and (optionally) observations. Other descriptive information
such as latitude, longitude, and station type may be included as well.
ensembleData( forecasts, dates = NULL, observations = NULL, ..., forecastHour, initializationTime, startupSpeed = NULL, exchangeable = NULL)ensembleData( forecasts, dates = NULL, observations = NULL, ..., forecastHour, initializationTime, startupSpeed = NULL, exchangeable = NULL)
forecasts |
A matrix or array (for vector quantities) with columns corresponding to forecasts from individual members of an ensemble and rows corresponding to forecasts for the same date. |
dates |
A numeric or character vector or factor specifying the valid dates
for the forecasts. If numeric, it is interpreted as a Julian date
if it has an |
observations |
Optional vector (or matrix for vector quantities) of observed weather conditions corresponding to the forecasts. Must be supplied if the data is to be used for BMA modeling. |
... |
A named list of additional attributes such as latitude, longitude, and startupSpeed for wind speed. |
forecastHour |
A number giving the forecast hour, the time interval between the initialization and forecast times, in units of hours. |
initializationTime |
A number or character string giving the initialization time. |
startupSpeed |
A numeric value specifying a value below which the anemometer readings for wind speed will be recorded as zero. This value is used for all stations when the startup speed is not explicity specified as part of the data. |
exchangeable |
A numeric or character vector or factor indicating groups of ensemble members that are exchangeable (indistinguishable). The models fit will have equal weights and parameters within each group. The same names/labels should be used as for the forecasts. The default assumes that none of the ensemble members are exhangeable. |
For use with batch processing modeling functions (ensembleBMA
etc), instances ensembleData object are assumed
the same forecast hour and initialization time, which should be
specified as part of the object.
Methods for ensembleData objects include ensembleSize,
ensembleForecasts, ensembleValidDates.
Subsetting is possible, but in the case of columns it applies only to
the ensemble forecasts.
For vector wind computations, the velocity should be in the first
column and the direction in the second.
An ensembleData object, incorporating forecasts and
(optionally) observations with the associated valid dates.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
ensembleBMA,
ensembleBMAgamma,
ensembleBMAgamma0,
ensembleBMAnormal
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) ## End(Not run) obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTestFit <- ensembleBMAgamma(winsTestData, trainingDays = 30) ## End(Not run)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTestFit <- ensembleBMAgamma0( prcpTestData, trainingDays = 30) ## End(Not run) obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTestFit <- ensembleBMAgamma(winsTestData, trainingDays = 30) ## End(Not run)
Fits a Bayesian Modeling Averaging mixture model to a given training set.
fitBMA( ensembleData, control = NULL, model = NULL, exchangeable = NULL)fitBMA( ensembleData, control = NULL, model = NULL, exchangeable = NULL)
ensembleData |
An |
control |
A list of control values for the fitting functions.
The default is |
model |
A character string describing the BMA model to be fit.
Current choices are |
exchangeable |
A numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The model fit will have equal weights and parameters
within each group.
The default determines exchangeability from |
This function fits a BMA model to a training data set.
Methods available for fitBMA objects (the output of fitBMA)
include: cdf, quantileForecast, and
modelParameters.
A list with the following output components:
... |
One or more components corresponding to the coeffcients of the model. |
weights |
The fitted BMA weights for the mixture components for each ensemble member. |
nIter |
The number of EM iterations. |
power |
A scalar value giving the power (if any) by which the data was transformed
for modeling.
The untransformed forecast is used to fit the variance model.
This is input as part of |
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
ensembleData,
ensembleBMA,
fitBMAgamma,
fitBMAgamma0,
fitBMAnormal,
cdf,
quantileForecast,
modelParameters,
controlBMAgamma,
controlBMAgamma0,
controlBMAnormal
data(ensBMAtest) ensNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00") tempTrain <- trainingData( tempTestData, trainingDays = 30, date = "2008010100") tempTrainFit <- fitBMA( tempTrain, model = "normal") ## equivalent to ## tempTrainFit <- fitBMAnormal( tempTrain) set.seed(0); exch <- sample(1:length(ens),replace=TRUE) tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], exchangeable = exch, observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00")data(ensBMAtest) ensNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00") tempTrain <- trainingData( tempTestData, trainingDays = 30, date = "2008010100") tempTrainFit <- fitBMA( tempTrain, model = "normal") ## equivalent to ## tempTrainFit <- fitBMAnormal( tempTrain) set.seed(0); exch <- sample(1:length(ens),replace=TRUE) tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], exchangeable = exch, observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00")
Fits a Bayesian Modeling Averaging mixture of gammas. Intended for wind speed forecasts.
fitBMAgamma( ensembleData, control = controlBMAgamma(), exchangeable = NULL)fitBMAgamma( ensembleData, control = controlBMAgamma(), exchangeable = NULL)
ensembleData |
An |
control |
A list of control values for the fitting functions. The defaults are
given by the function |
exchangeable |
An optional numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The model fit will have equal weights and parameters within each group.
If supplied, this argument will override any specification of
exchangeability in |
This function fits a BMA model to a training data set.
It is called by ensembleBMAgamma, which can produce a sequence
of fits over a larger precipitation data set.
Methods available for the output of fitBMA include:
cdf, quantileForecast, and
modelParameters.
A list with the following output components:
biasCoefs |
The fitted coefficients in the model for the mean of nonzero observations for each member of the ensemble (used for bias correction). |
varCoefs |
The fitted coefficients for the model for the variance of nonzero observations (these are the same for all members of the ensemble). |
weights |
The fitted BMA weights for the gamma components for each ensemble member. |
nIter |
The number of EM iterations. |
power |
A scalar value giving to the power by which the data was transformed
to fit the models for the point mass at 0 and the bias model.
The untransformed forecast is used to fit the variance model.
This is input as part of |
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
ensembleData,
controlBMAgamma,
ensembleBMAgamma,
cdf,
quantileForecast,
modelParameters
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], startupSpeed = 1, forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTrain <- trainingData( winsTestData, trainingDays = 30, date = "2008010100") ## End(Not run) # for quick run only; use more training days for forecasting winsTrain <- trainingData( winsTestData, trainingDays = 10, date = "2008010100") winsTrainFit <- fitBMAgamma( winsTrain) ## equivalent to ## winsTrainFit <- fitBMA( winsTrain, model = "gamma")data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("MAXWSP10","obs", sep = ".") ens <- paste("MAXWSP10", ensMemNames, sep = ".") winsTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], startupSpeed = 1, forecastHour = 48, initializationTime = "00") ## Not run: # R check winsTrain <- trainingData( winsTestData, trainingDays = 30, date = "2008010100") ## End(Not run) # for quick run only; use more training days for forecasting winsTrain <- trainingData( winsTestData, trainingDays = 10, date = "2008010100") winsTrainFit <- fitBMAgamma( winsTrain) ## equivalent to ## winsTrainFit <- fitBMA( winsTrain, model = "gamma")
Fits a Bayesian Modeling Averaging mixture of gammas with a point mass at 0 to a given training set. Intended for precipitation forecasts.
fitBMAgamma0( ensembleData, control = controlBMAgamma0(), exchangeable = NULL)fitBMAgamma0( ensembleData, control = controlBMAgamma0(), exchangeable = NULL)
ensembleData |
An |
control |
A list of control values for the fitting functions. The defaults are
given by the function |
exchangeable |
An optional numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The model fit will have equal weights and parameters within each group.
If supplied, this argument will override any specification of
exchangeability in |
This function fits a BMA model to a training data set.
It is called by ensembleBMAgamma0, which can produce a sequence
of fits over a larger precipitation data set.
Methods available for the output of fitBMA include:
cdf, quantileForecast, and
modelParameters.
A list with the following output components:
prob0coefs |
The fitted coefficients in the model for the point mass at 0 (probability of zero precipitation) for each member of the ensemble. |
biasCoefs |
The fitted coefficients in the model for the mean of nonzero observations for each member of the ensemble (used for bias correction). |
varCoefs |
The fitted coefficients for the model for the variance of nonzero observations (these are the same for all members of the ensemble). |
weights |
The fitted BMA weights for the gamma components for each ensemble member. |
nIter |
The number of EM iterations. |
power |
A scalar value giving to the power by which the data was transformed
to fit the models for the point mass at 0 and the bias model.
The untransformed forecast is used to fit the variance model.
This is input as part of |
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
ensembleData,
controlBMAgamma0,
ensembleBMAgamma0,
cdf,
quantileForecast,
modelParameters
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTrain <- trainingData( prcpTestData, trainingDays = 30, date = "2008010100") ## End(Not run) # quick run only; use more training days for forecasting prcpTrain <- trainingData( prcpTestData, trainingDays = 10, date = "2008010100") prcpTrainFit <- fitBMAgamma0( prcpTrain) ## equivalent to ## prcpTrainFit <- fitBMA( prcpTrain, model = "gamma0")data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("PCP24","obs", sep = ".") ens <- paste("PCP24", ensMemNames, sep = ".") prcpTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check prcpTrain <- trainingData( prcpTestData, trainingDays = 30, date = "2008010100") ## End(Not run) # quick run only; use more training days for forecasting prcpTrain <- trainingData( prcpTestData, trainingDays = 10, date = "2008010100") prcpTrainFit <- fitBMAgamma0( prcpTrain) ## equivalent to ## prcpTrainFit <- fitBMA( prcpTrain, model = "gamma0")
Fits a Bayesian Model Averaging mixture of normals to a given training set.
fitBMAnormal( ensembleData, control = controlBMAnormal(), exchangeable = NULL)fitBMAnormal( ensembleData, control = controlBMAnormal(), exchangeable = NULL)
ensembleData |
An |
control |
A list of control values for the fitting functions. The defaults are
given by the function |
exchangeable |
An optional numeric or character vector or factor indicating groups of
ensemble members that are exchangeable (indistinguishable).
The models have equal weights and parameters within each group.
If supplied, this argument will override any specification of
exchangeability in |
This function fits a BMA model to a training data set.
It is called by ensembleBMAnormal, which can produce a sequence
of fits over a larger data set.
Methods available for the output of fitBMAnormal include:
cdf, quantileForecast, and modelParameters.
A list with the following output components:
biasCoefs |
The fitted bias-correction coefficients. |
sd |
The fitted standard deviations for the mixture of normals model
(equal or varying across components according to the |
weights |
The fitted BMA weights for the normal components for each ensemble member. |
nIter |
The number of EM iterations. |
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian Model Averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
ensembleData,
controlBMAnormal,
ensembleBMAnormal,
cdf,
quantileForecast,
modelParameters
data(ensBMAtest) ensNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00") tempTrain <- trainingData( tempTestData, trainingDays = 30, date = "2008010100") tempTrainFit <- fitBMAnormal( tempTrain)data(ensBMAtest) ensNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00") tempTrain <- trainingData( tempTestData, trainingDays = 30, date = "2008010100") tempTrainFit <- fitBMAnormal( tempTrain)
Converts Julian dates to YYYYMMDDHH or YYYYMMDD character format.
julTOymdh( julianDates, origin = NULL, dropHour = NULL)julTOymdh( julianDates, origin = NULL, dropHour = NULL)
julianDates |
A numeric vector specifying Julian dates. |
origin |
A named vector specifying the month, day, and year for the
origin of the Julian dates. The default is
|
dropHour |
A logical value indicating whether of not the hour information
should be drop from the specifiation of the dates if none of the
Julian dates are fractional.
The default is |
Requires the chron library.
A character vector or numeric equivalent of dates in the form YYYYMMDDHH or YYYYMMDD, in which YYYY specifies the year, MM the month, DD the day, and (optionally) HH the hour corresponding to the Julian input.
data(ensBMAtest) julianIdates <- ymdhTOjul(ensBMAtest$idate) all.equal( julTOymdh(julianIdates), as.character(ensBMAtest$idate)) all.equal( ymdhTOjul(ensBMAtest$vdate), julianIdates+2)data(ensBMAtest) julianIdates <- ymdhTOjul(ensBMAtest$idate) all.equal( julTOymdh(julianIdates), as.character(ensBMAtest$idate)) all.equal( ymdhTOjul(ensBMAtest$vdate), julianIdates+2)
Computes the mean absolute error (MAE) for ensemble forecasting models.
MAE( fit, ensembleData, dates=NULL, ...)MAE( fit, ensembleData, dates=NULL, ...)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
dates |
The dates for which the CRPS and MAE will be computed.
These dates must be consistent with |
... |
Included for generic function compatibility. |
This method is generic, and can be applied to all ensemble forecasting
models.
Note the model may have been applied to a power transformation of the data,
but that information is included in the input fit, and
the output is transformed appropriately.
A vector giving the MAE for the deterministic forecasts associated with the raw ensemble and for the ensemble forecasting model. This is the mean absolute difference of the raw ensemble medians and the observations, and the mean absolute difference of the median forecast and the observations (as in Sloughter et al. 2007). \ Note that Raftery et al. 2005 uses the mean absolute difference of the raw ensemble means and the observations, and the mean absolute difference of the BMA predictive mean and the observations.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised in 2010).
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) MAE( tempTestFit, tempTestData)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) MAE( tempTestFit, tempTestData)
Extracts model parameters for ensemble forecasting models.
modelParameters( fit, ...)modelParameters( fit, ...)
fit |
A model fit to ensemble forecasting data. |
... |
For ensemble fits involving dates, there is an additional |
A list of parameters (including weights) corresponding to the ensemble forecasting model for the specified dates. The list may also include a power by which the forecasts were transformed to obtain the model parameters.
ensembleBMAgamma,
ensembleBMAgamma0,
ensembleBMAnormal,
fitBMAgamma,
fitBMAgamma0,
fitBMAnormal
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) modelParameters( tempTestFit, date = "2008010100") tempTrain <- trainingData( tempTestData, date = "2008010100", trainingDays = tempTestFit$training$days) tempTrainFit <- fitBMAnormal( tempTrain) modelParameters( tempTrainFit)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) modelParameters( tempTestFit, date = "2008010100") tempTrain <- trainingData( tempTestData, date = "2008010100", trainingDays = tempTestFit$training$days) tempTrainFit <- fitBMAnormal( tempTrain) modelParameters( tempTrainFit)
Computes the probabilty integral transform (PIT) of a BMA ensemble forecasting model at observation locations.
pit( fit, ensembleData, dates = NULL, randomizeATzero=FALSE, ...)pit( fit, ensembleData, dates = NULL, randomizeATzero=FALSE, ...)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
dates |
The dates for which the CDF will be computed.
These dates must be consistent with |
randomizeATzero |
For the |
... |
Included for generic function compatibility. |
Most often used for computing PIT histograms to assess calibration of
forecasts, in which case the observations in ensembleData would
be those used in modeling fit.
Instances in ensembleData without verifying observations
are ignored.
Note the model may have been applied to a power transformation of the data,
but that information is included in the input fit, and
the output is transformed appropriately.
The PIT is a continuous analog of the verification rank.
The value of the BMA cumulative distribution function CDF
corresponding to the fit at the observed values in ensembleData.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
T. Gneiting, F. Balabdaoui and A. Raftery, Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society, Series B 69:243–268, 2007.
J. M. Sloughter, A. E. Raftery, T Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
pitHist,
verifRankHist,
ensembleBMA,
fitBMA,
quantileForecast
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) tempTestForc <- quantileForecast( tempTestFit, tempTestData) range(tempTestForc) tempTestPIT <- pit( tempTestFit, tempTestData)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) tempTestForc <- quantileForecast( tempTestFit, tempTestData) range(tempTestForc) tempTestPIT <- pit( tempTestFit, tempTestData)
Computes the probability integral transform of the obervations relative to the BMA forecast, and plots its histogram.
pitHist( fit, ensembleData, dates=NULL)pitHist( fit, ensembleData, dates=NULL)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
dates |
The dates for which the CDF will be computed.
These dates must be consistent with |
PIT histograms are used to assess calibration of
forecasts, in which case the observations in ensembleData would
be those used in modeling fit.
Instances in ensembleData without verifying observations
are ignored.
In the case of the gamma0 model for precipitation,
observations of zero precipitation are randomized within their
probabilistics range to avoid a false
impression of bias.
Note the model may have been applied to a power transformation of the data,
but that information is included in the input fit, and
the output is transformed appropriately.
The PIT is a continuous analog of the verification rank.
The value of the BMA cumulative distribution function CDF
corresponding to the fit at the observed values in ensembleData.
The corresponding histogram is also plotted.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155-1174, 2005.
T. Gneiting, F. Balabdaoui and A. Raftery, Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society, Series B 69:243–268, 2007.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
ensembleData,
pit,
verifRankHist.
data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") ## Not run: # this takes time # the PIT should be evaluated over relatively long periods srftFITall <- ensembleBMA( srftData, model = "normal", trainingDays = 25) srftPIT <- pitHist( srftFITall, srftData) ## End(Not run)data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") ## Not run: # this takes time # the PIT should be evaluated over relatively long periods srftFITall <- ensembleBMA( srftData, model = "normal", trainingDays = 25) srftPIT <- pitHist( srftFITall, srftData) ## End(Not run)
Plots the Predictive Distribution Function (PDF) of an ensemble forecasting model.
## S3 method for class 'ensembleBMAgamma' plot( x, ensembleData, dates=NULL, ask=TRUE, ...) ## S3 method for class 'ensembleBMAgamma0' plot( x, ensembleData, dates=NULL, ask=TRUE, ...) ## S3 method for class 'ensembleBMAnormal' plot( x, ensembleData, dates=NULL, ask=TRUE, ...) ## S3 method for class 'fitBMAgamma' plot( x, ensembleData, dates=NULL, ...) ## S3 method for class 'fitBMAgamma0' plot( x, ensembleData, dates=NULL, ...) ## S3 method for class 'fitBMAnormal' plot( x, ensembleData, dates=NULL, ...)## S3 method for class 'ensembleBMAgamma' plot( x, ensembleData, dates=NULL, ask=TRUE, ...) ## S3 method for class 'ensembleBMAgamma0' plot( x, ensembleData, dates=NULL, ask=TRUE, ...) ## S3 method for class 'ensembleBMAnormal' plot( x, ensembleData, dates=NULL, ask=TRUE, ...) ## S3 method for class 'fitBMAgamma' plot( x, ensembleData, dates=NULL, ...) ## S3 method for class 'fitBMAgamma0' plot( x, ensembleData, dates=NULL, ...) ## S3 method for class 'fitBMAnormal' plot( x, ensembleData, dates=NULL, ...)
x |
A model fit to ensemble forecasting data. |
ensembleData |
An |
dates |
The dates for which the PDF will be computed.
These dates must be consistent with |
ask |
A logical value indicating whether or not the user should be prompted for the next plot. |
... |
Included for generic function compatibility. |
This method is generic, and can be applied to any ensemble forecasting
model.
The colored curves are the weighted PDFs of the ensemble members,
and the bold curve is the overall PDF. The vertical black line represents
the median forecast, and the dotted back lines represent the .1 and .9
quartiles. The vertical orange line is the verifying observation (if
any).
Exchangeable members are represented in the plots by the weighted
group sum rather than by the indivdual weighted PDFs of each member.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190-202, 2010.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) plot(tempTestFit, tempTestData) ## End(Not run)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) plot(tempTestFit, tempTestData) ## End(Not run)
Produces contour, image, or perspective plot of a forecast using loess prediction on a grid.
plotProbcast( forecast, longitude, latitude, nGrid = 65, type = c("image", "contour", "persp"), ..., interpolate = FALSE, span = 0.75, maps = NULL)plotProbcast( forecast, longitude, latitude, nGrid = 65, type = c("image", "contour", "persp"), ..., interpolate = FALSE, span = 0.75, maps = NULL)
forecast |
Numeric vector of forecasts. |
longitude |
Numeric vector giving the longitude of each forecast location. |
latitude |
Numeric vector giving the latitude of each forecast location. |
nGrid |
Number of grid points for |
type |
A character string indicating the desired plot type.
Should be one of either |
... |
Additional arguments to be passed to the plotting method. |
interpolate |
A logical variable indicating whether or not a |
span |
Smoothing parameter for |
maps |
A logical value indicating whether or not to include
a map outline. The default is to include an outline
if |
If the fields library is loaded, a legend (and optionally
a map outline) will be included in image plots.
An image, contour, or perspective plot of the forecast.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") ## Not run: # R check bmaFit <- ensembleBMA( srftData, date = "2004012900", trainingDays = 25, model = "normal") bmaForc <- quantileForecast( bmaFit, srftData, date = "2004012900", quantiles = c(.1, .5, .9)) obs <- srftData$date == "2004012900" lat <- srftData$latitude[obs] lon <- srftData$longitude[obs] plotProbcast( bmaForc[,"0.5"], lat, lon, type = "contour", interpolate = TRUE) title("Median Forecast") plotProbcast( srftData$obs[obs], lat, lon, type = "contour", interpolate = TRUE) title("Observed Surface Temperature") data(srftGrid) memberLabels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftGridData <- ensembleData(forecasts = srftGrid[,memberLabels], latitude = srftGrid[,"latitude"], longitude = srftGrid[,"longitude"], forecastHour = 48, initializationTime = "00") gridForc <- quantileForecast( bmaFit, srftGridData, date = "2004021400", quantiles = c( .1, .5, .9)) library(fields) plotProbcast(gridForc[,"0.5"],lon=srftGridData$lon, lat=srftGridData$lat,type="image",col=rev(rainbow(100,start=0,end=0.85))) title("Median Grid Forecast for Surface Temperature", cex = 0.5) probFreeze <- cdf( bmaFit, srftGridData, date = "2004021400", value = 273.15) plotProbcast(probFreeze, lon=srftGridData$lon, lat=srftGridData$lat, type="image",col=gray((32:0)/32)) title("Probability of Freezing", cex = 0.5) ## End(Not run)data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") ## Not run: # R check bmaFit <- ensembleBMA( srftData, date = "2004012900", trainingDays = 25, model = "normal") bmaForc <- quantileForecast( bmaFit, srftData, date = "2004012900", quantiles = c(.1, .5, .9)) obs <- srftData$date == "2004012900" lat <- srftData$latitude[obs] lon <- srftData$longitude[obs] plotProbcast( bmaForc[,"0.5"], lat, lon, type = "contour", interpolate = TRUE) title("Median Forecast") plotProbcast( srftData$obs[obs], lat, lon, type = "contour", interpolate = TRUE) title("Observed Surface Temperature") data(srftGrid) memberLabels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftGridData <- ensembleData(forecasts = srftGrid[,memberLabels], latitude = srftGrid[,"latitude"], longitude = srftGrid[,"longitude"], forecastHour = 48, initializationTime = "00") gridForc <- quantileForecast( bmaFit, srftGridData, date = "2004021400", quantiles = c( .1, .5, .9)) library(fields) plotProbcast(gridForc[,"0.5"],lon=srftGridData$lon, lat=srftGridData$lat,type="image",col=rev(rainbow(100,start=0,end=0.85))) title("Median Grid Forecast for Surface Temperature", cex = 0.5) probFreeze <- cdf( bmaFit, srftGridData, date = "2004021400", value = 273.15) plotProbcast(probFreeze, lon=srftGridData$lon, lat=srftGridData$lat, type="image",col=gray((32:0)/32)) title("Probability of Freezing", cex = 0.5) ## End(Not run)
A subset of daily 48 hour forecasts of 24 hour accumulated
precipitation over the US Pacific Northwest region
from December 2002 to January 2003 based
on a 9 member version of the University of Washington mesoscale ensemble
(Grimit and Mass 2002; Eckel and Mass 2005).
Precipitation amounts are quantized to hundredths of an inch.
Note that forecasts are not available for some of the interim dates.
A data frame with 175 rows and 15 columns: CENT,AVN,CMCG,ETA,GASP,JMA,NGAPS,TCWB,UKMO
forecasts from the 9 members of the ensemble (numeric).observation the observed accumulated precipitation (numeric).date the date of each forecast/observation,
format YYYYMMDDHH (categorical). station weather station identifier (categorical). latitude the latitude of each weather station (numeric).longitude the longitude of each weather station (numeric). elevation the elevation of each weather station (numeric).
This dataset is a small subset of the data used in Sloughter et al. (2006), provided for the purposes of testing. Typically forecasting would be performed on much larger datasets.
E. P. Grimit and C. F. Mass, Initial results of a mesoscale short-range ensemble forecasting system over the Pacific Northwest, Weather and Forecasting 17:192–205, 2002.
F. A. Eckel and C. F. Mass, Effective mesoscale, short-range ensemble forecasting, Weather and Forecasting 20:328–350, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3309–3320, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
## Not run: # R check data(prcpDJdata) data(prcpFit) prcpForc <- quantileForecast( prcpFit, prcpDJdata, date = "20030113", quantiles = c( .1, .5, .9)) ## End(Not run)## Not run: # R check data(prcpDJdata) data(prcpFit) prcpForc <- quantileForecast( prcpFit, prcpDJdata, date = "20030113", quantiles = c( .1, .5, .9)) ## End(Not run)
The ensembleBMAgamma0 model fit with a 30 day training period to the
precipitation data set from
http://www.stat.washington.edu/MURI,
which gives daily daily 48 hour forecasts of 24 hour accumulated
precipitation over the US Pacific Northwest region from December 12, 2002
through March 31, 2005 on a 9 member version of the University of Washington
mesoscale
ensemble (Grimit and Mass 2002; Eckel and Mass 2005).
Precipitation amounts are quantized to hundredths of an inch.
A list with the following arguments:
dateTableA named vector in which the names are the dates and the entries are the number of observations for each date.
trainingRuleThe training rule used to compute the model fits.
prob0coefsThe coefficients in the logistic regression for probability of zero precipitation.
biasCoefsThe coefficients in the linear regression for bias correction.
varCoefsThe variance coefficients of the models.
weightsThe BMA weights for the models.
powerAn scalar value giving the power by which the forecasts are transformed for the BMA fitting.
E. P. Grimit and C. F. Mass, Initial results of a mesoscale short-range ensemble forecasting system over the Pacific Northwest, Weather and Forecasting 17:192–205, 2002.
F. A. Eckel and C. F. Mass, Effective mesoscale, short-range ensemble forecasting, Weather and Forecasting 20:328–350, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3309–3320, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
## Not run: # R check data(prcpFit) modelParameters(prcpFit, date = "20030113") data(prcpGrid) prcpGridData <- ensembleData(forecasts = prcpGrid[,1:9], latitude = prcpGrid[,"latitude"], longitude = prcpGrid[,"longitude"], forecsatHour = 48, initializationTime = "00") # probability of precipitation 1 - cdf( prcpFit, prcpGridData, value = 0) # probability of precipitation above 0.25 in 1 - cdf( prcpFit, prcpGridData, date = "20030115", value = 25) ## End(Not run)## Not run: # R check data(prcpFit) modelParameters(prcpFit, date = "20030113") data(prcpGrid) prcpGridData <- ensembleData(forecasts = prcpGrid[,1:9], latitude = prcpGrid[,"latitude"], longitude = prcpGrid[,"longitude"], forecsatHour = 48, initializationTime = "00") # probability of precipitation 1 - cdf( prcpFit, prcpGridData, value = 0) # probability of precipitation above 0.25 in 1 - cdf( prcpFit, prcpGridData, date = "20030115", value = 25) ## End(Not run)
This data set gives 48-hour forecasts of 24 hour accumulated precipitation on a grid of locations in the US Pacific Northwest initialized on January 13, 2003 OOZ and valid on January 15, 2003 OOZ. The ensemble forecasts come from a nine member version of the University of Washington Mesoscale Ensemble (Grimit and Mass 2002; Eckel and Mass 2005). Precipitation amounts are quantized to hundredths of an inch.
A data frame with 8188 rows and 11 columns: avn/gfs,cent,cmcg,eta,gasp,jma,ngps,tcwb,ukmo
forecasts from the 9 members of the ensemble (numeric). latitude the latitude of each forecast (numeric). longitude the longitude of each forecast (numeric).
E. P. Grimit and C. F. Mass, Initial results of a mesoscale short-range ensemble forecasting system over the Pacific Northwest, Weather and Forecasting 17:192–205, 2002.
F. A. Eckel and C. F. Mass, Effective mesoscale, short-range ensemble forecasting, Weather and Forecasting 20:328–350, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Ensemble Forecasting
using Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2009.
## Not run: # R check data(prcpGrid) prcpGridData <- ensembleData(forecasts = prcpGrid[,1:9], latitude = prcpGrid[,"latitude"], longitude = prcpGrid[,"longitude"], forecastHour = 48, initilaizationTime = "00") data(prcpFit) # median forecast for Jan 15, 2003 at the grid points quantileForecast( prcpFit, prcpGridData, date = "20030115") ## End(Not run)## Not run: # R check data(prcpGrid) prcpGridData <- ensembleData(forecasts = prcpGrid[,1:9], latitude = prcpGrid[,"latitude"], longitude = prcpGrid[,"longitude"], forecastHour = 48, initilaizationTime = "00") data(prcpFit) # median forecast for Jan 15, 2003 at the grid points quantileForecast( prcpFit, prcpGridData, date = "20030115") ## End(Not run)
Computes quantiles for the probability distribution function (PDF) for ensemble forecasting models.
quantileForecast( fit, ensembleData, quantiles = 0.5, dates=NULL, ...)quantileForecast( fit, ensembleData, quantiles = 0.5, dates=NULL, ...)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
quantiles |
The vector of desired quantiles for the PDF of the BMA mixture model. |
dates |
The dates for which the quantile forecasts will be computed.
These dates must be consistent with |
... |
Included for generic function compatibility. |
This method is generic, and can be applied to any ensemble forecasting
model.
Note the model may have been applied to a power transformation of the data,
but that information is included in the input fit, and
the output is transformed appropriately.
This can be used to compute prediction intervals for the PDF.
For the bivariate normal model for wind speed and direction, the
CRPS is computed for the marginal wind speed distribution.
A vector of forecasts corresponding to the desired quantiles.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155–1174, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) tempTestForc <- quantileForecast( tempTestFit, tempTestData) ## Not run: # R check data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") srftFit <- ensembleBMAnormal(srftData, date = "2004013100", trainingDays = 25) data(srftGrid) srftGridData <- ensembleData(forecasts = srftGrid[ ,labels], latitude = srftGrid$lat, longitude = srftGrid$lon, forecastHour = 48, initializationTime = "00") srftGridForc <- quantileForecast( srftFit, srftGridData, date = "2004013100") ## End(Not run)data(ensBMAtest) ensMemNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensMemNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], dates = ensBMAtest[,"vdate"], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], forecastHour = 48, initializationTime = "00") ## Not run: # R check tempTestFit <- ensembleBMAnormal( tempTestData, trainingDays = 30) ## End(Not run) tempTestForc <- quantileForecast( tempTestFit, tempTestData) ## Not run: # R check data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") srftFit <- ensembleBMAnormal(srftData, date = "2004013100", trainingDays = 25) data(srftGrid) srftGridData <- ensembleData(forecasts = srftGrid[ ,labels], latitude = srftGrid$lat, longitude = srftGrid$lon, forecastHour = 48, initializationTime = "00") srftGridForc <- quantileForecast( srftFit, srftGridData, date = "2004013100") ## End(Not run)
This data set gives 48-hour forecasts of 2-m surface temperature and the
associated observations for the US Pacific Northwest from January 1, 2004
to February 28, 2004. The ensemble forecasts come from an eight-member
version of the University of Washington Mesoscale Ensemble
(Grimit and Mass 2002; Eckel and Mass 2005).
Temperatures are measured in kelvins.
Note that forecasts are not available for some of the interim dates.
A data frame with 36826 rows and 15 columns: CMCG,ETA,GASP,GFS,JMA,NGAPS,TCWB,UKMO
forecasts from the 8 members of the ensemble (numeric). observation the observed surface temperature (numeric). date the date of each forecast/observation set,
in the format YYYYMMDDHH (categorical). latitude the latitude of each forecast (numeric). longitude the longitude of each forecast (numeric). station weather station identifier (categorical). type weather station type (categorical).
F. A. Eckel and C. F. Mass, Effective mesoscale, short-range ensemble forecasting, Weather and Forecasting 20:328–350, 2005.
E. P. Grimit and C. F. Mass, Initial results of a mesoscale short-range ensemble forecasting system over the Pacific Northwest, Weather and Forecasting 17:192–205, 2002.
V. J. Berrocal, A. E. Raftery and T. Gneiting, Combining spatial and ensemble information in probabilistic weather forecasts, Monthly Weather Review 133:1386–1402, 2007.
V. J. Berrocal, A. E. Raftery, T. Gneiting and R. C. Steed, Probabilistic Weather Forecasting for Winter Road Maintenance, Journal of the American Statistical Association, 2010 (to appear).
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
## Not run: # R check data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") srftFit <- ensembleBMAnormal( srftData, date = "2004013100", trainingDays = 25) ## End(Not run)## Not run: # R check data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") srftFit <- ensembleBMAnormal( srftData, date = "2004013100", trainingDays = 25) ## End(Not run)
This data set gives 48-hour forecasts of 2-m surface temperature
on a grid of locations in the US Pacific Northwest
initialized on January 29, 2004 00UTC and valid on January 31, 2004 00UTC.
The ensemble forecasts come from an eight member version of the
University of Washington Mesoscale Ensemble
(Grimit and Mass 2002; Eckel and Mass 2005).
Temperatures are measured in kelvins.
Note that forecasts are not available for some of the interim dates.
A data frame with 10098 rows and 10 columns: CMCG,ETA,GASP,GFS,JMA,NGAPS,TCWB,UKMO
forecasts from the 8 members of the ensemble (numeric). latitude the latitude of each forecast (numeric). longitude the longitude of each forecast (numeric).
F. A. Eckel and C. F. Mass, Effective mesoscale, short-range ensemble forecasting, Weather and Forecasting 20:328–350, 2005.
E. P. Grimit and C. F. Mass, Initial results of a mesoscale short-range ensemble forecasting system over the Pacific Northwest, Weather and Forecasting 17:192–205, 2002.
V. J. Berrocal, A. E. Raftery and T. Gneiting, Combining spatial and ensemble information in probabilistic weather forecasts, Monthly Weather Review 133:1386–1402, 2007.
V. J. Berrocal, A. E. Raftery, T. Gneiting and R. C. Steed, Probabilistic Weather Forecasting for Winter Road Maintenance, Journal of the American Statistical Association, 2010 (to appear).
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
## Not run: # R check data(srft) data(srftGrid) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") srftFit <- ensembleBMAnormal( srftData, date = "2004013100", trainingDays = 25) srftGridData <- ensembleData( forecasts = srftGrid[ ,labels], latitude = srftGrid$lat, longitude = srftGrid$lon, forecastHour = 48, initializationTime = "00") CRPS( srtGridData, srftFit) ## End(Not run)## Not run: # R check data(srft) data(srftGrid) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") srftFit <- ensembleBMAnormal( srftData, date = "2004013100", trainingDays = 25) srftGridData <- ensembleData( forecasts = srftGrid[ ,labels], latitude = srftGrid$lat, longitude = srftGrid$lon, forecastHour = 48, initializationTime = "00") CRPS( srtGridData, srftFit) ## End(Not run)
Extracts a subset of an ensembleData object corresponding
to a given date and number of training days.
trainingData( ensembleData, trainingDays, date)trainingData( ensembleData, trainingDays, date)
ensembleData |
An |
trainingDays |
An integer specifying the number of days in the training period. |
date |
The date for which the training data is desired. |
The most recent days are used for training regardless of whether or not they are consecutive.
An ensembleData object corresponding to the training data for
the given date relative to ensembleData.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155-1174, 2005.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3309–3320, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensembles and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
data(ensBMAtest) ensNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00") tempTrain <- trainingData( tempTestData, trainingDays = 30, date = "2008010100") tempTrainFit <- fitBMAnormal( tempTrain)data(ensBMAtest) ensNames <- c("gfs","cmcg","eta","gasp","jma","ngps","tcwb","ukmo") obs <- paste("T2","obs", sep = ".") ens <- paste("T2", ensNames, sep = ".") tempTestData <- ensembleData( forecasts = ensBMAtest[,ens], observations = ensBMAtest[,obs], station = ensBMAtest[,"station"], dates = ensBMAtest[,"vdate"], forecastHour = 48, initializationTime = "00") tempTrain <- trainingData( tempTestData, trainingDays = 30, date = "2008010100") tempTrainFit <- fitBMAnormal( tempTrain)
Computes the median, 10th and 90th percentile forecasts, and plots the corresponding observations.
verifPlot( fit, ensembleData, dates = NULL)verifPlot( fit, ensembleData, dates = NULL)
fit |
A model fit to ensemble forecasting data. |
ensembleData |
An |
dates |
The dates for which the CDF will be computed.
These dates must be consistent with |
A matrix giving the median, 10th and 90th percentile forecasts for the ensemble data at the specified dates. If observations are available, they are plotted along with the forecasts in order of increasing 90th percentile forecast.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155-1174, 2005.
T. Gneiting, F. Balabdaoui and A. Raftery, Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society, Series B 69:243–268, 2007.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
data(prcpFit) data(prcpDJdata) forc <- verifPlot( prcpFit, prcpDJdata, date = "20030113")data(prcpFit) data(prcpDJdata) forc <- verifPlot( prcpFit, prcpDJdata, date = "20030113")
Computes the rank of verifying observations relative to the corresponding ensemble forecasts and plots its histogram.
verifRankHist( forecasts, observations)verifRankHist( forecasts, observations)
forecasts |
A matrix of ensemble forecasts, in which the rows corresponds to locations and times and the columns correspond to the individual ensemble members. |
observations |
A vector of observations corresponding to the locations and times of the forecasts. |
The verification rank is used to assess calibration of a forecast ensemble. A more uniform verification rank histogram indicates better calibartion.
A vector giving the rank of verifying observations relative to the corresponding ensemble forecasts. The verification rank historgram is plotted.
A. E. Raftery, T. Gneiting, F. Balabdaoui and M. Polakowski, Using Bayesian model averaging to calibrate forecast ensembles, Monthly Weather Review 133:1155-1174, 2005.
T. Gneiting, F. Balabdaoui and A. Raftery, Probabilistic forecasts, calibration and sharpness. Journal of the Royal Statistical Society, Series B 69:243–268, 2007.
J. M. Sloughter, A. E. Raftery, T. Gneiting and C. Fraley, Probabilistic quantitative precipitation forecasting using Bayesian model averaging, Monthly Weather Review 135:3209–3220, 2007.
C. Fraley, A. E. Raftery, T. Gneiting and J. M. Sloughter,
ensembleBMA: An R Package for Probabilistic Forecasting
using Ensemble and Bayesian Model Averaging,
Technical Report No. 516R, Department of Statistics, University of
Washington, 2007 (revised 2010).
C. Fraley, A. E. Raftery, T. Gneiting, Calibrating Multi-Model Forecast Ensembles with Exchangeable and Missing Members using Bayesian Model Averaging, Monthly Weather Review 138:190–202, 2010.
J. M. Sloughter, T. Gneiting and A. E. Raftery, Probabilistic wind speed forecasting using ensembles and Bayesian model averaging, Journal of the American Statistical Association, 105:25–35, 2010.
data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") use <- ensembleValidDates(srftData) >= "2004013000" verifRankHist( ensembleForecasts(srftData[use,]), dataVerifObs(srftData[use,]))data(srft) labels <- c("CMCG","ETA","GASP","GFS","JMA","NGPS","TCWB","UKMO") srftData <- ensembleData( forecasts = srft[ ,labels], dates = srft$date, observations = srft$obs, latitude = srft$lat, longitude = srft$lon, forecastHour = 48, initializationTime = "00") use <- ensembleValidDates(srftData) >= "2004013000" verifRankHist( ensembleForecasts(srftData[use,]), dataVerifObs(srftData[use,]))
Converts YYYYMMDDHH or YYYYMMDD dates to Julian dates.
ymdhTOjul( YYYYMMDDHH, origin = c(month = 1, day = 1, year = 2000))ymdhTOjul( YYYYMMDDHH, origin = c(month = 1, day = 1, year = 2000))
YYYYMMDDHH |
A character vector (or its factor equivalent) of dates in the form YYYYMMDDHH or YYYYMMDD, in which YYYY specifies the year, MM the month, DD the day, and (optionally) HH the hour. |
origin |
A named vector specifying the month, day, and year for the
origin of the Julian dates. The default is
|
Requires the chron library.
A vector of Julian dates corresponding to YYYYMMDDHH.
The vector has "origin" and "dropHour"
attributes which give the origin for the Julian output and
indicate whether or not the original format included the hour.
data(ensBMAtest) julianVdates <- ymdhTOjul(ensBMAtest$vdate) all.equal( julTOymdh(julianVdates), as.character(ensBMAtest$vdate)) all.equal( ymdhTOjul(ensBMAtest$idate), julianVdates-2)data(ensBMAtest) julianVdates <- ymdhTOjul(ensBMAtest$vdate) all.equal( julTOymdh(julianVdates), as.character(ensBMAtest$vdate)) all.equal( ymdhTOjul(ensBMAtest$idate), julianVdates-2)