<jupyter_start><jupyter_text>Assume that we have two categorical features with 9 values each, so there are 81 possible feature combinations. The target y depends on the features x as y = mean(x) + eps(x), where mean(x) is some unknown fixed value and eps(x) is a normally distributed noise with mean 0 and variance var(x). We let mean(x) to be randomly generated and var(x) to have two values - 0.01 and 0.04. Also, there will be an "unknown" region without any training examples. To indicate this region on the figure, we set var(x) = 0 for these feature combinations.<jupyter_code># function for generating mean and variance
def gen_parameters(noise=0.01, seed=0):
    
    np.random.seed(seed)
    
    mean = np.random.rand(9, 9)
    
    figure = np.array([
        [1, 1, 1, 1, 1, 1, 1, 1, 1],
        [1, 1, 4, 1, 1, 1, 4, 1, 1],
        [1, 4, 0, 4, 1, 4, 0, 4, 1],
        [4, 0, 0, 0, 4, 0, 0, 0, 4],
        [1, 4, 0, 0, 0, 0, 0, 4, 1],
        [1, 1, 4, 0, 0, 0, 4, 1, 1],
        [1, 1, 1, 4, 0, 4, 1, 1, 1],
        [1, 1, 1, 1, 4, 1, 1, 1, 1],
        [1, 1, 1, 1, 1, 1, 1, 1, 1]
    ])

    var = figure*noise
    
    return mean, var
# generate parameters and plot variance
mean, var = gen_parameters()
print("True data uncertainty (variance). The white region will not contain any traing data.")
sns.heatmap(var, cmap="Reds", vmax=0.05)
plt.show()
# function for generating train and validation sets
def generate_training_data(n_samples, mean, var, num_cat=9, seed=0):
    np.random.seed(seed)
    train = []
    target = []
    val = []
    val_target = []
    for i in range(num_cat):
        for j in range(num_cat):
            if var[i, j] == 0:
                continue
            for _ in range(n_samples):
                train.append([i, j])
                val.append([i, j])
                target.append(np.random.normal(mean[i, j], math.sqrt(var[i, j])))
                val_target.append(np.random.normal(mean[i, j], math.sqrt(var[i, j])))
                
    train = np.asarray(train)
    target = np.asarray(target)

    return train, target, val, val_target<jupyter_output><empty_output><jupyter_text>Now we generate the dataset using the obtained distribution. In each cell, except for the white ones, we have several training examples. The test set consists of all possible feature combinations.<jupyter_code># generate train and validation datasets
train, target, val, val_target = generate_training_data(1000, mean, var)

train_pool = Pool(train, target, cat_features = [0, 1])
val_pool = Pool(val, val_target, cat_features = [0, 1])

# generate test, consisting of all possible feature combinations
num_cat = 9
test = np.asarray([[i, j] for i in range(num_cat) for j in range(num_cat)])<jupyter_output><empty_output><jupyter_text>If we optimize the RMSE loss, then we predict mean(x). If we also want to estimate var(x) (data uncertainty), then we have to use a probabilistic regression model that predicts mean and variance. For this purpose, we can use the loss function called RMSEWithUncertainty.<jupyter_code># predict mean value and data uncertainty

model = CatBoostRegressor(iterations=1000, learning_rate=0.2, loss_function='RMSEWithUncertainty',
                          verbose=False, random_seed=0)

model.fit(train_pool, eval_set=val_pool)
print("best iteration =", model.get_best_iteration())
preds = model.predict(test)

mean_preds = preds[:, 0] #the first prediction is the estimated mean value
var_preds = preds[:, 1] #the second prediction is the estimated variance<jupyter_output>best iteration = 978<jupyter_text>Let us see how well we can predict the mean values with this loss. Below we compare the true mean values with the estimated ones. For easier visual comparison, we use a mask to hide all feature combinations for which we do not have any training data.<jupyter_code># define mask to hide elements inside the heart:
mask = np.array([
        [1, 1, 1, 1, 1, 1, 1, 1, 1],
        [1, 1, 1, 1, 1, 1, 1, 1, 1],
        [1, 1, 0, 1, 1, 1, 0, 1, 1],
        [1, 0, 0, 0, 1, 0, 0, 0, 1],
        [1, 1, 0, 0, 0, 0, 0, 1, 1],
        [1, 1, 1, 0, 0, 0, 1, 1, 1],
        [1, 1, 1, 1, 0, 1, 1, 1, 1],
        [1, 1, 1, 1, 1, 1, 1, 1, 1],
        [1, 1, 1, 1, 1, 1, 1, 1, 1]
    ])
print("True mean values (with mask):")
sns.heatmap(mean*mask, cmap="Reds")
plt.show()

