https://github.com/cran/ensembleBMA
Tip revision: fcc61119d42912dce1599f2208b60a6a719e36aa authored by Chris Fraley on 14 August 2008, 00:00:00 UTC
version 3.0-4.1
version 3.0-4.1
Tip revision: fcc6111
fitBMAnormal.R
`fitBMAnormal` <-
function(ensembleData, control = controlBMAnormal(),
exchangeable = NULL)
{
if (is.null(exchangeable)) exchangeable <- ensembleGroups(ensembleData)
if (length(unique(exchangeable)) == length(exchangeable))
exchangeable <- NULL
if (!(nullX <- 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)
}
# remove instances missing all forecasts or obs
M <- apply(ensembleForecasts(ensembleData), 1, function(z) all(is.na(z)))
M <- M | is.na(ensembleVerifObs(ensembleData))
ensembleData <- ensembleData[!M,]
if (is.null(obs <- ensembleVerifObs(ensembleData)))
stop("verification observations required")
nObs <- ensembleNobs(ensembleData)
ensMemNames <- ensembleMemberLabels(ensembleData)
nForecasts <- length(ensMemNames)
if(is.null(sd <- control$start$sd)) sd <- sd(obs)
weights <- if (is.null(control$start$weights) || any(is.na(control$start.weights))) 1 else control$start$weights
if (length(weights) == 1) weights <- rep(weights,nForecasts)
weights <- weights/sum(weights)
weights <- pmax(weights,1.e-4)
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
# bias correction (constraints)
switch(control$biasCorrection,
regression = {
if (nullX) {
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 (nullX) {
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
}
)
## may have missing values
RSQ <- RSQ^2
z <- matrix(0, ncol=nForecasts, nrow=nObs)
nIter <- 0
loglik <- 0
if (control$equalVariance) sd <- rep(sd, nForecasts)
## cat("\n")
while (TRUE) # EM
{
## MEAN may have missing values
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)
z <- sweep( z, MARGIN = 1, FUN = "/", STATS = zsum1)
old <- loglik
loglik <- sum(log(zsum1))
zsum2 <- apply(z, 2, sum, na.rm = TRUE)
weights <- zsum2/sum(zsum2)
if (nullX) {
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 <- sapply(split(weights,namEX), mean)[namEX]
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)
##
## sd <- sqrt(apply(z*RSQ,2,sum,na.rm=TRUE)/apply(z,2,sum,na.rm=TRUE))
## sd <- mean(sd)
##
}
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))
## cat("", nIter)
## cat(" ", c(loglik, ERROR, min(weights)), "\n")
if (nIter > 1 && ERROR < control$eps) break
if (nIter >= control$maxIter) break
}
names(weights) <- ensMemNames
structure(
list(biasCoefs = biasCoefs, sd = sd, weights = weights,
exchangeable = exchangeable, nIter = nIter,
loglikelihood = loglik),
class = "fitBMAnormal")
}