<jupyter_start><jupyter_text>Marginal Gaussianization* Author: <NAME>* Email: <EMAIL>In this demonstration, we will show how we can do the marginal Gaussianization on a 2D dataset using the Histogram transformation and Inverse CDF Gaussian distribution.<jupyter_code>import os, sys
cwd = os.getcwd()
# sys.path.insert(0, f"{cwd}/../")
sys.path.insert(0, "/home/emmanuel/code/rbig")

from rbig.data import ToyData
from rbig.transform.gaussianization import MarginalGaussianization
# from rbig.transform.gaussianization import HistogramGaussianization, KDEGaussianization
from rbig.transform import InverseGaussCDF

import numpy as np
from scipy import stats

# Plot Functions
import matplotlib.pyplot as plt
import seaborn as sns

sns.reset_defaults()
#sns.set_style('whitegrid')
#sns.set_context('talk')
sns.set_context(context='talk',font_scale=0.7)
%matplotlib inline

%load_ext autoreload
%autoreload 2<jupyter_output><empty_output><jupyter_text>DataFor this example, we are looking at a 2D dataset.<jupyter_code>def plot_2d_joint(data, color='blue', title='Original Data'):
    
    fig  = plt.figure(figsize=(5, 5))

    g = sns.jointplot(x=data[:, 0], y=data[:, 1], kind='hex', color=color)
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.suptitle(title)
    plt.tight_layout()
    plt.show()
    
def plot_prob(data, probs, title='Probabilities'):
    
    fig, ax = plt.subplots()

    h = ax.scatter(data[:, 0], data[:, 1], s=1, c=probs, cmap='Reds')

    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    cbar = plt.colorbar(h, )
    ax.set_title(title)
    plt.show()
seed = 123
rng = np.random.RandomState(seed=seed)
dataset = 'rbig'
n_samples = 10_000
n_features = 2
noise = 0.25
random_state=1
clusters = 2


data = ToyData(
    dataset=dataset, 
    n_samples=n_samples, 
    n_features=n_features, 
    noise=noise, 
    random_state=random_state,
    clusters=clusters,
).generate_samples()

X = data[:, 0]
Y = data[:, 1]

plot_2d_joint(data, title='Original Data')<jupyter_output><empty_output><jupyter_text>Uniformization Transformation<jupyter_code>from rbig.transform.uniformization import HistogramUniformization, KDEUniformization, MarginalUniformization
# from rbig.density.histogram import ScipyHistogram, QuantileHistogram
# from rbig.den<jupyter_output><empty_output><jupyter_text>Initialize Uniformization Algorithm<jupyter_code># INITIALIZE UNIFORMIZATION ALGORITHM
#===
# uniform_clf = HistogramUniformization(bins=100, support_extension=10, alpha=1e-4, n_quantiles=None)
uniform_clf = KDEUniformization(n_quantiles=50, method='fft')
# density_clf = KDEScipy(n_quantiles=50, bw_method='scott', support_extension=10)
# density_clf = KDESklearn(n_quantiles=100, support_extension=10)<jupyter_output><empty_output><jupyter_text>Add it to Marginal Transformation Algorithm<jupyter_code>mg_uniformizer = MarginalUniformization(uniform_clf)


mg_uniformizer.fit(data)
X_trans = mg_uniformizer.transform(data)

plot_2d_joint(X_trans, title='Transformed Data')
data_approx = mg_uniformizer.inverse_transform(X_trans)

plot_2d_joint(data_approx, title='Transformed Data')
X_ldj = mg_uniformizer.log_abs_det_jacobian(data)

plot_2d_joint(X_ldj, title='Transformed Data')
plot_2d_joint(np.exp(X_ldj), title='Transformed Data')
plot_prob(data, X_ldj.sum(-1), title='Log Probabilities')
plot_prob(data, np.exp(X_ldj.sum(-1)), title='Probabilities')<jupyter_output><empty_output><jupyter_text>Marginal Gaussinization<jupyter_code>from rbig.transform.uniformization import HistogramUniformization, KDEUniformization, MarginalUniformization
from rbig.transform.gaussianization import MarginalGaussianization
uniform_clf = HistogramUniformization(bins=100, support_extension=10, alpha=1e-4, n_quantiles=None)
uniform_clf = KDEUniformization(n_quantiles=50, method='fft', )

mg_gaussianizer = MarginalGaussianization(uniform_clf)


mg_gaussianizer.fit(data)
X_trans = mg_gaussianizer.transform(data)

plot_2d_joint(X_trans, title='Transformed Data')
data_approx = mg_gaussianizer.inverse_transform(X_trans)

plot_2d_joint(data_approx, title='Transformed Data')
X_ldj = mg_gaussianizer.log_abs_det_jacobian(data)

