<jupyter_start><jupyter_text>Importance WeightingLast time out we detected that there was a [covariate shift](5.1-covariate-shift.ipynb) between the training and validation sets. Recall, that we are using the validation set as a proxy for the test set. Here we will try and correct for the covariate shift by reweighting training samples according to their **importance**, which is the ratio of test to training densities. *Tsuboi et al.* explain in [Direct density ratio estimation for large-scale covariate shift adaptation](https://pdfs.semanticscholar.org/6b0b/12fbd2455c37131ac1cedc5137a9fd7c3e53.pdf) that "Under covariate shift, standard learning methods such as maximum likelihood estimation are no longer consistent, i.e., they do not produce the optimal solution even when the number of training samples tends to infinity. Thus, there exists an estimation bias induced by covariate shift. It has been shown that the bias can be asymptotically canceled by weighting the log likelihood terms according to the **importance**:\begin{equation}w(x) = \frac{p_{test}(x)}{p_{train}(x)}\end{equation}where $p_{test}(x)$ and $p_{train}(x)$ are test and training input densities." Technically, there is a more rigorous definition of importance that accounts for class imbalance between the training and test sets, which multiplies the above ratio by the ratio of number of training samples to the number of test samples. For a derivation, see section 2.2.1 of [Density-ratio matching under the Bregman divergence:a unified framework of density-ratio estimation](https://www.ism.ac.jp/editsec/aism/pdf/10463_2011_Article_343.pdf) by *Sugiyama et al*. Or alternatively, read the [density ratio trick](http://blog.shakirm.com/2018/01/machine-learning-trick-of-the-day-7-density-ratio-trick/) on **The Spectator: Shakir's machine learning blog**.Importance weighting may still not be clear to you. *<NAME>*, Math Ph.D. who works in machine learning, [elucidates exactly what importance weighting does from a mathematical perspective](https://srome.github.io/Covariate-Shift,-i.e.-Why-Prediction-Quality-Can-Degrade-In-Production-and-How-To-Fix-It/). Essentially the sample weights applied to each training example, change the distribution with respect to which the expected training loss is taken, from the training distribution to the test distribution:\begin{equation}\DeclareMathOperator{\E}{\mathbb{E}}\E_{x \sim w(x)p_{train}(x)}loss(x) = \int loss(x)w(x)p_{train}(x)dx = \int loss(x)p_{test}(x)dx = \E_{x \sim p_{test}(x)}loss(x)\end{equation}As he puts it, we can see that the expected loss on the test distribution $p_{test}(x)$ is the same as using the sample weights $w(x)$ on the training set! (Notice, this is trivially the case when $p_{train}(x)=p_{test}(x)$.) Once these sample weights are known, they can be supplied during the training phase via the `sample_weight` parameter in the `fit()` method of an *sklearn* learn estimator.The importance is rarely known in practice, so the central question becomes, how to accurately estimate it from the training and test samples? It is a difficult problem that we will answer shortly. If you are still unconvinced as to why estimating the importance is critical when covariate shift occurs, read the excellent *bias-variance tradeoff* argument given in the footnote of page 1, in the paper by *Tsuboi et al*.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from src.externals.pykliep.pykliep import DensityRatioEstimator
from src.features.features_utils import convert_categoricals_to_numerical

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Reading in the DataFirst let's read in both sets of training and validation features and convert the categorical fields to a numerical form that is suitable for building machine learning models. We will also read in the target for the training data as it will be useful for visualization purposes.<jupyter_code>train_features = pd.read_csv('../data/processed/train-features.csv')
X_train = convert_categoricals_to_numerical(train_features)
X_train.head()
train_features_topics = pd.read_csv('../data/processed/train-features-topics.csv')
X_train_topics = convert_categoricals_to_numerical(train_features_topics)
X_train_topics.head()
train_target = pd.read_csv('../data/processed/train-target.csv')
train_target.head()
validation_features = pd.read_csv('../data/processed/validation-features.csv')
X_validation = convert_categoricals_to_numerical(validation_features)
X_validation.head()
validation_features_topics = pd.read_csv('../data/processed/validation-features-topics.csv')
X_validation_topics = convert_categoricals_to_numerical(validation_features_topics)
X_validation_topics.head()<jupyter_output><empty_output><jupyter_text>Kullback-Leibler Importance Estimation Procedure (KLIEP)The importance can be estimated by estimating the test and training densities separately and taking their ratio. However, this does not work well in practice as the ratio is sensitive. *Sugiyama et al.* developed a more stable method known as [Kullback-Leibler Importance Estimation Procedure (KLIEP)](https://www.ism.ac.jp/editsec/aism/pdf/060_4_0699.pdf) that directly estimates the importance without the need of density estimation. They find an importance estimate $\hat{w}(x)$ such that the [Kullback–Leibler divergence](https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence) from the true test input density $p_{test}(x)$ to its estimate $\hat{p}_{test}(x) = \hat{w}(x)p_{train}(x)$ is minimized. The optimization problem involved in KLIEP is convex and so the unique global solution can be obtained. Technical details of the KLIEP method are described in full in section 2 of the *Sugiyama et al.* paper. A succint summary, that we will paraphrase here, is given by *<NAME>*, Math Ph.D. who works in machine learning. The KLIEP algorithm defines an estimator of the importance as:\begin{equation}\hat{w}(x) = \sum_{x_{i} \in D_{test}} \alpha_i K(x, x_i, \sigma)\end{equation}where the $\alpha_i \in \mathbb{R}$ are parameters to be learned from data samples and $D_{test} \subset X_{test}$ is a random sample of test set vectors $X_{test}$. *Sugiyama et al.* choose $K$ to be the **Gaussian kernel**:\begin{equation}K(x, x_i, \sigma) = \exp \bigg( -\frac{|x-x_i|^2}{2 \sigma^{2}} \bigg)\end{equation}where the hyperparameter $\sigma \in \mathbb{R}$ is called the **width** of the **kernel**. The authors then define the score:\begin{equation}J = \frac{1}{|X_{test}|} \sum_{x \in X_{test}} \log \hat{w}(x)\end{equation}and prove that maximizing this value is equivalent to minimizing the Kullback-Leibler divergence between $p_{test}(x)$ and $\hat{w}(x)p_{train}(x)$. Basically, maximizing $J$ leads to the optimal sample weights for training. The KLIEP algorithm uses [gradient descent](https://en.wikipedia.org/wiki/Gradient_descent) to find the $\alpha_i$ that maximize $J$ given $D_{test}$ and $\sigma$ subject to the following constraints:\begin{equation}\alpha_i \ge 0 \\\int \hat{w}(x)p_{train}(x)dx = 1\end{equation}The first constraint arises from the fact that the sample weights must be non-negative and the second constraint from the definition of $\hat{w}(x)p_{train}(x)$ as a [probability density function](https://en.wikipedia.org/wiki/Probability_density_function).  Model SelectionThere are essentially two hyperparameters that need to be chosen in the KLIEP algorithm, the **kernel width** $\sigma$ and the **ratio** $\frac{D_{test}}{X_{test}}$. Model selection can be carried out in a principled manner through the **likelihood cross-validation** (LCV) procedure that is described by *Sugiyama et al*. It is very similar to the grid-search cross-validation that most of us are familiar with, except for a few notable differences:1. Split the test set $D_{test}$ into $k$-folds.2. For each test fold, train a model on the union of the training set and $k$-1 test folds. Evaluate the score $J$ on the $k$-th test fold.3. Select the parameters based on the best averaged score and re-train the model on the full training and test sets.The small test set size means that the computational costs of LCV are low, so we will set $\frac{D_{test}}{X_{test}} = 1$, which places a Gaussian kernel at every test point. This leaves only $\sigma$ to be selected by LCV. Let's go ahead and use 5-fold LCV in [pykliep](https://github.com/srome/pykliep) to select the optimal value of $\sigma$ for both the original features and the topics features.<jupyter_code>test_ratio = 1.0
sigmas = list(np.linspace(0.1, 1.0, 21))
cv = 5
kliep = DensityRatioEstimator(num_params=test_ratio, cv=cv, sigmas=sigmas, random_state=0)
kliep.fit(X_train.values, X_validation.values)
print('((num_params, sigma), max j score) = ', kliep._j_scores[0])
kliep_topics = DensityRatioEstimator(num_params=test_ratio, cv=cv, sigmas=sigmas, random_state=1)
kliep_topics.fit(X_train_topics.values, X_validation_topics.values)
print('((num_params, sigma), max j score) = ', kliep_topics._j_scores[0])<jupyter_output><empty_output><jupyter_text>Sample WeightsUsing the fitted models, we can now estimate the sample weights $\hat{w}(x)$, for both the original features and the topics features.<jupyter_code>weights = kliep.predict(X_train.values)
weights = pd.DataFrame({'full_name': X_train.index, 'weight': weights})
assert(len(weights) == len(train_features))
weights.head()
weights_topics = kliep_topics.predict(X_train_topics.values)
weights_topics = pd.DataFrame({'full_name': X_train_topics.index, 'weight': weights_topics})
assert(len(weights_topics) == len(train_features_topics))
weights_topics.head()<jupyter_output><empty_output><jupyter_text>It is interesting to visualize bar charts of the sample weights for the two sets of features. We identify the physicists by their index, rather than their name, for legibility of the figure. Furthermore, we distinguish between laureates and non-laureates and place the figures on the same scale to ease comparison.<jupyter_code>def plot_sample_weights(weights, weights_topics, train_target):

    """Plot bar charts of the sample weights of original features and topic features
    for the training data.

    Args:
        weights (array, shape (n_samples,)): Sample weights for training features.
        weights_topics (array, shape (n_samples,)): Sample weights for training topic features.
        train_target (pandas.DataFrame): Training target.

    Returns:
        tuple of matplotlib.axes.Axes: axes.
    """
    
    fig, (ax1, ax2) = plt.subplots(nrows=2, sharex=True, sharey=True, figsize=(10, 8))
    
    laureates_indices, non_laureate_indices = _get_laureate_non_laureate_indices(
        weights, train_target)
    _plot_axis(ax1, laureates_indices, non_laureate_indices, weights)
    ax1.set_title('Sample weights for original features')
    
    laureates_indices, non_laureate_indices = _get_laureate_non_laureate_indices(
        weights_topics, train_target)
    _plot_axis(ax2, laureates_indices, non_laureate_indices, weights_topics)
    ax2.set_xlabel('Physicist index')
    ax2.set_title('Sample weights for topic features')
    
    fig.tight_layout()
    return ax1, ax2

    
def _plot_axis(axes, laureates_indices, non_laureate_indices, weights):
    axes.bar(laureates_indices, weights.loc[laureates_indices, 'weight'], label='laureates',
            width=1.5, color='orange')
    axes.bar(non_laureate_indices, weights.loc[non_laureate_indices, 'weight'],
            label='non laureates', width=1.5, color='C0')
    axes.set_ylabel('Weight / importance')
    axes.set_ylim(0, 18)
    axes.set_xlim(0, len(weights))
    axes.legend()
    

def _get_laureate_non_laureate_indices(weights, target):
    physicists = pd.merge(weights, target, on='full_name')
    laureates_indices = physicists[physicists.physics_laureate == 'yes'].index
    non_laureate_indices = physicists[physicists.physics_laureate == 'no'].index
    return laureates_indices, non_laureate_indices
plot_sample_weights(weights, weights_topics, train_target);<jupyter_output><empty_output><jupyter_text>We can make the following observations about the two figures:1. The sample weights, on average, are higher for the original features than the topics features. This is not surprising as the covariate shift was much larger for the original features.2. For the original features, most of the weights are very low for the laureates. The laureates are not given much importance. This is worrying as we suspect that using these weights during learning may actually lead to a worse performing machine learning model! We will see what happens during model building. We suspect that this bevahior may be caused by overfitting during the LCV in KLIEP.3. For the topics features, the weights for laureates are far more comparable to those of non-laureates. They are certainly given importance that is of significance.4. In both charts, the highest weights are for non-laureate. This is most likely caused by the class imbalance in the data. There are far more non-laureates than laureates.   Persisting the Sample WeightsNow we have the sample weights for the original features and the topic features, let's persist them for future use in model building.<jupyter_code>weights.to_csv('../models/train-features-sample-weights.csv', index=False)
weights_topics.to_csv('../models/train-features-topics-sample-weights.csv', index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09395923225710459, 'fraction_numerical': 0.008406487129891385, 'mean_word_length': 5.47855421686747, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 27, 'https://': 8, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '404852', 'license': 'MIT', 'n_tokens_mistral': 4013, 'n_tokens_neox': 3862, 'n_words': 1444, 'path': 'nobel_physics_prizes/notebooks/5.2-importance-weighting.ipynb', 'repo_name': 'covuworie/nobel-physics-prizes', 'size': '13450'}