added hpb and bsf functions

Author: David Dyjack

NAMESPACE

@@ -1,5 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
+export(hpb)
 export(hpe)
+importFrom(stats,dist)
 importFrom(stats,median)
 importFrom(stats,quantile)

R/.bsf.R.swp

@@ -0,0 +1,34 @@
+#' Bootstrap function (bsf) for H+ Estimate
+#'
+#' Estimates H+ for a random subset of the supplied data
+#' to be replicated within the main hpb() function
+#' Samples the observations (rows) of d and 
+#' the entries of l to generate a single bootstrap estimate.
+#'
+#' @param d (data, nxm)
+#' @param l (identities, nx1)
+#' @param r (numeric, 1x1)
+#' @param q (list, pre-computed by hpb) 
+#' @return One iteration of H+ bootstrap estimate
+#' @importFrom stats dist
+#' @keywords internal
+
+bsf <- function(d,l,s,q){
+  #is this a general solution to sampling from classes in a balanced way?
+  #q <- table(l)/length(l)
+  #p <- as.vector(sapply(names(q), function(x) sample(which(l==x),round(s/length(q)))))
+  p <- as.vector(sapply(q, function(x) sample(x,round(s/length(q)))))
+  lab <- l[p]
+  dis <- as.matrix(dist(d[p,]))
+  dis <- dis[upper.tri(dis)]
+  ind <- sapply(lab, function(x) x==lab)
+  ind <- ind[upper.tri(ind)]
+  iw <- which(ind)
+  ib <- which(!ind)
+  dw <- dis[iw]
+  db <- dis[ib]
+  sp <- sum(sapply(dw, function(x) sum(x>db)))
+  hp <- sp / (as.numeric(length(dw))*as.numeric(length(db)))
+  return(hp)
+
+}

R/.hpb.R.swp

@@ -0,0 +1,51 @@
+#' @title H+ discordance estimation with bootstrapping
+#'
+#' @description Estimates the H+ discordance metric 
+#' given data Dat and a label vector L. 
+#' Bootstrapping is performed over r iterations with s points sampled per iterate.
+#'
+#' @param D numeric matrix or data frame with observations in rows (nxm)
+#' @param L numeric vector containing a set of length n
+#' @param r numeric number of bootstrap iterations
+#' @param s numeric  pre-bootstrap sample size 
+#' 
+#' @return list, h is the estimated H+ value.
+#' @return gamma1 and gamma2 are plausible ranges for what % of A (or Dw)
+#' @return are strictly greater than B (or Db)
+#' @export
+#' 
+#' @importFrom stats median quantile
+#' 
+#' @examples
+#' a <- sapply(1:500, function(i) rnorm(n=50, mean=0))
+#' b <- sapply(1:500, function(i) rnorm(n=50, mean=0))
+#' x <- t(cbind(a,b))
+#' l <- c(rep(0,500), rep(1,500))
+#' h <- hpb(D=x, L=l, r=30,s=50)
+#' 
+hpb <- function(D, L, r=30, s) {
+  dlflg <- missing(D) | missing(L)
+  if(dlflg){
+    stop("please provide both D and L")
+  }
+  tyflg <- !any(class(D) %in% c('matrix','data.frame')) & is.numeric(D)
+  if(tyflg){
+   stop("please ensure D is a numeric matrix or data.frame")
+  }
+  dmflg <- ! nrow(D) == length(L)
+  if(dmflg){
+    stop("Dimension mismatch for D and L")
+  }
+  srflg <- ! class(r) %in% c('numeric','integer') & class(s) %in% c('numeric','integer')
+  if(srflg){
+    stop("please ensure r and s are numeric or integer")
+  }
+  snflg <- s >= nrow(D)
+  if(snflg){
+    stop("s is too large, try ~1-5% of observations")
+  }
+  lbtab <- table(L)/length(L) #pre-compute this table / list
+  lblis <- lapply(names(lbtab), function(x) which(L==x))
+  hpb <- mean(replicate(r,bsf(d=D,l=L,s=s,q=lblis),T))
+  return(hpb)
+}

R/bsf.R

@@ -39,7 +39,7 @@ hpe <- function(A, B, D, L, p = 101, alg = "brute_force",alphas=F) {
   if (abflg & dlflg) {
     stop("please provide either (A and B) or (D and L)")
   } else if (!abflg) {
-    nmflg <- (!is.numeric(A) & !is.numeric(A))
+    nmflg <- (!is.numeric(A) & !is.numeric(B))
     if (nmflg) {
       stop("please ensure A B are numeric")
     } else {