plot_2d_joint(X_ldj, title='Transformed Data')
plot_2d_joint(np.exp(X_ldj), title='Transformed Data')
X_lprob = mg_gaussianizer.score_samples(data)

plot_prob(data, X_lprob, title='Log Probabilities')
plot_prob(data, np.exp(X_lprob), title='Probabilities')<jupyter_output><empty_output><jupyter_text>Negative Log Likelihood<jupyter_code>X_nll = mg_gaussianizer.score(data,)

print(f"Negative Log-Likelihood Score: {X_nll:.4f}")<jupyter_output>Negative Log-Likelihood Score: -2.8415<jupyter_text>Marginal Histogram TransformationSo, for this transformation, we are going to transform our data from the current distribution to a marginally Gaussian distribution and then perform a rotation. In theory, if we do enough of these, we will eventually convert to a Gaussian distribution.<jupyter_code># parameters
nbins = 1_000   # number of bins to do the histogram transform
alpha = 1e-05   # adds some regularization (noise)
support_extension = 10


# initialize the transformer
mg_transformer = HistogramGaussianization(
    nbins=nbins,
    alpha=alpha
)

# fit the transformer to the data
mg_transformer.fit(data);<jupyter_output><empty_output><jupyter_text>1. Forward TransformationFor this transformation, we will be applying the following:$$\Psi(\mathbf{x}) = \Phi^{-1}(\mathbf{x})$$where $\Phi^{-1}(\cdot)$ is the inverse CDF of the Gaussian distribution.<jupyter_code>data_trans = mg_transformer.transform(data)

plot_2d_joint(data_trans, title='Transformed Data')<jupyter_output><empty_output><jupyter_text>So clearly we can see that the transformation works. Both of the marginals are Gaussian distributed..  2. Inverse TransformationFor this step, we will apply the inverse transformation:$$\Psi^{-1}(\mathbf{x}) = \Phi \left( \mathbf{x} \right)$$where $\Phi(\cdot)$ is the CDF of the Gaussian distribution.<jupyter_code>data_approx = mg_transformer.inverse_transform(data_trans)

# check that its more or less equal
np.testing.assert_array_almost_equal(data_approx, data, decimal=1e-5)<jupyter_output><empty_output><jupyter_text>We see that this transformation is very close to the original. In fact, it's close to approximately 1e-5 decimal places. The errors will definitely stem from the boundaries.<jupyter_code># Plot results
plot_2d_joint(data_approx, title='Inverse Transformed Data')<jupyter_output><empty_output><jupyter_text>Log Absolute Determinant JacobianUsing the derivative of inverse-functions theorem, we can calculate the derivative like so:$$\nabla_\mathbf{x} \Phi^{-1}(\mathbf{x}) = \frac{1}{\phi (\Phi^{-1} (x)) }$$where $\phi(\cdot)$ is the PDF of the Gaussian distribution. Taking the log of these terms gives us:$$ \log \nabla_\mathbf{x} \Phi^{-1}(\mathbf{x}) = - \log \phi (\Phi^{-1} (x))$$<jupyter_code>X_slogdet = mg_transformer.log_abs_det_jacobian(data)

print(X_slogdet.min(), X_slogdet.max())
print(np.exp(X_slogdet).min(), np.exp(X_slogdet).max())
# plot the gradients
plot_2d_joint(np.exp(X_slogdet), title='Jacobian Data')<jupyter_output><empty_output><jupyter_text>Log Probability$$\log p_\theta(\mathbf{x}) = \log p_\theta \left( \mathbf{z} \right) + \log \left| \nabla_\mathbf{x} \mathbf{z} \right|$$where $\mathbf{z} = \Psi(\mathbf{x})$<jupyter_code># score samples 
log_prob = mg_transformer.score_samples(data)
# score samples 
log_prob = mg_transformer.score_samples(data)

plot_prob(data, log_prob, title='Log Probabilities')<jupyter_output><empty_output><jupyter_text>ProbabilityThis is the same as above but without the log scale:$$p_\theta(\mathbf{x}) = p_\theta \left( \mathbf{z} \right) \left| \nabla_\mathbf{x} \mathbf{z} \right|$$where $\mathbf{z} = \Psi(\mathbf{x})$<jupyter_code>plot_prob(data, np.exp(log_prob), title='Probabilities')<jupyter_output><empty_output><jupyter_text>Negative Log-LikelihoodWe need to take the expected value (mean) of all log probabilities.$$\text{nll} = \frac{1}{N} \sum_{n=1}^{N} \log p_\theta(\mathbf{x})$$<jupyter_code>score = mg_transformer.score(data)

print(f"Negative Log-Likelihood Score: {score:.4f}")<jupyter_output>Negative Log-Likelihood Score: -2.0724
StarcoderdataJupyterStructuredCleanDedup