R/utilities.r

In aws: Adaptive Weights Smoothing

########## risk()

setMethod("risk", "aws", function(y, u = 0) {
  riskyhat(extract(y, "yhat")$yhat, u)
})
setMethod("risk", "awssegment", function(y, u = 0) {
  riskyhat(extract(y, "yhat")$yhat, u)
})


setMethod("risk", "kernsm", function(y, u = 0) {
  riskyhat(y@yhat, u)
})
setMethod("risk", "ICIsmooth", function(y, u = 0) {
  riskyhat(y@yhat, u)
})


setMethod("risk", "numeric", function(y, u = 0) {
  riskyhat(y, u)
})
setMethod("risk", "array", function(y, u = 0) {
  riskyhat(y, u)
})

########## extract()

setMethod("extract", "aws", function(x, what = "y") {
  what <- tolower(what)
  z <- list(NULL)
  if ("y" %in% what)
    z$y <- x@y
  if ("yhat" %in% what)
    z$yhat <- if (x@degree == 0)
      drop(x@theta)
  else
    switch(length(x@dy), x@theta[, 1], x@theta[, , 1], x@theta[, , , 1])
  if ("x" %in% what)
    z$x <- x@x
  if ("sigma2" %in% what)
    z$sigma2 <- x@sigma2
  if ("ni" %in% what)
    z$ni <- x@ni
  if ("mask" %in% what)
    z$mask <- x@mask
  if (length(z) > 1)
    z <- z[-1]
  invisible(z)
})
setMethod("extract", "awssegment", function(x, what = "y") {
  what <- tolower(what)
  z <- list(NULL)
  if ("y" %in% what)
    z$y <- x@y
  if ("yhat" %in% what)
    z$yhat <- drop(x@theta)
  if ("segment" %in% what)
    z$segment <- x@segment
  if ("x" %in% what)
    z$x <- x@x
  if ("sigma2" %in% what)
    z$sigma2 <- x@sigma2
  if ("ni" %in% what)
    z$ni <- x@ni
  if ("mask" %in% what)
    z$mask <- x@mask
  if (length(z) > 1)
    z <- z[-1]
  invisible(z)
})
setMethod("extract", "kernsm", function(x, what = "y") {
  what <- tolower(what)
  z <- list(NULL)
  if ("y" %in% what)
    z$y <- x@y
  if ("yhat" %in% what)
    z$yhat <- x@yhat
  if ("vred" %in% what)
    z$vred <- x@vred
  if ("vhat" %in% what)
    z$vhat <- (median(abs(diff(x@y)) / .9538)) ^ 2 / x@vred
  if (length(z) > 1)
    z <- z[-1]
  invisible(z)
})
setMethod("extract", "ICIsmooth", function(x, what = "y") {
  what <- tolower(what)
  z <- list(NULL)
  if ("y" %in% what)
    z$y <- x@y
  if ("yhat" %in% what)
    z$yhat <- x@yhat
  if ("vhat" %in% what)
    z$vhat <- x@vhat
  if ("vred" %in% what)
    z$vred <- x@sigma ^ 2 / x@vhat
  if ("hbest" %in% what)
    z$hbest <- x@hbest
  if (length(z) > 1)
    z <- z[-1]
  invisible(z)
})
################################################################
#                                                              #
# Section for summary(), print(), plot() functions (generic)   #
#                                                              #
################################################################

setMethod("show", "aws",
          function(object) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Max. Bandwidth       :",
                signif(object@hmax, 3),
                "  degree=",
                object@degree,
                "\n")
            cat("  Lambda               :", object@lambda, "  (ladjust=", object@ladjust, ")\n")
            cat("  Slots", slotNames(object), "\n")
            invisible(NULL)
          })
setMethod("print", "aws",
          function(x) {
            cat("  Object of class", class(x), "\n")
            cat("  Generated by calls   :\n")
            print(x@call)
            cat("  Dimension            :", paste(x@dy, collapse = "x"), "\n")
            cat("  Max. Bandwidth       :",
                signif(x@hmax, 3),
                "  degree=",
                x@degree,
                "\n")
            cat("  Lambda               :", x@lambda, "  (ladjust=", x@ladjust, ")\n")
            cat("  Slots", slotNames(x), "\n")
            invisible(NULL)
          })
