R/xtms2eqns.R

In ipeglim: Imprecise Inferential Framework on Statistical Reasoning

#' @title 
#' Converting convex extreme points to a set of linear inequality constraints
#'
#' @description 
#' Once the information about the extreme points in a convex hull, 
#' a set of unique linear inequality constraints exist.
#' The function \code{xtms2eqns} searches for this set.
#' 
#' @param obj a list of extreme points
#' @param show a logical for visualization 
#'  
#' @author Chel Hee Lee <\email{gnustats@@gmail.com}>
#' @export 
xtms2eqns <- function(obj, show=FALSE, type=c("eqns", "circle")){
	stopifnot(is.list(obj))
	if(missing(type)) type <- "eqns"
	
	# object is originally taken from 'cmfit'
	# thus, xy <- obj$xtms
	xy1 <- xy <- obj$xtms
	xy0 <- do.call(rbind, xy)
	xy1[[length(xy)+1]] <- xy[[1]]
	xy1 <- do.call(rbind, xy1)
	cc <- colMeans(xy0)
	ui <- ci <- list()
	
	xlim <- range(xy0[,1])
	ylim <- range(xy0[,2])
	if(show) plot(0, xlim=c(xlim[1], xlim[2]), ylim=c(ylim[1], ylim[2]))

	for(i in 1:length(xy)){
		x1 <- xy1[i,1]
		y1 <- xy1[i,2]
		x2 <- xy1[i+1,1]
		y2 <- xy1[i+1,2]
		b <- (y2-y1)/(x2-x1)
		a <- -b*x1 + y1
		val <- a+b*cc[1]
		
		if(is.na(val) | is.nan(val)){
		
			if(x1 == x2){
			if(x1 <= cc[1]){
				ui[[i]] <- c(1, 0)
				ci[[i]] <- x1
				if (show) abline(v=x1, col="green")
				if (show) polygon(rbind(c(x1,y1), c(x2,y2), cc), col="green")
			} else {
				ui[[i]] <- c(-1, 0)
				ci[[i]] <- -x1
				if (show) abline(v=x1, col="orange")
				if (show) polygon(rbind(c(x1,y1), c(x2,y2), cc), col="orange")
			}
			next
			}
			
			if(y1==y2){
			if(y1 <= cc[2]){
				ui[[i]] <- c(0, -1)
				ci[[i]] <- y1
				if (show) abline(v=x1, col="green")
				if (show) polygon(rbind(c(x1,y1), c(x2,y2), cc), col="green")
			} else {
				ui[[i]] <- c(0, 1)
				ci[[i]] <- y1
				if (show) abline(v=x1, col="orange")
				if (show) polygon(rbind(c(x1,y1), c(x2,y2), cc), col="orange")
			}
			next
			}
		}
		if ( val <= cc[2] ) {
			ui[[i]] <- c(1,-b)
			ci[[i]] <- a
			if (show) abline(a=a,b=b, col="red")
			if (show) polygon(rbind(c(x1,y1), c(x2,y2), cc), col=alpha("red", 0.3))
		} else {
			if(type=="circle") { 
				ui[[i]] <- -c(1,-b) # gamma prior- (circle is working; however, general polygon is not working)
				ci[[i]] <- -a
			} else {
				ui[[i]] <- c(-1,-b) # normal prior - polygon works
				ci[[i]] <- -a
			}
			if (show) abline(a=a,b=b, col="blue")
			if (show) polygon(rbind(c(x1,y1), c(x2,y2), cc), col=alpha("blue", 0.3))
		}
		if (show) polygon(xy0[chull(xy0),], col=alpha("azure2", 0.2), border=alpha("darkblue",0.5))
		if (show) points(x=cc[1], y=cc[2])
		if (show) points(xy0)
	}
	rvals <- list(cc=cc, ui=do.call(rbind, ui), ci=do.call(c, ci))
}
NULL