diff --git a/R/whittaker.mortalityTable.R b/R/whittaker.mortalityTable.R
index 335a3748e76f4b8fb462034dd5de2985a62942eb..f6f09982bca117f4f43db026989407fa9165daa0 100644
--- a/R/whittaker.mortalityTable.R
+++ b/R/whittaker.mortalityTable.R
@@ -41,20 +41,39 @@ whittaker.mortalityTable = function(table, lambda = 10, d = 2, name.postfix = ",
table@name = paste0(table@name, name.postfix)
}
+# browser()
probs = table@deathProbs
+ orig.probs = probs
ages = table@ages
if (missing(weights) || is.null(weights)) {
if (is.na(table@exposures) || is.null(table@exposures)) {
- weights = table@exposures
- } else {
weights = rep(1, length(ages))
+ } else {
+ weights = table@exposures
}
}
- # Missing values are always interpolated, i.e. assigned weight 0:
- weights = weights * is.na(probs)
+ # Missing values are always interpolated, i.e. assigned weight 0; Similarly,
+ # ignore zero probabilities (cause problems with log)
+ weights = weights * (!is.na(probs) & (probs > 0))
+ weights[is.na(weights)] = 0
+ if (sum(probs > 0, na.rm = TRUE) < d) {
+ warning("Table '", table@name, "' does not have at least ", d, " finite, non-zero probabilities. Unable to graduate. The original probabilities will be retained.")
+ return(table)
+ }
+
+ # Normalize the weights to sum 1
+ weights = weights / sum(weights)
+
+ # We cannot pass NA to whittaker, as this will result in all-NA graduated values.
+ # However, if prob==NA, then weight was already set to 0, anyway
+ probs[is.na(probs)] = 0
+ probs.smooth = exp(whittaker.interpolate(log(probs), lambda = lambda, d = d, weights = weights))
- probs.smooth = exp(whittaker.interpolate(log(probs[!NApos]), lambda = lambda, d = d, weights = weights))
+ # Do not extrapolate probabilities, so set all ages below the first and
+ # above the last raw probability to NA
+ probsToClear = (cumsum(!is.na(orig.probs)) == 0) | (rev(cumsum(rev(!is.na(orig.probs)))) == 0)
+ probs.smooth[probsToClear] = NA_real_
table@deathProbs = probs.smooth
table
@@ -64,10 +83,12 @@ whittaker.mortalityTable = function(table, lambda = 10, d = 2, name.postfix = ",
whittaker.interpolate = function(y, lambda = 1600, d = 2, weights = rep(1, length(y))) {
m <- length(y)
E <- eye(m)
+ weights = weights * is.finite(y) # non-finite or missing values in y get zero weight
+ y[!is.finite(y)] = 0
W <- diag(weights)
D <- diff(E, lag = 1, differences = d)# - r * diff(E, lab = 1, differences = d - 1)
B <- W + (lambda * t(D) %*% D)
- z <- solve(B, W %*% y)
+ z <- solve(B, as.vector(W %*% y))
return(z)
}
#