setMethod("summary", "aws",
          function(object, ...) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Max. Bandwidth       :",
                signif(object@hmax, 3),
                "  degree=",
                object@degree,
                "\n")
            cat("  Lambda               :", object@lambda, "  (ladjust=", object@ladjust, ")\n")
            cat("  mean sum of weights  :", mean(object@ni), "\n")
            cat("\n")
            invisible(NULL)
          })
setMethod("plot", "aws",
          function(x,
                   what = "yhat",
                   col = grey(0:255 / 255),
                   zind = NULL) {
            if (length(what) > 1) {
              warning("what should be a single character")
              return(invisible(NULL))
            }
            if (!(what %in% c("data", "yhat", "ni", "mask"))) {
              warning("illegal value of what")
              return(invisible(NULL))
            }
            dy <- x@dy
            d <- length(dy)
            img <- extract(x, what)[[1]]
            if (d == 3) {
              if (is.null(zind))
                zind <- 1:x@dy[3]
              nview <- length(zind)
              nv1 <- as.integer(sqrt((nview + 2) / 1.5))
              nv2 <- ((nview + nv1 - 1) %/% nv1)
              par(
                mfrow = c(nv1, nv2),
                mar = c(2, 2, 2, .1),
                mgp = c(1, 1, 0)
              )
              for (i in zind) {
                image(
                  1:dy[1],
                  1:dy[2],
                  img[, , i],
                  col = col,
                  zlim = range(img[, , zind]),
                  xlab = "",
                  ylab = ""
                )
                title(paste("slice", i))
              }
            }
            if (d == 2) {
              image(1:dy[1],
                    1:dy[2],
                    img,
                    col = col,
                    xlab = "",
                    ylab = "")
              title(switch(
                what,
                "data" = "data",
                "yhat" = paste("aws hmax=", signif(x@hmax, 3)),
                "ni" = "sum of weights",
                "mask" = "mask"
              ))
            }
            if (d == 1) {
              if (what %in% c("data", "yhat")) {
                plot(x@y)
                lines(x@theta[, 1], col = 2)
                title(paste("Data and AWS for hmax=", signif(x@hmax, 3)))
              } else {
                plot(img)
                title(switch(what, "ni" = "sum of weights", "mask" = "mask"))
              }
            }
            invisible(NULL)
          })
setMethod("show", "awssegment",
          function(object) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Max. Bandwidth       :", signif(object@hmax, 3), "\n")
            cat("  Lambda               :", object@lambda, "  (ladjust=", object@ladjust, ")\n")
            cat("  Slots", slotNames(object), "\n")
            invisible(NULL)
          })
setMethod("print", "awssegment",
          function(x) {
            cat("  Object of class", class(x), "\n")
            cat("  Generated by calls   :\n")
            print(x@call)
            cat("  Dimension            :", paste(x@dy, collapse = "x"), "\n")
            cat("  Max. Bandwidth       :", signif(x@hmax, 3), "\n")
            cat("  Lambda               :", x@lambda, "  (ladjust=", x@ladjust, ")\n")
            cat("  Slots", slotNames(x), "\n")
            invisible(NULL)
          })
setMethod("summary", "awssegment",
          function(object, ...) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Max. Bandwidth       :", signif(object@hmax, 3), "\n")
            cat("  Lambda               :", object@lambda, "  (ladjust=", object@ladjust, ")\n")
            cat("  Size of segments:     ", paste(c("-1:", "0:", "1:"), table(object@segment)), "\n")
            cat("  mean sum of weights  :", mean(object@ni), "\n")
            cat("\n")
            invisible(NULL)
          })


