https://github.com/cran/ensembleBMA
Tip revision: f80c44bfcea4b355579f751ff30ba3c239b90084 authored by Chris Fraley on 23 January 2009, 00:00:00 UTC
version 4.1.1
version 4.1.1
Tip revision: f80c44b
cdf.fitBMAgamma.R
`cdf.fitBMAgamma` <-
function(fit, ensembleData, values, dates=NULL, ...)
{
powfun <- function(x,power) x^power
powinv <- function(x,power) x^(1/power)
weps <- 1.e-4
if (!is.null(dates)) warning("dates ignored")
M <- matchEnsembleMembers(fit,ensembleData)
nForecasts <- ensembleSize(ensembleData)
if (!all(M == 1:nForecasts)) ensembleData <- ensembleData[,M]
# remove instances missing all forecasts
M <- apply(ensembleForecasts(ensembleData), 1, function(z) all(is.na(z)))
ensembleData <- ensembleData[!M,]
nObs <- nrow(ensembleData)
CDF <- matrix( NA, nrow = nObs, ncol = length(values))
dimnames(CDF) <- list(ensembleObsLabels(ensembleData), as.character(values))
nForecasts <- ensembleSize(ensembleData)
ensembleData <- ensembleForecasts(ensembleData)
WEIGHTS <- fit$weights
if (!all(Wmiss <- is.na(WEIGHTS))) {
for (i in 1:nObs) {
f <- ensembleData[i,]
M <- is.na(f) | Wmiss
VAR <- (fit$varCoefs[1] + fit$varCoefs[2]*f)^2
fTrans <- sapply(f, powfun, power = fit$power)
MEAN <- apply(rbind(1, fTrans) * fit$biasCoefs, 2, sum)
W <- WEIGHTS
if (any(M)) {
W <- W + weps
W <- W[!M]/sum(W[!M])
}
CDF[i,] <- sapply( sapply( values, powfun, power = fit$power),
cdfBMAgamma, WEIGHTS = W, MEAN = MEAN[!M], VAR = VAR[!M])
}
}
CDF
}