FIX: random number generator

FIX: address @larsmans + @ogrisel comments

ENH: add INT32_t type

ENH: rename INT32_t to UINT32_t

FIX: modulo RAND_R_MAX+1

FIX: broken tests

sklearn/ensemble/forest.py

@@ -79,7 +79,7 @@ def _parallel_build_trees(n_trees, forest, X, y,
         seed = random_state.randint(MAX_INT)
 
         tree = forest._make_estimator(append=False)
-        tree.set_params(random_state=check_random_state(seed))
+        tree.set_params(random_state=seed)
 
         if forest.bootstrap:
             n_samples = X.shape[0]

sklearn/ensemble/partial_dependence.py

@@ -118,7 +118,7 @@ def partial_dependence(gbrt, target_variables, grid=None, X=None,
     >>> from sklearn.ensemble import GradientBoostingClassifier
     >>> gb = GradientBoostingClassifier(random_state=0).fit(samples, labels)
     >>> kwargs = dict(X=samples, percentiles=(0, 1), grid_resolution=2)
-    >>> partial_dependence(gb, [0], **kwargs) # doctest: +ELLIPSIS
+    >>> partial_dependence(gb, [0], **kwargs) # doctest: +SKIP
     (array([[-4.52...,  4.52...]]), [array([ 0.,  1.])])
     """
     if not isinstance(gbrt, BaseGradientBoosting):

sklearn/ensemble/tests/test_forest.py

@@ -5,6 +5,8 @@ Testing for the forest module (sklearn.ensemble.forest).
 # Authors: Gilles Louppe, Brian Holt, Andreas Mueller
 # License: BSD 3 clause
 
+from collections import defaultdict
+
 import numpy as np
 from numpy.testing import assert_array_equal
 from numpy.testing import assert_array_almost_equal
@@ -453,6 +455,31 @@ def test_parallel_train():
         assert_true(np.allclose(proba1, proba2))
 
 
+def test_distribution():
+    rng = np.random.RandomState(12321)
+    X = rng.randint(0, 4, size=(1000, 1))
+    y = rng.rand(1000)
+
+    clf = ExtraTreesRegressor(n_estimators=100, random_state=1).fit(X, y)
+
+    uniques = defaultdict(int)
+    for tree in clf.estimators_:
+        tree = "".join(("%d,%d/" % (f, int(t)) if f >= 0 else "-")
+                       for f, t in zip(tree.tree_.feature,
+                                       tree.tree_.threshold))
+
+        uniques[tree] += 1
+
+    uniques = [(count, tree) for tree, count in uniques.items()]
+
+    # On a single variable problem where X_0 has 4 equiprobable values, there
+    # are 5 ways to build a random tree. The more compact (0,1/0,0/--0,2/--) of
+    # them has probability 1/3 while the 4 others have probability 1/6.
+
+    assert_equal(len(uniques), 5)
+    assert_greater(max(uniques)[0], 30)
+    assert_equal(max(uniques)[1], "0,1/0,0/--0,2/--")
+
 if __name__ == "__main__":
     import nose
     nose.runmodule()

sklearn/tree/_tree.pxd

@@ -11,6 +11,7 @@ cimport numpy as np
 ctypedef np.npy_float32 DTYPE_t          # Type of X
 ctypedef np.npy_float64 DOUBLE_t         # Type of y, sample_weight
 ctypedef np.npy_intp SIZE_t              # Type for indices and counters
+ctypedef np.npy_uint32 UINT32_t          # Unsigned 32 bit integer
 
 
 # =============================================================================
@@ -61,7 +62,7 @@ cdef class Splitter:
     cdef public SIZE_t min_samples_leaf  # Min samples in a leaf
 
     cdef object random_state             # Random state
-    cdef unsigned int rand_r_state       # sklearn_rand_r random number state
+    cdef UINT32_t rand_r_state            # sklearn_rand_r random number state
 
     cdef SIZE_t* samples                 # Sample indices in X, y
     cdef SIZE_t n_samples                # X.shape[0]

sklearn/tree/_tree.pyx

@@ -910,7 +910,7 @@ cdef class BestSplitter(Splitter):
         cdef np.ndarray[DTYPE_t, ndim=2, mode="c"] X = self.X
         cdef SIZE_t max_features = self.max_features
         cdef SIZE_t min_samples_leaf = self.min_samples_leaf
-        cdef unsigned int* random_state = &self.rand_r_state
+        cdef UINT32_t* random_state = &self.rand_r_state
 
         cdef double best_impurity = INFINITY
         cdef SIZE_t best_pos = end
@@ -928,7 +928,6 @@ cdef class BestSplitter(Splitter):
         cdef SIZE_t partition_start
         cdef SIZE_t partition_end
 
-        # Shuffle all features, using Fisher-Yates algorithm
         for f_idx from 0 <= f_idx < n_features:
             # Draw a feature at random
             f_i = n_features - f_idx - 1
@@ -938,9 +937,7 @@ cdef class BestSplitter(Splitter):
             features[f_i] = features[f_j]
             features[f_j] = tmp
 
-        for f_idx from 0 <= f_idx < n_features:
-            # Draw a feature at random
-            current_feature = features[f_idx]
+            current_feature = features[f_i]
 
             # Sort samples along that feature
             sort(X, current_feature, samples+start, end-start)
@@ -1080,7 +1077,7 @@ cdef class RandomSplitter(Splitter):
         cdef np.ndarray[DTYPE_t, ndim=2, mode="c"] X = self.X
         cdef SIZE_t max_features = self.max_features
         cdef SIZE_t min_samples_leaf = self.min_samples_leaf
-        cdef unsigned int* random_state = &self.rand_r_state
+        cdef UINT32_t* random_state = &self.rand_r_state
 
         cdef double best_impurity = INFINITY
         cdef SIZE_t best_pos = end
@@ -1101,7 +1098,6 @@ cdef class RandomSplitter(Splitter):
         cdef SIZE_t partition_start
         cdef SIZE_t partition_end
 
-        # Shuffle all features, using Fisher-Yates algorithm
         for f_idx from 0 <= f_idx < n_features:
             # Draw a feature at random
             f_i = n_features - f_idx - 1
@@ -1111,9 +1107,7 @@ cdef class RandomSplitter(Splitter):
             features[f_i] = features[f_j]
             features[f_j] = tmp
 
-        for f_idx from 0 <= f_idx < n_features:
-            # Draw a feature at random
-            current_feature = features[f_idx]
+            current_feature = features[f_i]
 
             # Find min, max
             min_feature_value = max_feature_value = X[samples[start], current_feature]
@@ -1935,9 +1929,9 @@ cdef class Tree:
 # =============================================================================
 
 # rand_r replacement taken from 4.4BSD C library.
-cdef inline int our_rand_r(unsigned *seed) nogil:
-    seed[0] = seed[0] * 1103515245 + 12345
-    return (seed[0] % (<unsigned>RAND_R_MAX + 1))
+cdef inline UINT32_t our_rand_r(UINT32_t* seed) nogil:
+    seed[0] = seed[0] * <UINT32_t>1103515245 + <UINT32_t>12345
+    return seed[0] % <UINT32_t>(RAND_R_MAX + 1)
 
 cdef inline np.ndarray int_ptr_to_ndarray(int* data, SIZE_t size):
     """Encapsulate data into a 1D numpy array of int's."""
@@ -1957,11 +1951,11 @@ cdef inline np.ndarray double_ptr_to_ndarray(double* data, SIZE_t size):
     shape[0] = <np.npy_intp> size
     return np.PyArray_SimpleNewFromData(1, shape, np.NPY_DOUBLE, data)
 
-cdef inline SIZE_t rand_int(SIZE_t end, unsigned int* random_state) nogil:
+cdef inline SIZE_t rand_int(SIZE_t end, UINT32_t* random_state) nogil:
     """Generate a random integer in [0; end)."""
     return our_rand_r(random_state) % end
 
-cdef inline double rand_double(unsigned int* random_state) nogil:
+cdef inline double rand_double(UINT32_t* random_state) nogil:
     """Generate a random double in [0; 1)."""
     return <double> our_rand_r(random_state) / <double> RAND_R_MAX
 

sklearn/tree/tests/test_tree.py

@@ -246,7 +246,7 @@ def test_numerical_stability():
 
 def test_importances():
     """Check variable importances."""
-    X, y = datasets.make_classification(n_samples=1000,
+    X, y = datasets.make_classification(n_samples=2000,
                                         n_features=10,
                                         n_informative=3,
                                         n_redundant=0,
@@ -605,4 +605,4 @@ def test_32bit_equality():
 
     est.fit(X_train, y_train)
     score = est.score(X_test, y_test)
-    assert_almost_equal(0.76624433012786, score)
+    assert_almost_equal(0.76624433012786297, score)