setMethod("plot", "awssegment",
          function(x,
                   what = "segment",
                   col = grey(0:255 / 255),
                   zind = NULL) {
            if (length(what) > 1) {
              warning("what should be a single character")
              return(invisible(NULL))
            }
            if (!(what %in% c("data", "yhat", "ni", "mask", "segment"))) {
              warning("illegal value of what")
              return(invisible(NULL))
            }
            dy <- x@dy
            d <- length(dy)
            img <- extract(x, what)[[1]]
            if (d == 3) {
              if (is.null(zind))
                zind <- 1:x@dy[3]
              nview <- length(zind)
              nv1 <- as.integer(sqrt((nview + 2) / 1.5))
              nv2 <- ((nview + nv1 - 1) %/% nv1)
              par(
                mfrow = c(nv1, nv2),
                mar = c(2, 2, 2, .1),
                mgp = c(1, 1, 0)
              )
              for (i in zind) {
                image(
                  1:dy[1],
                  1:dy[2],
                  img[, , i],
                  col = col,
                  zlim = range(img[, , zind]),
                  xlab = "",
                  ylab = ""
                )
                title(paste("slice", i))
              }
            }
            if (d == 2) {
              image(1:dy[1],
                    1:dy[2],
                    img,
                    col = col,
                    xlab = "",
                    ylab = "")
              title(switch(
                what,
                "data" = "data",
                "yhat" = paste(
                  "awssegm hmax=",
                  signif(x@hmax, 3),
                  "ni" = "sum of weights",
                  "mask" = "mask",
                  "segment" = "segments"
                )
              ))
            }
            if (d == 1) {
              if (what %in% c("data", "yhat")) {
                plot(x@y)
                lines(x@theta[, 1], col = 2)
                title(paste("Data and AWS for hmax=", signif(x@hmax, 3)))
              } else {
                plot(img)
                title(switch(
                  what,
                  "ni" = "sum of weights",
                  "yhat" = "yhat",
                  "mask" = "mask",
                  "segment" = "segments"
                ))
              }
            }
            invisible(NULL)
          })
setMethod("show", "kernsm",
          function(object) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Bandwidth            :", object@h, "\n")
            cat("  Slots", slotNames(object), "\n")
            invisible(NULL)
          })
setMethod("print", "kernsm",
          function(x) {
            cat("  Object of class", class(x), "\n")
            cat("  Generated by calls   :\n")
            print(x@call)
            cat("  Dimension            :", paste(x@dy, collapse = "x"), "\n")
            cat("  Bandwidth            :", x@h, "\n")
            cat("  Slots", slotNames(x), "\n")
            invisible(NULL)
          })
setMethod("summary", "kernsm",
          function(object, ...) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Bandwidth            :", object@h, "  derivatives=", object@m, "\n")
            cat("  Variance reduction   :", mean(object@vred), "\n")
            cat("\n")
            invisible(NULL)
          })
setMethod("plot", "kernsm",
          function(x,
                   what = "yhat",
                   col = grey(0:255 / 255),
                   zind = NULL) {
            if (length(what) > 1) {
              warning("what should be a single character")
              return(invisible(NULL))
            }
            if (!(what %in% c("data", "yhat", "vred"))) {
              warning("illegal value of what")
              return(invisible(NULL))
            }
            dy <- x@dy
            d <- length(dy)
            img <- extract(x, what)[[1]]
            if (d == 3) {
              if (is.null(zind))
                zind <- 1:x@dy[3]
              nview <- length(zind)
              nv1 <- as.integer(sqrt((nview + 2) / 1.5))
              nv2 <- ((nview + nv1 - 1) %/% nv1)
              par(
                mfrow = c(nv1, nv2),
                mar = c(2, 2, 2, .1),
                mgp = c(1, 1, 0)
              )
              for (i in zind) {
                image(
                  1:dy[1],
                  1:dy[2],
                  img[, , i],
                  col = col,
                  zlim = range(img[, , zind]),
                  xlab = "",
                  ylab = ""
                )
                title(paste("slice", i))
              }
            }
            if (d == 2) {
              image(1:dy[1],
                    1:dy[2],
                    img,
                    col = col,
                    xlab = "",
                    ylab = "")
              title(switch(
                what,
                "data" = "data",
                "yhat" = paste("ksmooth h=", x@h),
                "vred" = "variance reduction"
              ))
            }
            if (d == 1) {
              if (what %in% c("data", "yhat")) {
                plot(x@y)
                lines(x@yhat, col = 2)
                title(paste("Data and kernsm for h=", x@h))
              } else {
                plot(img)
                title(switch(what, "vred" = "variance reduction"))
              }
            }
            invisible(NULL)
          })
