R/iprior.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

#' @title
#' Imprecise Prior Specification
#'
#' @description 
#' The function name \code{iprior} stands for \code{imprecise prior}.
#' A set of $k$ linear inequality constraints is used for modelling 
#' \code{iprior}.
#' The details of this specification are given under \sQuote{Details}.
#' 
#' @param obj an object of class \code{imprecise} produced from 
#' \code{\link{model}}.
#'
#' @param eqns 
#' a list with components of \code{lhs} and \code{rhs};
#' See \sQuote{Details} for more information.
#' 
#' @param circle 
#' a list with components of \code{x}, \code{y}, \code{z}, \code{r}, and
#' \code{len}.  See \sQuote{Details} for more information.
#'
#' @param x
#' a numeric vector.
#'
#' @param mat
#' a numeric matrix.
#'
#' @param verbose
#' a logical value; Be more verbose about the process by displaying messages.
#' Defaults to \code{FALSE}.
#' 
#' @details
#' Modelling a prior ignorance is the second stage on the imprecise inferential
#' framework. 
#' \code{stage} has the value \code{"iprior"} that is the name of 
#' environment called.
#' 
#' A convex hull is defined by \code{lhs \%*\% x >= rhs}, where \code{lhs} is
#' a constraint matrix \eqn{k \times p}{k x p} and \code{rhs} is a constraint 
#' vector of length \eqn{k} from \eqn{k} linear inequality constraints. 
#' The function \code{iprior} searches for the solutions \code{xtms} of 
#' this system which are the extreme points of that convex hull.
#'
#' Depending on the numbers of dimensions \eqn{p} and constraints \eqn{k},
#' a shape of this convex hull is determined. 
#' For example, a box-constrained convex hull is made on \eqn{p=2} dimensional
#' space with \eqn{k=4} constraints. 
#' When \eqn{k} has an infinite number, a nutural shape defaults to a 
#' \code{circle} on \eqn{p=2} dimensional space.
#'
#' @note
#' Conventionally, the current version classifies a shape of convex hull into
#' five categories: \code{eqns2d}, \code{eqns3d}, \code{eqns4more}, 
#' \code{circle2d}, and \code{sphere3d} for visualization.
#' 
#' It is a difficult task to sarch for extreme points in more than three
#' dimensional cases.
#' The function \code{convhulln} in the package \code{geometry} is a good 
#' replacement of \code{iprior}.
#' (See more details about \code{convhulln} at 
#' \url{http://geometry.r-forge.r-project.org/}).
#'
#' @return
#' An object of class \code{imprecise} with the components:
#' \item{stage}{the environment name called.}
#' \item{xtms}{the list of solutions of the system supplied.}
#' 
#' @examples
#' \dontrun{
#' lc0 <- list(lhs=rbind(diag(2), -diag(2)), rhs=c(rep(-1,4)))
#' cmfit <- iprior(obj=m2fit, eqns=lc0)
#' plot(cmfit)
#' }
#' 
#' @keywords 
#' convex hull, geometry, extreme points, constrained optimization
#'
#' @seealso 
#' \code{\link{chull}}, \code{\link{convhulln}}, 
#'
#' @references
#' Lee (2013) ``Imprecise inferential framework'', PhD thesis.
#'
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>
#' @export
iprior <- function(obj, eqns=list(), circle=list(), x, mat, verbose=FALSE, ...) {
	UseMethod(generic="iprior")
}
NULL

