Conformal prediction - wstęp i przypomnienie¶
Conformal prediction to metoda, która umożliwia konstruowanie prognoz z predefinowaną gwarancją dokładności. Metoda może być stosowana jako uzupełnienie dowolnego tradycyjnego algorytmu predykcyjnego zarówno do klasyfikacji, jak i regresji.
- Rozszerzymy prosty model SVM, który będzie zwracał zbiór (przedział).
- Zbudujemy klasyfikator (conformal predictor) krok po kroku.
- Potwierdzimy, że skonstruowany model daje takie same wyniki jak funkcje z biblioteki
nonconformist.
Idea regresji (conformal regression) jest taka sama jak w przypadku klasyfikacji (konformal classification).
Konfig¶
Biblioteki
import pandas as pd
import numpy as np
from sklearn.svm import SVC
import matplotlib.pyplot as plt
#plt.style.use('dark_background')
Dataset¶
Jako danych, na których będziemy operować użyjemy syntetycznego zbioru danych z 4 klasami i 2000 punktami w każdej klasie. Punkty klasy będą rozmieszczone wokół odpowiedniego środka klasy na płaszczyźnie 2D.
np.random.seed(20)
# centers of 4 classes
centers = [
[ 1, 0],
[ 0, 1],
[-1, 0],
[ 0, -1],
]
# colors to plot members of different classes
color_arr = ['r', 'b', 'g', 'y']
# number of instances per class
n_points = 2000
# these arrays will store information about datapoints
data_x = []
data_y = []
data_class = []
data_color = []
# standard deviation to generate the class instances
sigma = 0.5
# data generation
for class_val in range(0, len(centers)):
x, y = centers[class_val]
data_class.extend([class_val for j in range(0, n_points)])
data_color.extend([color_arr[class_val] for j in range(0, n_points)])
data_x.extend(np.random.normal(x, sigma, size=n_points))
data_y.extend(np.random.normal(y, sigma, size=n_points))
# putting everything into a dataframe
data_df = pd.DataFrame({
'x': data_x,
'y': data_y,
'class': data_class,
'color': data_color,
})
# plotting the dataset
data_df.plot(
kind='scatter',
x='x',
y='y',
c=data_df['color'],
s=1,
grid=True,
figsize=(10,10),
)
# showing the centers in orange
plt.scatter(np.array(centers).T[0], np.array(centers).T[1], s=100, c='red', label='Centers')
plt.legend()
plt.show()
Prosty klasyfikator¶
Zbudujemy tradycyjny (punktowy) klasyfikator przy użyciu algorytmu SVM (Support Vector Machine).
Najpierw dzielimy dane na zbiór treningowy (90\%) i testowy (10\%).
# fraction of the dataset to be used for testing
test_frac = 0.1
np.random.seed(2)
# performing random permutation of the dataset
idx = np.random.permutation(len(data_df))
# constucting training and test datasets
test_size = int(len(data_df) * 0.1)
train_size = len(data_df) - test_size
idx_train = idx[:train_size]
idx_test = idx[train_size:]
print('Size of training set: {}'.format(len(idx_train)))
print('Size of test set: {}'.format(len(idx_test)))
Teraz jesteśmy gotowi do zbudowania modelu SVM i przetestowania go.
model = SVC(probability=True)
model.fit(data_df[['x', 'y']].values[idx_train, :], data_df['class'].values[idx_train])
prediction = model.predict(data_df[['x', 'y']].values[idx_test, :])
acc = sum(prediction[:] == data_df['class'].values[idx_test][:]) / len(idx_test)
print('Accuracy: {}%'.format(acc * 100))
Osiągnęliśmy dokładność około $65\%$.
Analizując rozkład punktów danych (pokazany poniżej dla odniesienia), możemy stwierdzić, że prawdopodobnie żaden inny algorytm klasyfikacji nie poprawi znacząco tego wyniku. Wiele punktów danych nakłada się w naszej przestrzeni 2D i nie ma wystarczających informacji, aby je rozróżnić.
x_grid = np.linspace(-4,4,200)
y_grid = np.linspace(-4,4,200)
xx, yy = np.meshgrid(x_grid, y_grid)
r1, r2 = xx.flatten(), yy.flatten()
r1, r2 = r1.reshape((len(r1), 1)), r2.reshape((len(r2), 1))
grid = np.hstack((r1,r2))
yhat = model.predict(grid)
zz = yhat.reshape(xx.shape)
data_df.plot(kind='scatter', x='x', y='y', c=data_df['color'], s=1, grid=True, figsize=(10,10),)
# showing the centers in orange
plt.scatter(np.array(centers).T[0], np.array(centers).T[1], s=70, c='red', label='Centers')
plt.contourf(xx, yy, zz, levels = [-1,0,1,2,3], colors=color_arr, alpha=.2)
plt.legend()
plt.show()
Powiedzmy, że potrzebujemy większej dokładności w naszej aplikacji. Jesteśmy gotowi poświęcić precyzję naszych prognoz, aby osiągnąć wymaganą dokładność. Na przykład jesteśmy gotowi zaakceptować zdanie "$a$ należy do klasy 1 lub 3 z prawdopodobieństwem $0,95$" zamiast "$a$ należy do klasy 3 z prawdopodobieństwem $0,65$".
Przez precyzję rozumiemy tutaj przewidywaną liczbę etykiet klas przypadających na obserwację. Tradycyjny klasyfikator (punktowy) przewiduje tylko jedną etykietę klasy, jest ona możliwie najdokładniejsza. Tutaj jesteśmy gotowi poświęcić tę precyzję. Jesteśmy gotowi dopuścić 0, 1 lub więcej etykiet klas na punkt danych, aby zagwarantować wymaganą dokładność. W rezultacie będziemy tworzyć prognozy zbioru (przedziału).
Poniżej kilka przykładów możliwych zastosowań:
- medycyna
- ubezpieczenie
- bankowość
- ...
Czy jest więc sposób na osiągnięcie, np dokładności $95\%$ przy dopuszczeniu większej ilości klas?
Jak uzyskać dokładność $95\%$ używając prostego modelu?¶
Jeden z możliwych sposobów osiągnięcia wymaganej dokładności może opierać się na analizie prawdopodobieństw wygenerowanych przez podstawowy model klasyfikacji. Użyjemy do tego stworzonego już modelu SVM.
prediction_proba = model.predict_proba(data_df[['x', 'y']].values[idx_test, :])
print('Probabilities of 4 possible class labels generated by our SVM classifier \n\
for the first 5 datapoints in the test set:\n')
print(' class_1 class_2 class_3 class_4')
for i in range(0, 5):
print('{}: {}'.format(i + 1, prediction_proba[i]))
Predyktory punktu standardowego wyprowadzają etykietę klasy skojarzoną z najwyższą wartością prawdopodobieństwa. Wyniki będą wyglądać tak (True wskazuje przewidywaną etykietę klasy):
prediction_point = np.zeros_like(prediction_proba)
prediction_point[np.arange(len(prediction_proba)), prediction_proba.argmax(1)] = 1
np.array(prediction_point, dtype=bool)
for i in range(0, 5):
print('{}: {}'.format(i + 1, np.array(prediction_point[i], dtype=bool)))
Jeśli musimy osiągnąć dokładność $95\%$, możemy na wyjściu wypluć wszystkie etykiety klas z powiązanym prawdopodobieństwem większym niż $\geq 1 - 0,95 = 0,05$. W ten sposób wszystkie odrzucone etykiety klas nie powinny powodować więcej niż $5\%$ błędów. Będziemy odnosić się do tej wartości jako dozwolony poziom błędów, istotność lub $\epsilon$.
significance = 0.05
prediction_region_naive = prediction_proba > significance
for i in range(0, 5):
print('{}: {}'.format(i + 1, prediction_region_naive[i]))
Zgodnie z oczekiwaniami, teraz przewidujemy się więcej niż 1 etykietę klasy dla danej obserwacji, otrzymaliśmy prognozę w postaci zbioru.
Oszacujmy teraz, jak dobre jest przewidywanie zbudowanego zbioru.
- Najpierw chcemy potwierdzić, że osiągnęliśmy pożądaną dokładność 95\%$.
- Po drugie, byłoby interesujące ocenić, jak daleko jesteśmy teraz od precyzyjnego predyktora punktowego. W szczególności możemy obliczyć wartości następujących metryk wydajności:
- $oneC$ - ułamek przewidywań zbioru jednoelementowego. W przypadku predyktora punktowego wszystkie predykcje są singletonami. Oczywiście chcielibyśmy, aby ta wartość była bliska 1; czyli preferujemy większe wartości.
- $avgC$ - średni rozmiar zbioru predykcyjnego. Dla predyktora punktowego mamy $avgC=1$. W przypadku predyktora regionu skonstruowanego powyżej wiele punktów danych ma więcej niż 1 przewidywaną etykietę. Dlatego preferowane są mniejsze wartości $avgC$.
Metryki do oceny klasyfikatorów:
# Accuracy for a region prediction:
# A region prediction produces an error if the resulting prediction set does not
# contain the true value
def get_accuracy(prediction, real_class):
correct = 0
N = len(prediction)
for i in range(0, N):
if prediction[i][real_class[i]]:
correct += 1
return correct / N
# calculating metrics: oneC & avgC
def get_oneC_avgC(prediction):
arr = np.array(prediction)
oneC = 0
avgC = 0
for i in range(0, len(arr)):
# number of predicted lables
num_predicted = arr[i].sum()
avgC += num_predicted
# is it a singleton?
if num_predicted == 1:
oneC += 1
pass
oneC /= len(arr)
avgC /= len(arr)
return oneC, avgC
Oszacujmy teraz skuteczność zbioru predykcyjnego skonstruowanego powyżej:
acc = get_accuracy(prediction_region_naive, data_df['class'].values[idx_test])
oneC, avgC = get_oneC_avgC(prediction_region_naive)
print('Accuracy of the region predictor: {}%'.format(acc * 100))
print('Efficiency:\n\t oneC = {}\n\t avgC = {}'.format(oneC, avgC))
Osiągnęliśmy jeszcze lepszą dokładność niż wymagana: $96,6 \% $.
Być może, wykluczając niektóre z przewidywanych etykiet, moglibyśmy poprawić skuteczność przewidywania naszego regionu. Oznacza to, że możemy osiągnąć lepsze wartości $oneC$ i $avgC$. Właśnie w tym mogą nam pomóc conformal prediction.
Conformal classifier krok po kroku¶
Podzielmy zbiór danych na 3 części: treningowy, testowy i kalibracyjny.
Zbiory uczące i testowe będą wykorzystywane w taki sam sposób jak dotychczas: do trenowania modelu predykcyjnego i oceny jego efektywności. Użyjemy części dotyczącej kalibracji, aby oszacować, kiedy konkretna etykieta powinna zostać uwzględniona w ostatecznym przewidywaniu regionu.
Tak jak poprzednio, zbiór testowy będzie zawierał $10\%$ dostępnych danych. Reszta zostanie podzielona na treningowy, $80\%$ i kalibracyjny, $20\%$.
# fraction of the dataset to be used for testing
test_frac = 0.1
# fraction of the remaining dataset to be used for calibration
calib_frac = 0.2
np.random.seed(2)
# random permulation
idx = np.random.permutation(len(data_df))
# constucting test, training and calibration datasets
test_size = int(len(data_df) * test_frac)
calib_size = int(len(data_df) * (1 - test_frac) * calib_frac)
train_size = len(data_df) - test_size - calib_size
idx_train, idx_cal, idx_test = idx[:train_size], idx[train_size:train_size + calib_size], idx[train_size + calib_size:]
print('Test size: {}'.format(test_size))
print('Calibration size: {}'.format(calib_size))
print('Train size: {}'.format(train_size))
Trenujemy model SVM na zbiorze treningowym:
model.fit(data_df[['x', 'y']].values[idx_train, :], data_df['class'].values[idx_train])
Intuition: measuring conformity of predictions on the calibration dataset¶
Przeanalizujmy przewidywania dotyczące zestawu kalibracyjnego. Zobaczmy, które wartości prawdopodobieństwa są skojarzone z prawdziwymi etykietami klas i jak są dystrybuowane. Oznacza to, że przeanalizujemy, jak konformalny jest każdy punkt danych z zestawu kalibracyjnego (funkcja zgodności).
# generating predictions for the calibration dataset
predictions_cal = model.predict_proba(X=data_df[['x', 'y']].values[idx_cal, :])
# extracting the predicted probability for the correct class labesl
calib_conformal_vals = []
for i in range(0, len(idx_cal)):
# correct class label
cl = data_df['class'].values[idx_cal][i]
# associated probability
calib_conformal_vals.append(predictions_cal[i][cl])
calib_conformal_vals[:5]
plt.hist(calib_conformal_vals)
plt.show()
np.quantile(calib_conformal_vals, 0.05)
Posortujmy teraz otrzymane wartości i narysujmy "dystrybuantę empiryczną"
calib_conformal_vals = np.sort(calib_conformal_vals)
plt.plot(calib_conformal_vals)
plt.grid(True)
plt.ylabel('Conformity value')
plt.title('Distribution of conformity values')
Przeanalizujmy teraz, gdzie w tym rozkładzie znajdują się przewidywane wartości prawdopodobieństwa dla testowego zbioru danych i jaki jest związany z nim empiryczny poziom błędu nieuwzględnienia odpowiedniej etykiety klasy zbiorze predykcyjnym.
predictions_test = model.predict_proba(X=data_df[['x', 'y']].values[idx_test, :])
Wykonamy analizę dla pierwszej obserwacje testowej. Oto wygenerowane prawdopodobieństwa związane z 4 możliwymi etykietami klas:
predictions_test[0]
Czy 1 klasa powinna być w zbiorze predykcyjnym?¶
print('Distribution of probabilities for the first datapoint in the test set:')
predictions_test[0]
Przeanalizujmy prawdopodobieństwo związane z klasą 1 w naszym rozkładzie empiryczny,:
plt.plot(calib_conformal_vals)
plt.axhline(y=predictions_test[0][0], color='r', linestyle='-')
plt.legend(['calib dist', 'class_1'])
plt.grid(True)
plt.ylabel('Conformity value')
plt.title('Distribution of conformity values')
Jeżeli do rozkładu empirycznego uwzględnimy nowo obliczoną wartość, to po lewej stronie znajduje się następujący ułamek punktów (w tym nowo dodany):
val = sum(calib_conformal_vals <= predictions_test[0][0]) + 1
print('Number of points on the left hand side (<=): {}'.format(val))
print('Fraction of points on the left hand side: {}%'.format(np.round(val / (len(calib_conformal_vals) + 1) * 100)))
And the following fraction of the points are situated on the right-hand side:
val = (sum(calib_conformal_vals > predictions_test[0][0]) + 1)
print('Number of points on the right hand side (>): {}'.format(val))
print('Fraction of points on the right hand side: {}%'.format(np.round(val / (len(calib_conformal_vals) + 1) * 100)))
Przypomnijmy, że ten rozkład został skonstruowany dla wartości prawdopodobieństwa prawdziwych etykiet klas. Oznacza to, że ponad $43\%$ zaobserwowanych prawdopodobieństw właściwych klas miało niższe wartości niż analizowane. Stąd możemy wywnioskować, co następuje. Jeśli nie uwzględnimy powiązanej etykiety klasy (w tym przypadku 1), to z obserwowanego rozkładu prawdopodobieństwo popełnienia błędu przy tej akcji wynosi $43\%$.
Jest wyższa niż dopuszczalna stopa błędów $\epsilon = 5\%$. Oznacza to, że należy uwzględnić tę etykietę klasy.
Czy 2 klasa powinna być w zbiorze predykcyjnym?¶
Zróbyme tę samą analizę dla kklasy 2
print('Distribution of probabilities for the first datapoint in the test set:')
predictions_test[0]
plt.plot(calib_conformal_vals)
plt.axhline(y=predictions_test[0][1], color='r', linestyle='-')
plt.grid(True)
plt.legend(['calib dist', 'class_2'])
plt.ylabel('Conformity value')
plt.title('Distribution of conformity values')
Jeżeli do rozkładu empirycznego uwzględnimy nowo obliczoną wartość, to po lewej stronie znajduje się następujący odsetek obserwacji (wraz z nowo dodanym):
val = sum(calib_conformal_vals <= predictions_test[0][1]) + 1
print('Number of points on the left hand side (<=): {}'.format(val))
print('Fraction of points on the left hand side: {}%'.format(np.round(val / (len(calib_conformal_vals) + 1) * 100, 2)))
A następująca część obserwacji znajduje się po prawej stronie:
val = (sum(calib_conformal_vals > predictions_test[0][1]) + 1)
print('Number of points on the right hand side (>): {}'.format(val))
print('Fraction of points on the right hand side: {}%'.format(np.round(val / (len(calib_conformal_vals) + 1) * 100, 2)))
Using the same reasoning, we can conclude that if this class label is not included in the region prediction, the probability of making a mistake is $3.33\%$. It is less than our allowed error rate $\epsilon = 5\%$. Thereby, this class label should not be included.
Generalizacja na cały zbiór testowy¶
Użyjemy tego samego rozumowania, aby zdecydować, czy każda etykieta klasy powinna zostać uwzględniona w końcowej prognozie regionu dla każdego punktu danych z naszego testowego zbioru danych. Korzystając z empirycznego rozkładu wartości zgodności obliczonych na zbiorze danych kalibracyjnych, obliczymy prawdopodobieństwo popełnienia błędu przez nieuwzględnienie określonej etykiety klasy. Ta wartość jest określana jako p-wartość. Jeśli $p\_value < \epsilon$, odpowiednia etykieta klasy jest odrzucana. W przeciwnym razie jest uwzględniany w prognozie regionu.
p_vals_inverse = np.zeros((len(idx_test), len(centers)))
for i in range(0, len(idx_test)):
for j in range(0, len(centers)):
p_vals_inverse[i,j] = (sum(calib_conformal_vals <= predictions_test[i,j]) + 1) / (len(calib_conformal_vals) + 1)
p_wartości dla pierwszej obserwacji:
p_vals_inverse[0]
Zbiór predykcyjny:
p_vals_inverse[0] > significance
Wyniki dla całego zbioru testowego:
prediction = p_vals_inverse > significance
prediction
Oszacujmy teraz dokładność wynikowego zbioru predykcyjnego:
print('Conformal classification based on conformity values:')
acc = get_accuracy(prediction, data_df['class'].values[idx_test])
oneC, avgC = get_oneC_avgC(prediction)
print('Accuracy of region predictor: {}%'.format(acc * 100))
print('Efficience:\n\t oneC = {}\n\t avgC = {}'.format(oneC, avgC))
Poprzednie podejście (nie korzystamy z technik conformal):
acc = get_accuracy(prediction_region_naive, data_df['class'].values[idx_test])
oneC, avgC = get_oneC_avgC(prediction_region_naive)
print('Accuracy of the region predictor: {}%'.format(acc * 100))
print('Efficiency:\n\t oneC = {}\n\t avgC = {}'.format(oneC, avgC))
plt.hist(np.sum(prediction, axis=1))
plt.show()
Adaptacyjność¶
Własność conformal prediction (marginal coverage):
$$ P [Y_{test} \in 𝓣 (X_{test})] \ge 1 - \alpha $$Własność pokrycia warunkowego (conditional coverage):
$$ P [Y_{test} \in 𝓣 (X_{test}) | X_{test}] \ge 1 - \alpha $$| Poprzednia metoda | Nowa metoda |
|---|---|
| najmniejsza możliwa średnia moc zbiorów C | zwykle większy średna moc zbiorów C |
| nieadaptacyjna | adaptacyjna |
| używa tylko wyjścia dla poprawnej klasy | używa informacji ze wszystkich wyjść klas |
# generating predictions for the calibration dataset
predictions_cal = model.predict_proba(X=data_df[['x', 'y']].values[idx_cal, :])
# extracting the predicted probability for the correct class labesl
calib_conformal_vals = []
for i in range(0, len(idx_cal)):
#for i in range(0, 10):
# correct class label
cl = data_df['class'].values[idx_cal][i]
# associated probability
calib_conformal_vals.append(predictions_cal[i][cl])
threshold = predictions_cal[i][cl]
calib_conformal_vals.append(sum(predictions_cal[i][predictions_cal[i] >= threshold]))
plt.hist(calib_conformal_vals)
plt.show()
calib_conformal_vals = np.sort(calib_conformal_vals)
plt.plot(calib_conformal_vals)
plt.grid(True)
plt.ylabel('Conformity value')
plt.title('Distribution of conformity values')
q_hat = np.quantile(calib_conformal_vals, 1-significance)
q_hat
idx_sort = np.argsort(-predictions_test, axis=1)
predictions = np.zeros(prediction.shape, dtype = bool)
i=0
for row in idx_sort:
cumsum = 0
true_idx = []
for idx in row:
true_idx.append(idx)
cumsum += predictions_test[i][idx]
if(cumsum >= q_hat):
for t_idx in true_idx:
predictions[i][t_idx] = True
break
i+=1
predictions
acc = get_accuracy(predictions, data_df['class'].values[idx_test])
oneC, avgC = get_oneC_avgC(predictions)
print('Accuracy of region predictor: {}%'.format(acc * 100))
print('Efficience:\n\t oneC = {}\n\t avgC = {}'.format(oneC, avgC))
plt.hist(np.sum(predictions, axis=1))
plt.show()
Jak mierzyć zdolności adaptacyjne metody?¶
Metryka SSC (size-stratified coverage):
$$ SSC = \min_{g \in \{1, \dots, G \}} \sum_{i \in I_g } \mathbb{1} \{ Y_{i}^{(val)} \in 𝓣 (X_{i}^{(val)}) \} $$Metody adaptacyjne zapewniają, że SSC $> 1-\alpha$
def get_ssc(prediction, real_class):
ssc = 1
for size in [1,2,3,4]:
idx_size = np.argwhere(np.sum(prediction, axis=1) == size).reshape(-1)
if(idx_size.shape[0] != 0):
acc = get_accuracy(prediction[idx_size], real_class[idx_size])
else:
acc=1
if (acc < ssc):
ssc = acc
print(f"Dla mocy zbioru = {size} mamy dokładność {acc}")
return(ssc)
print("SSC dla podstawowej metody conformal prediction:")
round(get_ssc(prediction, data_df['class'].values[idx_test]),4)
print("SSC dla podstawowej metody adaptacyjnej:")
round(get_ssc(predictions, data_df['class'].values[idx_test]),4)