def nominal_splits(x, y, x_vals, y_dist, separate_max):
    """
    Function uses heuristic to find best binary split of nominal values. Heuristic is described in (1) and it is
    originally defined for binary classes. We extend it to work with multiple classes by comparing label with least
    samples to others.

    x: numpy array - nominal feature
    y: numpy array - label
    x_vals: numpy array - unique nominal values of x
    y_dist: dictionary - distribution of labels

    Reference:
    (1) Classification and Regression Trees by Breiman, Friedman, Olshen, and Stone, pages 101- 102. 
    """
    # select a label with least samples
    y_val = max(y_dist, key=y_dist.get) if separate_max else min(y_dist, key=y_dist.get)

    prior = y_dist[y_val] / float(len(y))  # prior distribution of selected label

    values, dist, splits = [], [], []
    for x_val in x_vals:  # for every unique nominal value
        dist.append(Counter(y[x == x_val]))  # distribution of labels at selected nominal value
        splits.append(x_val)
        suma = sum([prior * dist[-1][y_key] for y_key in y_dist.keys()])
        # estimate probability
        values.append(prior * dist[-1][y_val] / float(suma))
    indices = np.array(values).argsort()[::-1]

    # distributions and splits are sorted according to probabilities
    return np.array(dist)[indices], np.array(splits)[indices].tolist()