Nothing
fitBMAnormal <-
function (ensembleData, control = controlBMAnormal(), exchangeable = NULL)
{
ZERO <- 1e-100
if (is.null(exchangeable))
exchangeable <- ensembleGroups(ensembleData)
if (length(unique(exchangeable)) == length(exchangeable))
exchangeable <- NULL
if (!(nullEX <- is.null(exchangeable))) {
namEX <- as.character(exchangeable)
uniqueEX <- unique(namEX)
nEX <- length(uniqueEX)
splitEX <- split(seq(along = exchangeable), exchangeable)
matEX <- sapply(uniqueEX, function(i, x) {
x <- as.numeric(x == i)
x/sum(x)
}, x = exchangeable)
dimnames(matEX) <- list(NULL, uniqueEX)
}
nObs <- dataNobs(ensembleData)
if (!nObs)
stop("no observations")
ensembleData <- ensembleData[!dataNA(ensembleData, dates = FALSE),
]
nObs <- dataNobs(ensembleData)
if (!nObs)
stop("no observations")
obs <- dataVerifObs(ensembleData)
ensMemNames <- ensembleMembers(ensembleData)
nForecasts <- length(ensMemNames)
if (is.null(sd <- control$init$sd))
sd <- sd(obs)
weights <- if (is.null(control$init$weights) || any(is.na(control$init$weights)))
1
else control$init$weights
if (length(weights) == 1)
weights <- rep(weights, nForecasts)
weights <- weights/sum(weights)
weights <- pmax(weights, 1e-04)
weights <- weights/sum(weights)
if (!is.null(names(weights)))
weights <- weights[ensMemNames]
MISSmat <- is.na(ensembleForecasts(ensembleData))
MEAN <- RSQ <- matrix(NA, nObs, nForecasts)
dimnames(MEAN) <- dimnames(RSQ) <- NULL
switch(control$biasCorrection, regression = {
if (nullEX) {
meanFit <- apply(ensembleForecasts(ensembleData),
2, function(x, y) {
components <- c("coefficients", "fitted.values",
"residuals", "model")
lm(y ~ x, na.action = na.omit)[components]
}, y = obs)
biasCoefs <- lapply(meanFit, function(x) x$coefficients)
biasCoefs <- as.matrix(data.frame(biasCoefs))
dimnames(biasCoefs) <- list(NULL, ensMemNames)
FIT <- lapply(meanFit, function(x) x$fitted.values)
RES <- lapply(meanFit, function(x) x$residuals)
for (i in 1:nForecasts) {
I <- !MISSmat[, i]
MEAN[I, i] <- FIT[[i]]
RSQ[I, i] <- RES[[i]]
}
MEAN <- as.matrix(data.frame(MEAN))
dimnames(MEAN) <- NULL
RSQ <- as.matrix(data.frame(RSQ))
dimnames(RSQ) <- NULL
} else {
lmFunc <- function(x, y) {
x <- as.matrix(x)
n <- ncol(x)
x <- as.vector(x)
y <- rep(y, n)
components <- c("coefficients", "fitted.values",
"residuals", "model")
lm(y ~ x, na.action = na.omit)[components]
}
biasCoefs <- matrix(NA, 2, nForecasts)
dimnames(biasCoefs) <- list(NULL, ensMemNames)
for (labX in uniqueEX) {
I <- namEX == labX
fit <- lmFunc(ensembleForecasts(ensembleData)[,
I, drop = FALSE], obs)
biasCoefs[, I] <- fit$coefficients
M <- !MISSmat[, I, drop = FALSE]
MEAN[, I][M] <- fit$fitted.values
RSQ[, I][M] <- fit$residuals
}
}
bad <- biasCoefs[2, ] < 0
if (any(bad)) print("biasCoefs < 0")
}, additive = {
if (nullEX) {
intcpt <- apply(obs - ensembleForecasts(ensembleData),
2, mean, na.rm = TRUE)
biasCoefs <- rbind(intcpt, 1)
dimnames(biasCoefs) <- NULL
MEAN <- sweep(ensembleForecasts(ensembleData), MARGIN = 2,
FUN = "+", STATS = intcpt)
dimnames(MEAN) <- NULL
RSQ <- obs - MEAN
dimnames(RSQ) <- NULL
} else {
intFunc <- function(x, y) {
x <- as.matrix(x)
n <- ncol(x)
x <- as.vector(x)
y <- rep(y, n)
mean(y - x, na.rm = TRUE)
}
intcpt <- rep(NA, nForecasts)
for (labX in uniqueEX) {
I <- namEX == labX
intcpt[I] <- intFunc(ensembleForecasts(ensembleData)[,
I, drop = FALSE], obs)
}
biasCoefs <- rbind(intcpt, 1)
dimnames(biasCoefs) <- NULL
MEAN <- sweep(ensembleForecasts(ensembleData), MARGIN = 2,
FUN = "+", STATS = intcpt)
dimnames(MEAN) <- NULL
RSQ <- obs - MEAN
dimnames(RSQ) <- NULL
}
}, none = {
biasCoefs <- matrix(c(0, 1), 2, nForecasts)
MEAN <- ensembleForecasts(ensembleData)
dimnames(MEAN) <- NULL
RSQ <- obs - MEAN
dimnames(RSQ) <- NULL
})
RSQ <- RSQ^2
z <- matrix(0, ncol = nForecasts, nrow = nObs)
nIter <- 0
loglik <- 0
if (control$equalVariance)
sd <- rep(sd, nForecasts)
while (TRUE) {
z <- as.matrix(data.frame(lapply(1:nForecasts, function(i,
y, mu, sd) dnorm(y, mean = mu[, i], sd = sd[i]),
y = obs, mu = MEAN, sd = rep(sd, length = nForecasts))))
z <- sweep(z, MARGIN = 2, FUN = "*", STATS = weights)
dimnames(z) <- list(dimnames(z)[[1]], ensMemNames)
zsum1 <- apply(z, 1, sum, na.rm = TRUE)
if (any(zsum1 == 0))
stop("parameter estimation fails")
z <- sweep(z, MARGIN = 1, FUN = "/", STATS = zsum1)
z[z < ZERO] <- 0
old <- loglik
loglik <- sum(log(zsum1))
zsum2 <- apply(z, 2, sum, na.rm = TRUE)
weights <- zsum2/sum(zsum2)
weights[weights < ZERO] <- 0
if (nullEX) {
if (control$equalVariance) {
sd <- sqrt(sum(z * RSQ, na.rm = TRUE)/sum(z,
na.rm = TRUE))
}
else {
sd <- sqrt(apply(z * RSQ, 2, sum, na.rm = TRUE)/zsum2)
}
}
else {
weights <- drop(crossprod(weights, matEX)[, namEX])
if (control$equalVariance) {
sd <- sqrt(sum((z * RSQ) %*% matEX, na.rm = TRUE)/sum(z %*%
matEX, na.rm = TRUE))
z <- sweep(z, MARGIN = 2, FUN = "/", STATS = nEX)
}
else {
sd <- sqrt(apply(z * RSQ, 2, sum, na.rm = TRUE)/zsum2)
sd <- drop(crossprod(sd, matEX))[namEX]
}
}
nIter <- nIter + 1
ERROR <- abs(loglik - old)/(1 + abs(loglik))
if (nIter > 1 && ERROR < control$tol)
break
if (nIter >= control$maxIter)
break
}
names(weights) <- ensMemNames
structure(list(biasCoefs = biasCoefs, sd = sd, weights = weights,
exchangeable = exchangeable,
nIter = nIter, loglikelihood = loglik, call = match.call()),
class = c("fitBMAnormal", "fitBMA"))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.