setMethod("show", "ICIsmooth",
          function(object) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Bandwidth            :", object@hmax, " hinc:", object@hinc, "\n")
            cat("  Threshold            :", object@thresh, "\n")
            if (object@nsector > 1)
              cat(
                "  # of sectors          :",
                object@nsector,
                if (object@sector > 0)
                  paste("sector", object@sector),
                "symmetric:",
                object@symmetric,
                "\n"
              )
            cat("  Slots", slotNames(object), "\n")
            invisible(NULL)
          })
setMethod("print", "ICIsmooth",
          function(x) {
            cat("  Object of class", class(x), "\n")
            cat("  Generated by calls   :\n")
            print(x@call)
            cat("  Dimension            :", paste(x@dy, collapse = "x"), "\n")
            cat("  Bandwidth            :", x@hmax, " hinc:", x@hinc, "\n")
            cat("  Threshold            :", x@thresh, "\n")
            if (x@nsector > 1)
              cat(
                "  # of sectors          :",
                x@nsector,
                if (x@sector > 0)
                  paste("sector", x@sector),
                "symmetric:",
                x@symmetric,
                "\n"
              )
            cat("  Slots", slotNames(x), "\n")
            invisible(NULL)
          })
setMethod("summary", "ICIsmooth",
          function(object, ...) {
            cat("  Object of class", class(object), "\n")
            cat("  Generated by calls   :\n")
            print(object@call)
            cat("  Dimension            :", paste(object@dy, collapse = "x"), "\n")
            cat("  Bandwidth            :",
                object@hmax,
                "  hinc:",
                object@hinc,
                "  derivatives=",
                object@m,
                "\n")
            cat("  Threshold            :", object@thresh, "\n")
            if (object@nsector > 1)
              cat(
                "  # of sectors          :",
                object@nsector,
                if (object@sector > 0)
                  paste("sector",
                        object@sector),
                "symmetric:",
                object@symmetric,
                "\n"
              )
            cat("  Mean variance         :", mean(object@vhat), "\n")
            cat("\n")
            invisible(NULL)
          })
setMethod("plot", "ICIsmooth",
          function(x,
                   what = "yhat",
                   col = grey(0:255 / 255),
                   zind = NULL,
                   ...) {
            if (length(what) > 1) {
              warning("what should be a single character")
              return(invisible(NULL))
            }
            if (!(what %in% c("data", "yhat", "vhat", "vred", "hbest"))) {
              warning("illegal value of what")
              return(invisible(NULL))
            }
            dy <- x@dy
            d <- length(dy)
            img <- extract(x, what)[[1]]
            if (d == 3) {
              if (is.null(zind))
                zind <- 1:x@dy[3]
              nview <- length(zind)
              nv1 <- as.integer(sqrt((nview + 2) / 1.5))
              nv2 <- ((nview + nv1 - 1) %/% nv1)
              par(
                mfrow = c(nv1, nv2),
                mar = c(2, 2, 2, .1),
                mgp = c(1, 1, 0)
              )
              cat("displaying", switch(
                what,
                "data" = "data",
                "yhat" = paste("ICIsmooth hmax=", x@hmax),
                "vhat" = "variance",
                "vred" = "variance reduction",
                "hbest" = "optimal bandwidth"
              ), "for ICIsmooth")
              for (i in zind) {
                image(
                  1:dy[1],
                  1:dy[2],
                  img[, , i],
                  col = col,
                  zlim = range(img[, , zind]),
                  xlab = "",
                  ylab = ""
                )
                title(paste("slice", i))
              }
            }
            if (d == 2) {
              image(1:dy[1],
                    1:dy[2],
                    img,
                    col = col,
                    xlab = "",
                    ylab = "")
              title(switch(
                what,
                "data" = "data",
                "yhat" = paste("ICIsmooth hmax=", x@hmax),
                "vhat" = "variance",
                "vred" = "variance reduction",
                "hbest" = "optimal bandwidth"
              ))
            }
            if (d == 1) {
              if (what %in% c("data", "yhat")) {
                plot(x@y, ylab = what)
                lines(x@yhat, col = 2)
                title(paste("Data and ICIsmooth for hmax=", x@hmax))
              } else {
                plot(img, ylab = what)
                title(switch(
                  what,
                  "vhat" = "variance",
                  "vred" = "variance reduction",
                  "hbest" = "optimal bandwidth"
                ))
              }
            }
            invisible(NULL)
          })