title: isolation forest algorithm in python
body: <p>I am trying to reproduce the algorithm described in the Isolation Forest paper in python. <a href="http://cs.nju.edu.cn/zhouzh/zhouzh.files/publication/icdm08b.pdf?q=isolation" rel="nofollow">http://cs.nju.edu.cn/zhouzh/zhouzh.files/publication/icdm08b.pdf?q=isolation</a></p>
  
  <p>This is my current code:</p>
  
  <pre><code>import numpy as np
  import sklearn as sk
  import matplotlib.pyplot as plt
  import pandas as pd
  from sklearn.decomposition import PCA
  
  
  def _h(i):
      return np.log(i) + 0.5772156649 
  
  
  def _c(n):
      if n &gt; 2:
          h = _h(n-1)
          return 2*h - 2*(n - 1)/n
      if n == 2:
          return 1
      else:
          return 0
  
  
  def _anomaly_score(dict_scores, n_samples):
      score = np.array([np.mean(dict_scores[k]) for k in dict_scores.keys()])
      score = -score/_c(n_samples)
  
      return 2**score
  
  
  def _split_data(X):
      ''' split the data in the left and right nodes ''' 
      n_samples, n_columns = X.shape
      n_features = n_columns - 1
  
      feature_id = np.random.randint(low=0, high=n_features-1)
      feature = X[:, feature_id]
      split_value = np.random.choice(feature)
      left_X = X[feature &lt;= split_value]
      right_X = X[feature &gt; split_value]
      return left_X, right_X, feature_id, split_value
  
  
  def iTree(X, add_index=False, max_depth = np.inf):            
      ''' construct an isolation tree and returns the number of step required
      to isolate an element. A column of index is added to the input matrix X if  
      add_index=True. This column is required in the algorithm. ''' 
  
      n_split = {} 
      def iterate(X, count = 0):
  
          n_samples, n_columns = X.shape
          n_features = n_columns - 1
  
          if count &gt; max_depth:
              for index in X[:,-1]:
                  n_split[index] = count
              return
  
          if n_samples == 1:
              index = X[0, n_columns-1]
              n_split[index] = count
              return 
          else:
              lX, rX, feature_id, split_value = _split_data(X)
              # Uncomment the print to visualize a draft of 
              # the construction of the tree
              #print(lX[:,-1], rX[:,-1], feature_id, split_value, n_split)
              n_samples_lX, _ = lX.shape
              n_samples_rX, _ = rX.shape
              if n_samples_lX &gt; 0:
                  iterate(lX, count+1)
              if n_samples_rX &gt;0:
                  iterate(rX, count+1)
  
      if add_index:
          n_samples, _ = X.shape
          X = np.c_[X, range(n_samples)]
  
      iterate(X)
      return n_split
  
  
  class iForest():
      ''' Class to construct the isolation forest.
  
      -n_estimators: is the number of trees in the forest,
  
      -sample_size: is the bootstrap parameter used during the construction
      of the forest,
  
      -add_index: adds a column of index to the matrix X. This is required and 
      add_index can be set to False only if the last column of X contains 
      already indeces.
  
      -max_depth: is the maximum depth of each tree
      '''
      def __init__(self, n_estimators=20, sample_size=None, add_index = True, 
                   max_depth = 100):
          self.n_estimators = n_estimators
          self.sample_size = sample_size
          self.add_index = add_index
          self.max_depth = max_depth
          return 
  
      def fit(self, X):
          n_samples, n_features = X.shape
          if self.sample_size == None:
              self.sample_size = int(n_samples/2)
  
          if self.add_index:
              X = np.c_[X, range(n_samples)]
  
  
          trees = [iTree(X[np.random.choice(n_samples, 
                                            self.sample_size, 
                                            replace=False)],
                         max_depth=self.max_depth) 
                   for i in range(self.n_estimators)]
  
          self.all_anomaly_score_ = {k:None for k in range(n_samples)}
          for k in self.all_anomaly_score_.keys():
              self.all_anomaly_score_[k] = np.array([tree[k] 
                                                     for tree in trees 
                                                     if k in tree])
  
          self.anomaly_score_ = _anomaly_score(self.all_anomaly_score_, n_samples)
          return self
  </code></pre>
  
  <p>The main part of the code is the iTree function that returns a dictionary with the number of steps required to isolate each sample. </p>
  
  <p>A column of index is attached to the input matrix <code>X</code> in order to make easier to understand what are the samples in each node. </p>
  
  <p>When I compare the result obtained with my code and the ones obtained with the isolation forest available for R I get different results.</p>
  
  <p>Consider for example that stackloss dataset :</p>
  
  <pre><code>data = pd.read_csv("stackloss.csv")
  X = data.as_matrix()[:, 1:]
  
  max_depth = 100
  itree = iTree(X, add_index=True, max_depth=max_depth) #example of isolation tree 
  
  
  iforest = iForest(n_estimators=1, max_depth=max_depth, sample_size=21) # isolation forest 
  iforest.fit(X)
  
  sol = np.argsort(iforest.anomaly_score_)
  #correct sol = [10  5  4  8 12  9 11 17  6 19  7 14 13 15 18  3 20 16  2  1  0]
  </code></pre>
  
  <p><code>sol</code> is often different by correct solution obtained with the R software.
  <a href="https://r-forge.r-project.org/projects/iforest/" rel="nofollow">https://r-forge.r-project.org/projects/iforest/</a></p>
  
  <p>The correct solution in R has been obtained with:</p>
  
  <pre><code>&gt; tr = IsolationTrees(stackloss,ntree = 100000,hlim = 100, rFactor = 1)
  &gt; as = AnomalyScore(stackloss, tr)
  &gt; order(as$outF)
   [1] 11  6  5  9 13 10 12 18  7 20  8 15 14 16 19  4 21 17  3  2  1
  &gt; order(as$outF)-1
   [1] 10  5  4  8 12  9 11 17  6 19  7 14 13 15 18  3 20 16  2  1  0
  &gt; 
  </code></pre>
  
  <p>Where is the mistake ?</p>