#' @rdname iprior
#' @method iprior default
#' @S3method iprior default
iprior.default <- function(obj, eqns=list(), circle=list(), x, mat, verbose=FALSE, ...){

	mc <- match.call()
	obj <- list(stage=NA)
  stopifnot(any(!missing(eqns), !missing(circle)))
  stopifnot(missing(x), missing(mat))

  eq <- list(lhs=NULL, rhs=NULL)
  
  if (all(names(eqns) %in% names(eq))) {
    eq[names(eqns)] <- eqns
  } else {
    stop("Unknown names in ", sQuote("eqns"))
  }
  
  if (!missing(eqns)) {
    stopifnot(is.list(eqns))
    stopifnot(is.matrix(eq$lhs), is.numeric(eq$lhs))
    stopifnot(is.vector(eq$rhs), is.numeric(eq$rhs))
    obj$constraints <- eq
    
    if (ncol(eq$lhs)==2) {
      obj$m0shape <- "eqns2d" # box-constrained
    } else if (ncol(eq$lhs)==3) {
      obj$m0shape <- "eqns3d" # cube
    } else {
      obj$m0shape <- "eqns4more" # higher dimensions
    }    
  }
  
  ccl <- list(x=0, y=0, z=NULL, r=1, len=15)
  
  if (all(names(circle) %in% names(ccl))) {
    ccl[names(circle)] <- circle
  } else {
    stop("Unknown names in ", sQuote("circle"))
  }
  
  if (!missing(circle)) {
    stopifnot(is.list(circle))
    stopifnot(is.numeric(ccl$x), length(ccl$x)==1)
    stopifnot(is.numeric(ccl$y), length(ccl$y)==1)
    stopifnot(is.numeric(ccl$r), length(ccl$r)==1)
    stopifnot(is.numeric(ccl$len), length(ccl$len)==1)
    obj$constraints <- ccl

    if (is.null(ccl$z)) {
      m0shape <- "circle2d"
    } else {
      m0shape <- "sphere3d"
    }
    obj$m0shape <- m0shape
  }
  
  if (!missing(x)) {
    stopifnot(is.numeric(x), is.vector(x))
    obj$constraints <- x
  }
  
  if (!missing(mat)) {
    stopifnot(is.numeric(mat), is.matrix(mat))
    obj$constraints <- mat
  }
  
  if (verbose) {
    message(mc[[1]], ":\nInput sanity check .................... PASS!")
  }
  
  if (!missing(eqns) ) {
    input.type <- "eqns" # 2d? or 3d?
    xtms <- s4xtms.eqns(x=eq)
  } else if (!missing(circle)) {
    input.type <- "circle" # 2d" or 3d?
    xtms <- s4xtms.circle(x=ccl)
  } else if (!missing(x)) {
    input.type <- "scalar"
  } else {
    input.type <- "matrix"
  }
  if (verbose) {
    message("Input type: ", sQuote(input.type))
  }
  
  # make sure if a set of extreme points are in the form of 'list'
  xtms <- as.list(as.data.frame(t(xtms)))
  names(xtms) <- gsub("V", "x", names(xtms))
  xtms <- lapply(xtms, function(x) {
		names(x) <- paste("V", seq_len(length(x)), sep="")
		return(x)
	})
	obj$xtms <- xtms
	class(obj) <- c("imprecise")
  invisible(obj)
  
}
NULL

