Learning Multivariate Model Free Distributions

In this tutorial, we will go through the ideas, math and implementation that resulted in Joint Probability Trees (JPTs). JPTs are non-parametric, multivariate, smooth, deterministic and decomposable distributions.

Historic Origin

In 2020, Daniel Nyga started to think about what happens to decision trees when you take away more and more features. A decision tree usually solves a discriminating task; \(argmax_Y \,P(Y|X)\) where \(X = <x_1, \cdots x_M>\) is some feature vector.

Decision Trees partition the feature-space and assign the mode of Y to each partition. For a full feature vector, this kind of inference will always result in exactly on partition, where the prediction can be made.

As soon as features are missing, multiple partitions are possible. Hence, a mixture like voting is done. The first enhancement to this voting was to weight the votes of the partitions by the probability of being in such a partition. These probabilities can easily be obtained during learning by the frequencies of each partition in the training data. However, not every partition is equally likely under arbitrary values. Hence, the partitions were further enhanced by independent, univariate distributions in for each variable in each partition. For discrete variables, frequencies easily get these distributions again. The final problem was the functional form of the continuous distributions. Due to the infinite support of Gaussian distributions, these were unfeasible for finite partitions. Furthermore, not every variable is normally distributed. The flexible answer to this problem was the Nyga distribution, which is explained in Learning univariate, model-free and deterministic distributions

Without knowing it, Daniel had forged a robust, fast and scalable probabilistic circuit.

Theoretical Background

Learning in JPTs is done by a modified version of the C4.5 algorithm. The C4.5 algorithm is a decision tree induction algorithm, which optimizes the information gained by each split.

Let’s investigate the process more formally. Let \(\mathcal{D} = \{\boldsymbol{d_1}, \boldsymbol{d_2}, \cdots, \boldsymbol{d_n}\}\) be a dataset of \(n\) samples, where each sample is an \(m\)-dimensional vector.

For every \(m_{k}\) dimension, we evaluate every possible split by measuring the impurity of the resulting partitions on the remaining variables. The impurity is measured by the variance \(\mathbb{V}\) of the partition if the variable is continuous or ordinal, and by the Gini impurity \(G\) otherwise.

\[\begin{align*} G(\mathcal{D}_{m_j}) &= 1 - \sum_{state \, \in \, domain(m_{j})} P(state) ^ 2 \\ \mathbb{V}(\mathcal{D}_{m_j}) &= \mathbb{E}\left[ D_{m_j}^2 \right] - \mathbb{E}\left[ D_{m_j} \right]^2\\ Impurity(\mathcal{D}) &= \frac{1}{m} \sum_{m_{k} \in m} \begin{cases} \mathbb{V}(\mathcal{D}_{m_k}), \text{if variable is ordinal or continuous}\\ G(\mathcal{D}_{m_k}), \text{otherwise} \end{cases} \end{align*}\]

The information gain \(IG\) is the difference between the impurity of the dataset and the impurity of the partitioned dataset. The split with the highest information gain is chosen.

\[\begin{align*} IG_{m_k,condition} (m_k, split) &= Impurity(\mathcal{D}_{\setminus m_k}|split) + Impurity(\mathcal{D}_{\setminus m_k}| \neg split)\\ & argmax_{1 \cdots m, split} IG_{m_k,condition} (m_k, split) \end{align*}\]

Practical Application

First, let’s generate a random dataset, consisting of two bi-variate Gaussian’s.

import plotly
plotly.offline.init_notebook_mode()

import plotly.graph_objects as go
import numpy as np
import pandas as pd

np.random.seed(69)

distribution_1 = np.random.multivariate_normal(np.zeros((2,)), np.eye(2, 2), size=(2000,))
distribution_2 = np.random.multivariate_normal(np.zeros((2,)), np.eye(2, 2), size=(3000,))
distribution_2 += [3, 4]
distribution_2 = distribution_2 @ [[1, 2], [1, 0]]
dataset = np.concatenate((distribution_1, distribution_2))
dataset = pd.DataFrame(dataset, columns=["x", "y"])

fig = go.Figure()
fig.add_trace(go.Scatter(x=distribution_1[:, 0], y=distribution_1[:, 1], mode="markers", name="First Component"))
fig.add_trace(go.Scatter(x=distribution_2[:, 0], y=distribution_2[:, 1], mode="markers", name="Second Component"))
fig.update_layout(title="Dataset generated by a mixture of two Multivariate Gaussian's.")
fig.show()

Next, we have to overload the c45 method of the JPT class to see the learning progress.

from probabilistic_model.learning.jpt.variables import infer_variables_from_dataframe
from probabilistic_model.learning.jpt.jpt import JPT

class TutorialJPT(JPT):
    
    fig: go.Figure
    """
    The figure to plot the decision criteria into.
    """
    
    def c45(self, data: np.ndarray, start: int, end: int, depth: int):
        """
        Construct a DecisionNode or DecomposableProductNode from the data.

        :param data: The data to calculate the impurity from.
        :param start: Starting index in the data.
        :param end: Ending index in the data.
        :param depth: The current depth of the induction
        :return: The constructed decision tree node
        """

        number_of_samples = end - start
        root = self.root
        # if the inducing in this step results in inadmissible nodes, skip the impurity calculation
        if depth >= self.max_depth or number_of_samples < 2 * self.min_samples_leaf:
            max_gain = -float("inf")
        else:
            max_gain = self.impurity.compute_best_split(start, end)

        # if the max gain is insufficient
        if max_gain <= self.min_impurity_improvement:

            # create decomposable product node
            leaf_node = self.create_leaf_node(data[self.indices[start:end]])
            weight = number_of_samples / len(data)
            root.add_subcircuit(leaf_node, np.log(weight))

            if self.keep_sample_indices:
                leaf_node.sample_indices = self.indices[start:end]

            # terminate the induction
            return

        # if the max gain is sufficient
        split_pos = self.impurity.best_split_pos
        split_var_idx = self.impurity.best_var
        other_axis_index = 1 - split_var_idx
        split_value = (data[self.indices[start + split_pos], split_var_idx] + data[self.indices[start + split_pos + 1], split_var_idx]) / 2
        
        
        trace_data = np.zeros((2,2))
        minimal_value = data[self.indices[start:end], other_axis_index].min()
        maximal_value = data[self.indices[start:end], other_axis_index].max()

        trace_data[:, split_var_idx] = np.array([split_value, split_value]).T
        trace_data[:, other_axis_index] = np.array([minimal_value, maximal_value]).T
        self.fig.add_trace(go.Scatter(x=trace_data[:, 0], y=trace_data[:, 1], mode="lines", name="Decision Criterion"))
        

        # increase the depth
        new_depth = depth + 1

        # append the new induction steps
        self.c45queue.append((data, start, start + split_pos + 1, new_depth))
        self.c45queue.append((data, start + split_pos + 1, end, new_depth))

Next, let’s define, and fit the model. Plot the decision criteria and the resulting distribution.

variables = infer_variables_from_dataframe(dataset, scale_continuous_types=False, min_likelihood_improvement=2.5)
model = TutorialJPT(variables, min_impurity_improvement=0.05, min_samples_leaf=500)
model.fig = fig
pc = model.fit(dataset)
model.fig.show()

The lines labeled as “Decision Criterion” are the decision criteria that were learned by the C4.5 algorithm. They are ordered in the same way as they are induced in the distribution.

Finally, we plot resulting leaf distributions.

figure = go.Figure(pc.plot(1000, surface=True))
figure.show()

  0%|          | 0/7 [00:00<?, ?it/s]

100%|██████████| 7/7 [00:00<00:00, 92.02it/s]