!pip install ucimlrepo

from ucimlrepo import fetch_ucirepo 
  
# Fetch dataset 
dry_bean_dataset = fetch_ucirepo(id=602) 
  
# Load data (as pandas dataframes) 
X = dry_bean_dataset.data.features 
y = dry_bean_dataset.data.targets 

X

y

from sklearn.preprocessing import LabelEncoder

y = np.ravel(y)
# Convert string labels to integer labels

label_encoder = LabelEncoder()
y_encoded = label_encoder.fit_transform(y)
y_encoded

# Print the mapping of encoded integer labels to original string labels

print("Mapping of Encoded Labels to Original Labels:")
for encoded_label, original_label in zip(label_encoder.classes_, label_encoder.transform(label_encoder.classes_)):
    print(f"{encoded_label} -> {original_label}")

beans_names = list(label_encoder.classes_)
beans_names

# Split data into training, calibration and test sets

from sklearn.model_selection import train_test_split

X_train_cal, X_test, y_train_cal, y_test = train_test_split(X, y_encoded, test_size=100, random_state=11)
X_train, X_cal, y_train, y_cal = train_test_split(X_train_cal, y_train_cal, test_size=0.25, random_state=11)

print("Training data size:", X_train.shape, y_train.shape)
print("Calibration data size:", X_cal.shape, y_cal.shape)
print("Test data size:", X_test.shape, y_test.shape)

from sklearn.naive_bayes import GaussianNB
from sklearn import metrics

model = GaussianNB().fit(X_train, y_train)

# Check accuracy on the calibration set

y_pred = model.predict(X_cal) 
accuracy = metrics.accuracy_score(y_cal, y_pred)
print('Accuracy:', accuracy)

# Get predicted probabilities for each class
predictions = model.predict_proba(X_cal)  

idx = 23
# Extract the predicted probability for each class for a data point
probabilities_for_data_point = predictions[idx]
print('Probabilities for a sample data point:\n',probabilities_for_data_point)

# Get the true class of the selected data point
true_class = y_encoded[idx]
print('\nTrue Class of the sample data:', true_class)

plt.bar(x=range(len(probabilities_for_data_point)), height=probabilities_for_data_point)
plt.xlabel('Class')
plt.ylabel('Probability Estimate')
plt.title('Probability Estimate for a sample data point')
plt.show()

# Calculate adjusted quantile

alpha = 0.15 # means (1-alpha) coverage
print('alpha:', alpha)

n = X_cal.shape[0]
print('Calibration set size:', n)

# Get predicted probabilities for each class
y_cal_predictions = model.predict_proba(X_cal)  
print('Predict Probability shape:', predictions.shape)

# Get the probabaility of the true class
prob_true_class = y_cal_predictions[np.arange(len(y_cal)),y_cal] # The second index taken from y_cal is 1 only for the true class
print('Prob true class shape:', prob_true_class.shape)

y_cal_scores = 1-prob_true_class # Larger conformal scores encode worse aggrement between X and y

q_level = np.ceil((n+1)*(1-alpha))/n # finite-size correction
print("q-level", q_level)

qhat = np.quantile(y_cal_scores, q_level)
print("qhat", qhat)

# Flatten the array of conformal scores
flat_scores = y_cal_scores.flatten()

# Plot histogram of frequency of conformal scores
plt.hist(flat_scores, bins=10, color='skyblue', edgecolor='black')

# Add vertical dotted line denoting qhat
plt.axvline(x=qhat, color='red', linestyle='--', label=f'qhat = {qhat:.3f}')
plt.xlabel('Conformal Score')
plt.ylabel('Frequency')
plt.title('Histogram of Conformal Scores for the Calibration Set')
plt.legend()
plt.show()

id = 2
y_test_predictions = model.predict_proba(X_test)
print('X_test predicted scores:\n', y_test_predictions[id])

predictions_set = (1-y_test_predictions <= qhat)
print('\nPrediction set:\n', predictions_set[id])
print('\nTrue class:', y_test[id])

accepted_test_set = np.zeros(len(predictions_set[id]))
for i in range(len(predictions_set[id])):
    if predictions_set[id][i]:
        accepted_test_set[i] = y_test_predictions[id][i]


# Create a list to represent the heights of the bars of the true class
true_bar_heights = [0] * len(predictions_set[id])
true_bar_heights[y_test[id]] = 1

x=list(range(len(y_test_predictions[id])))

bar_width = 0.35

plt.figure(figsize=(10, 6))

plt.bar(x, y_test_predictions[id], width=bar_width, color='lightblue', label='Predicted Probabilities')
plt.bar(x, accepted_test_set, width=bar_width, color='lightblue', edgecolor='black', hatch='\\', label='Conformal Predicted Classes')
plt.axhline(y=1-qhat, color='red', linestyle='--', label=f'1-qhat = {1-qhat:.3f}')

plt.bar([i + bar_width for i in x], true_bar_heights, width=bar_width, color='red', label='True Class')

xticks_positions = [i + bar_width / 2 for i in x]
plt.xticks(xticks_positions, beans_names, rotation=90, fontsize=10)

plt.xlabel('Class')
plt.ylabel('Probability Estimates')
plt.title('Probability estimates for a test data point')
plt.legend()
plt.show()