print("Predicted mean values (with mask):")
sns.heatmap(mean_preds.reshape([num_cat, num_cat])*mask, cmap="Reds")
plt.show()<jupyter_output>True mean values (with mask):<jupyter_text>We see that CatBoost captures the mean values well. Now let us plot the predicted data uncertainty (variance). We can see that CatBoost is able to predict the variance outside the heart and on its border. Inside the heart we have no training data, so anything can be predicted there.<jupyter_code>print("Predicted data uncertainty:")
sns.heatmap(var_preds.reshape([num_cat, num_cat]), cmap="Reds", vmin=0, vmax=0.05)
plt.show()<jupyter_output>Predicted data uncertainty:<jupyter_text>Ok, we know how to estimate noise in the data. But how to measure knowledge uncertainty coming from the lack of training data in a particular region? What to do if we want to detect outliers?  Estimating knowledge uncertainty requires an ensemble of models. If all the models give similar predictions, then there is low knowledge uncertainty, and if they give diverse predictions (strongly disagree with each other), then there is high knowledge uncertainty.Let’s generate an ensemble:<jupyter_code>def ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2):
    ens_preds = []
    for seed in range(num_samples):
        model = CatBoostRegressor(random_seed=seed, iterations=iters, learning_rate=lr, 
                                  loss_function='RMSEWithUncertainty', posterior_sampling=True,
                                  verbose=False)
        model.fit(train_pool, eval_set=val_pool)
        ens_preds.append(model.predict(test))
    return ens_preds<jupyter_output><empty_output><jupyter_text>Note that the models are generated using the option posterior_sampling, since this allows the obtained (random) predictions to be nicely distributed (with good theoretical properties).For regression, knowledge uncertainty can be obtained by measuring the variance of the mean across multiple models. Below we also plot the average estimated data uncertainty (average predicted variance for the modes).<jupyter_code>ens_preds = ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2)
ens_preds = np.asarray(ens_preds)
data = np.mean(ens_preds, axis=0)[:, 1] # average estimated data uncertainty
knowledge = np.var(ens_preds, axis=0)[:, 0] # estimated knowledge uncertainty

print("Knowledge uncertainty via ensemble:")            
sns.heatmap(knowledge.reshape([num_cat, num_cat]), cmap="Reds")
plt.show()

print("Average predicted data uncertainty:")            
sns.heatmap(data.reshape([num_cat, num_cat]), cmap="Reds", vmin=0, vmax=0.05)
plt.show()<jupyter_output>Knowledge uncertainty via ensemble:<jupyter_text>The models perfectly detected knowledge uncertainty inside the heart (we can see no traces of the original heart border on the figure). This illustrates how by estimating knowledge uncertainty, we can detect anomalous inputs.In practice, training an ensemble of several CatBoost models can be too expensive. If you need a cheaper method - you can use a virtual ensemble, obtained from a single trained model:<jupyter_code>def virt_ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2, seed=0):
    ens_preds = []
    model = CatBoostRegressor(iterations=iters, learning_rate=lr, loss_function='RMSEWithUncertainty', 
                              posterior_sampling=True, verbose=False, random_seed=seed)
    model.fit(train_pool, eval_set=val_pool)
    ens_preds = model.virtual_ensembles_predict(test, prediction_type='VirtEnsembles', virtual_ensembles_count=num_samples)
    return ens_preds<jupyter_output><empty_output><jupyter_text>As before, knowledge uncertainty can be obtained by measuring the variance of the mean across the obtained models. Similarly, we can also compute the average estimated data uncertainty.<jupyter_code>ens_preds = virt_ensemble(train_pool, val_pool, num_samples=10, iters=1000, lr=0.2)

ens_preds = np.asarray(ens_preds)
ens_preds = ens_preds

data = np.mean(ens_preds, axis = 1)[:,1]
knowledge = np.var(ens_preds, axis = 1)[:,0]

print("Knowledge uncertainty via virtual ensemble:")
sns.heatmap(knowledge.reshape([num_cat,num_cat]), cmap="Reds")
plt.show()

print("Average predicted data uncertainty:")
sns.heatmap(data.reshape([num_cat,num_cat]), cmap="Reds", vmax = 0.05)
plt.show()<jupyter_output>Knowledge uncertainty via virtual ensemble:<jupyter_text>We can also use prediction_type='TotalUncertainty' and get the same result easier. For this prediction type, CatBoost calculates and returns the following statistics using the virtual ensemble: [Mean Predictions, Knowledge Uncertainty, Data Uncertainty]<jupyter_code>model = CatBoostRegressor(iterations=1000, learning_rate=0.2, loss_function='RMSEWithUncertainty', 
                          posterior_sampling=True, verbose=False, random_seed=0)
model.fit(train_pool, eval_set=val_pool)
preds = model.virtual_ensembles_predict(test, prediction_type='TotalUncertainty', virtual_ensembles_count=10)
mean_preds = preds[:,0] # mean values predicted by a virtual ensemble

knowledge = preds[:,1] # knowledge uncertainty
print("Knowledge uncertainty via virtual ensemble:")
sns.heatmap(knowledge.reshape([num_cat,num_cat]), cmap="Reds")
plt.show()

data = preds[:,2] # estimated data uncertainty
print("Average predicted data uncertainty:")
sns.heatmap(data.reshape([num_cat,num_cat]), cmap="Reds", vmin = 0, vmax = 0.05)
plt.show()<jupyter_output>Knowledge uncertainty via virtual ensemble:
{'fraction_non_alphanumeric': 0.10691056910569105, 'fraction_numerical': 0.027845528455284554, 'mean_word_length': 4.259754142169963, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '185759', 'license': 'Apache-2.0', 'n_tokens_mistral': 3354, 'n_tokens_neox': 3062, 'n_words': 1080, 'path': 'catboost/tutorials/uncertainty/uncertainty_regression.ipynb', 'repo_name': 'jochenater/catboost', 'size': '9840'}