#' @rdname iprior
#' @method iprior imprecise
#' @S3method iprior imprecise
iprior.imprecise <- function(obj, eqns=list(), circle=list(), x, mat, 
                             verbose=FALSE, ...){
                             
  # naming convention
  xreg <- obj$xreg
  stage <- obj$stage 
  mc <- match.call()
  
  if (stage != "model") {
    stop(sQuote("iprior()"), " is follwed by ", sQuote("model()"), "\n",
         "in the context of imprecise inferential framework.")
  }
	obj$stage <- "iprior"
  
  stopifnot(!missing(obj))
  stopifnot(any(!missing(eqns), !missing(circle), !missing(x), !missing(mat)))

  eq <- list(lhs=NULL, rhs=NULL)
  
  if (all(names(eqns) %in% names(eq))) {
    eq[names(eqns)] <- eqns
  } else {
    stop("Unknown names in ", sQuote("eqns"))
  }
  
  if (!missing(eqns)) {
  
    stopifnot(is.list(eqns))
    stopifnot(is.matrix(eq$lhs), is.numeric(eq$lhs))
    stopifnot(is.vector(eq$rhs), is.numeric(eq$rhs))
    obj$constraints <- eq
    
    if (ncol(eq$lhs)==2) {
      obj$m0shape <- "eqns2d" # box-constrained
    } else if (ncol(eq$lhs)==3) {
      obj$m0shape <- "eqns3d" # cube
    } else {
      obj$m0shape <- "eqns4more" # higher dimensions
    }
  }
  
  ccl <- list(x=0, y=0, z=NULL, r=1, len=15)
  
  if (all(names(circle) %in% names(ccl))) {
    ccl[names(circle)] <- circle
  } else {
    stop("Unknown names in ", sQuote("circle"))
  }
  
  if (!missing(circle)) {
    stopifnot(is.list(circle))
    stopifnot(is.numeric(ccl$x), length(ccl$x)==1)
    stopifnot(is.numeric(ccl$y), length(ccl$y)==1)
    stopifnot(is.numeric(ccl$r), length(ccl$r)==1)
    stopifnot(is.numeric(ccl$len), length(ccl$len)==1)
    obj$constraints <- ccl

    if (is.null(ccl$z)) {
      m0shape <- "circle2d"
    } else {
      m0shape <- "sphere3d"
    }
    obj$m0shape <- m0shape
  }
  
  if (!missing(x)) {
    stopifnot(is.numeric(x), is.vector(x))
    obj$constraints <- x
  }
  
  if (!missing(mat)) {
    stopifnot(is.numeric(mat), is.matrix(mat))
    obj$constraints <- mat
  }
  
  if (verbose) {
    message(mc[[1]], ":\nInput sanity check .................... PASS!")
  }
  
  if (!missing(eqns) ) {
    input.type <- "eqns"
  } else if (!missing(circle)) {
    input.type <- "circle"
  } else if (!missing(x)) {
    input.type <- "scalar"
  } else {
    input.type <- "matrix"
  }
  if (verbose) {
    message("Input type: ", sQuote(input.type))
  }
  
  if(FALSE){  ##########################################################
  fn.x <- function(...){
    obj$m0shape <- "scalar"
    xtms <- x
    return(xtms)
  }
  
  fn.mat <- function(...){
    obj$m0shape <- "matrix"
    
    stopifnot(ncol(mat)>=2)
    
    if (ncol(mat)==2) {
      xtms <- mat[chull(mat),]
    } 
    
    if (ncol(mat)>=3) {
      xtms <- mat[unique(geometry::convhulln(mat)),]
    }
    
    if (verbose) {
      message("Extreme points are found")
    }
    return(xtms)
  }
  } ####################################################################
  
  # search for extremes
  # xtms <- switch(input.type, 
  #  "eqns"=fn.eqns(...),
  #  "circle"=fn.circle(...),
  #  "scalar"=fn.x(...),
  #  "matrix"=fn.mat(...))

  if (!missing(eqns) ) {
    input.type <- "eqns" # 2d? or 3d?
    xtms <- s4xtms.eqns(x=eq)
  } else if (!missing(circle)) {
    input.type <- "circle" # 2d" or 3d?
    xtms <- s4xtms.circle(x=ccl)
  } else if (!missing(x)) {
    input.type <- "scalar"
  } else {
    input.type <- "matrix"
  }
  if (verbose) {
    message("Input type: ", sQuote(input.type))
  }
  
  # make sure if a set of extreme points are in the form of 'list'
  xtms <- as.list(as.data.frame(t(xtms)))
  names(xtms) <- gsub("V", "x", names(xtms))
  
  if (xreg) { 
    xtms <- lapply(xtms, function(x){ 
      names(x) <- colnames(obj$X)
      return(x)
    })
  }
  obj$xtms <- xtms
    
  invisible(obj)
}
NULL