Open-Category-Detection-with-PAC-Guaranteesp-Synthetic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Code for the Experiments in Open Category Detection with PAC Guarantees on Synthetic Datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(devtools)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(histogram)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "install_github(\"liusi2019/btloda\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(btloda)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(MASS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## fuction for getting quantile estimate tau^hat_q using raw ECDF method\n",
    "raw_cdf<- function(datab, datan, alpha, q){\n",
    "    trialv <- sort(c(datan, datab))\n",
    "    F.n <- ecdf(datan)\n",
    "    Fn  <- F.n(trialv)\n",
    "    F.b <- ecdf(datab)\n",
    "    Fb <- F.b(trialv)\n",
    "    Fa <- (Fn - (1-alpha)*Fb)/alpha \n",
    "    if (length(which(Fa <= q))==0){\n",
    "      index <- 1\n",
    "    }else{\n",
    "      index <-max(which(Fa <= q))\n",
    "    }\n",
    "    return(trialv[index])\n",
    "  }\n",
    "\n",
    "## function for getting quantile estimate tau^hat_q using isotonized ECDF method\n",
    "iso_cdf <- function(datab, datan, alpha, q){\n",
    "    trialv <- sort(c(datan, datab))\n",
    "    F.n <- ecdf(datan)\n",
    "    Fn  <- F.n(trialv)\n",
    "    F.b <- ecdf(datab)\n",
    "    Fb <- F.b(trialv)\n",
    "    Fa <- (Fn - (1-alpha)*Fb)/alpha \n",
    "    F.is = isoreg(Fa)$yf ## Computes the Isotonic Estimator of Fa\n",
    "    F.is[which(F.is<=0)]=0  \n",
    "    F.is[which(F.is>=1)]=1\n",
    "    if (length(which(F.is <= q))==0){\n",
    "      index.is <- 1\n",
    "    }else{\n",
    "      index.is <-max(which(F.is <= q))\n",
    "    }\n",
    "    return(trialv[index.is])\n",
    "  }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## take as input anomaly scores for nominal and mixture data sets, alpha, q, anomaly scores of the\n",
    "## big nominal and anomaly data set\n",
    "## provide as output the \n",
    "\n",
    "total_result <- function(datab, datan, alpha, q, score_nominal, score_anomaly){\n",
    "    ## estiamte from raw ECDF method\n",
    "    est1 <- raw_cdf(datab, datan, alpha, q)\n",
    "    ## estiamte from isotonized ECDF method\n",
    "    est2 <- iso_cdf(datab, datan, alpha, q)\n",
    "      \n",
    "    result = rep(0, 5)\n",
    "    result[1] <- mean(score_anomaly >= est1)## recall, raw ECDF\n",
    "    result[2] <- mean(score_nominal >= est1)## FPR, the proportion of nominal that are misclassified as alien, raw ECDF\n",
    "    result[3] <- mean(score_anomaly >= est2)## recall, isotonized ECDF\n",
    "    result[4] <- mean(score_nominal >= est2)## FPR, isotonized ECDF\n",
    "    result[5] <- mean(score_nominal >= quantile(score_anomaly, q)) ## ground truth about FPR\n",
    "    return(result)\n",
    "  }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The following cells are for the experiments on synthetic data sets.\n",    
    "### Note that for the two big datasets: big_normal and big_anomaly, we create them once and never regenerate them in the experiments;\n",
    "### But for the pair of small nominal and small mixture data sets:  data1 and data2, we regenerate them for every replicate of the experiments. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## create two big datsets: big_nominal, big_anomaly for calculating ground truth about FPR\n",
    "big_nominal = base::matrix(rnorm(20000*9, 0, 1), nrow = 20000, ncol = 9)\n",
    "big_anomaly = matrix(ncol = 9, nrow = 20000)\n",
    "vmat = matrix(0, ncol = 9, nrow = 9)\n",
    "diag(vmat) = 1\n",
    "for(i in (1:20000)){\n",
    "    center = rep(0, 9)\n",
    "    if(rbinom(1, 1, 0.4)==1){\n",
    "      center[sample(9, 3,replace = F)] = 3\n",
    "      }else{\n",
    "      center[sample(9, 4,replace = F)] = 3\n",
    "      }\n",
    "    big_anomaly[i,] = mvrnorm(1, center, vmat)\n",
    "    } \n",
    "\n",
    "## write them into .csv file for running iforest in linux\n",
    "write.csv(big_nominal, file = paste(\"big_nominal.csv\", sep = \"\"), row.names = FALSE)\n",
    "write.csv(big_anomaly, file = paste(\"big_anomaly.csv\", sep = \"\"), row.names = FALSE) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## set up the sample size n, alien proportion alpha, and q\n",
    "n = 1000\n",
    "vpro = 10\n",
    "alpha = vpro/100  ## anomaly proportion in the mixture data set\n",
    "q = 0.05 ## targetting on q quantile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##  generate nominal and mixture datasets, each of size n\n",
    "dat =  matrix(ncol = 9, nrow = 2*n)\n",
    "vmat = matrix(0, ncol = 9, nrow = 9)\n",
    "diag(vmat) = 1\n",
    "for(i in (1:n*alpha)){\n",
    "  center = rep(0, 9)\n",
    "  if(rbinom(1, 1, 0.4)==1){\n",
    "    center[sample(9, 3,replace = F)] = 3\n",
    "    }else{\n",
    "    center[sample(9, 4,replace = F)] = 3\n",
    "    }\n",
    "  dat[i,] = mvrnorm(1, center, vmat)\n",
    "  } \n",
    "    \n",
    "nnrow <- round(n * (2 - alpha))\n",
    "rvec <- rnorm(9 * nnrow, 0, 1)\n",
    "## the top alpha * n points are alien points \n",
    "dat[(round(n*alpha)+1):(2*n),] <- base::matrix(rvec, nrow = round(n*(2-alpha)), ncol = 9)\n",
    "    \n",
    "data1 <- dat[(n+1):(2*n),] #nominal data set\n",
    "data2 <- dat[1:n,] #mixture data set\n",
    "\n",
    "## write them into .csv file for running iforest using linux    \n",
    "write.csv(data1, file = paste(\"data1.csv\", sep = \"\"), row.names = FALSE)\n",
    "write.csv(data2, file = paste(\"data2.csv\", sep = \"\"), row.names = FALSE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The isolation forest implementation used here is provided by Tadesse Zemicheal and available from https://github.com/tadeze/osu_iforest .\n",
    "### The code in the following cell send command to run iforest in Linux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## grow iforest using the nominal data set, get out-of-bag anomaly scores for each nominal point\n",
    "## and save the iforest grown\n",
    "system(paste('./iforest','-i', paste('./data1.csv', sep = \"\"),'-o', paste('./depth1.csv', sep = \"\"),'-s', round(0.2*n), '-t 1000 -g -k -b', paste('./forest1.bin',sep = \"\")), wait = TRUE)\n",
    "score1 <- read.csv(paste('./depth1.csv', sep = \"\"), header = TRUE)\n",
    "datab <- as.numeric(unlist(score1))\n",
    "   \n",
    "## run the mixture data set through the saved iforest \n",
    "system(paste('./iforest','-i', paste('./data2.csv', sep = \"\") ,'-o',paste('./depth2.csv', sep = \"\"),'-g -f', paste('./forest1.bin',sep = \"\")), wait = TRUE)\n",
    "score2 <- read.csv(paste('./depth2.csv', sep = \"\"), header = TRUE)\n",
    "datan <- as.numeric(unlist(score2))\n",
    "\n",
    "### run the big nominal and anomaly datasets through the forest\n",
    "system(paste('./iforest','-i', paste('./big_anomaly','.csv', sep = \"\") ,'-o',paste('./depth_anomaly.csv', sep = \"\"),'-g -f', paste('./forest1.bin',sep = \"\")), wait = TRUE)\n",
    "depth_anomaly <- read.csv(paste('./depth_anomaly.csv', sep = \"\"), header = TRUE)\n",
    "score_anomaly <- as.numeric(unlist(depth_anomaly))\n",
    "    \n",
    "system(paste('./iforest','-i', paste('./big_nominal','.csv', sep = \"\") ,'-o',paste('./depth_nominal.csv', sep = \"\"),'-g -f', paste('./forest1.bin',sep = \"\")), wait = TRUE)\n",
    "depth_nominal <- read.csv(paste('./depth_nominal.csv', sep = \"\"), header = TRUE)\n",
    "score_nominal <- as.numeric(unlist(depth_nominal))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_iforest <- total_result(datab, datan, alpha, q, score_nominal, score_anomaly)\n",
    "output_iforest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This cell below uses R package btloda, which implements the bootstrap and out-of-bag version of LODA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## build projections and create histograms using bootstrap of the nominal data set\n",
    "## get out-of-bag anomaly scores for each nominal point\n",
    "bt_out = btloda(data1,sparsity=NA, maxk=1000, inf_replace = log(1e-09))\n",
    "\n",
    "bt_datab = bt_out$oob_nll \n",
    "bt_datan = get_neg_ll_all_hist(data2, bt_out$pvh$w, bt_out$pvh$hists, inf_replace = log(1e-09))\n",
    "\n",
    "## get anomaly scores for the big nominal and anomaly data sets using projections and histograms\n",
    "bt_score_nominal = get_neg_ll_all_hist(big_nominal, bt_out$pvh$w, bt_out$pvh$hists, inf_replace = log(1e-09))\n",
    "bt_score_anomaly = get_neg_ll_all_hist(big_anomaly, bt_out$pvh$w, bt_out$pvh$hists, inf_replace = log(1e-09))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "output_btloda = total_result(bt_datab, bt_datan, alpha, q, bt_score_nominal, bt_score_anomaly)\n",
    "output_btloda"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}