mirror of
https://github.com/ada-dmitry/msmp-mephi_python.git
synced 2026-09-24 00:50:14 +00:00
69 KiB
69 KiB
In [38]:
# 1) Названия файлов с исходными данными и пустого Excel-файла для вывода программы
# Если файлы лежат не в одной директории с программой, необходимо указать полный путь к файлам
# В целом, если конечный файл не создать, то программа его создаст сама, но лучше явно ей указать его.
INPUT_PATH = 'data.xlsx'
OUTPUT_PATH = 'output.xlsx'
# 2) Количество признаков в методе с включением
CNT_ENTER = 6
# 3) Количество признаков в методе с исключением
CNT_REMOVE = 6
# 4) Обучающая выборка в формате словаря, где ключ - номер класса, значение - список номеров строк, принадлежащих этому классу
TRAIN_SAMPLE = {
1: [4,23,28,53,59,69],
2: [13,21,34,36,39,67,70,83],
3: [14,20,24,40,45,61,65,68,74,78],
4: [2,54],
5: [18,31],
6: [10,62],
7: [47,79]
}
In [56]:
import pandas as pd
import numpy as np
from scipy.stats import f
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from scipy.spatial.distance import mahalanobis
data = pd.read_excel(INPUT_PATH,
usecols=[i for i in range(0,10)])
data = data.iloc[:85]
# Стандартизация только числовых колонок (начиная со второй)
numeric_cols = data.columns[1:] # Все колонки кроме первой
data[numeric_cols] = (data[numeric_cols] - data[numeric_cols].mean(axis=0)) / data[numeric_cols].std(axis=0)
print('data:')
print(data.head())
FEATURES = [f'X{i}' for i in range(1,10)]
# Сохраняем первую колонку (названия субъектов) и признаки X1-X9 для вывода в Excel
data_to_excel = data[[data.columns[0]] + FEATURES].copy()
data:
Наименование X1 X2 X3 \
0 Алтайский край -0.567424 -0.126713 -1.023246
1 Амурская область 0.167298 -0.608339 -0.667051
2 Архангельская область без автономного округа 0.246382 0.884700 -0.982147
3 Астраханская область -0.560670 0.354912 -0.105358
4 Белгородская область -0.116102 -0.030388 0.922129
X4 X5 X6 X7 X8 X9
0 -0.917052 -0.432397 1.006771 -0.250080 -0.730154 0.163593
1 -0.769444 0.622505 0.216527 -0.640826 -0.087726 1.593888
2 -0.707543 1.135777 1.441228 -0.657606 -0.427835 -0.629321
3 -0.221865 1.041252 -0.715829 -0.525759 0.894809 -0.162722
4 0.778060 -0.533539 -0.390494 1.025235 -0.163306 -0.827550
In [40]:
def get_train_data(data, features, train_samples=None): #объявление функции принимает 3 параметра последние не обязательно вводить (3 это обучающая выборка) фичерс это Х загаловки а дата это таблица с данными
train_data = pd.DataFrame() # создаётся переменная трейн дата( новый дата фрейм дата фрейм это таблица )
for cls, samples in train_samples.items(): # type: ignore # заходим в цикл переменная cls будет последовательно принемать значения номеров классов то есть каждый цикл это следующий номер класса, sempls бедет принимать список номеров строк которые принадлежат этому классу точка айтемс это пара значений (ключ и значение) для словаря train samples ключ это номер класса из обучающей выборки а значение это номера строк этого класса из обучающей выборки
train_samps = data[features].loc[samples] #train_samps это новая таблица в которую мы записываем значения(Х) которое входят в обучающаю выборку (значение номер класса и так до конца)
train_samps["Class"] = cls # это доп колонка номер класса в которую вносится номер класса к которому принадлежит внесённый X
train_data = pd.concat([train_data, train_samps]) # сохраняем полоученные строчки в train_data а train_samps в котором они были до этого очищаем
train_data = train_data.astype({"Class": 'int32'}) # в train_data присваиваем итоговую таблицу с обучающей выборкой и показываем что в колонке класс данные типа интедгер
return train_data
train_data = get_train_data(data, FEATURES, TRAIN_SAMPLE) # type: ignore
print(train_data)
# Создаём колонку с классами для всех строк (NaN для строк не из обучающей выборки)
train_sample_column = pd.Series(index=data.index, dtype='Int32')
train_sample_column[train_data.index] = train_data.Class
data_to_excel['Train sample'] = train_sample_column
X1 X2 X3 X4 X5 X6 X7 \
4 -0.116102 -0.030388 0.922129 0.778060 -0.533539 -0.390494 1.025235
23 0.191144 -0.800989 0.470034 0.663782 -0.446576 -0.809216 0.423536
28 -0.191100 -0.367526 1.018027 -0.064734 -0.955121 -0.997007 0.629696
53 -0.808551 -0.560176 1.524921 0.468559 0.398480 -0.924936 0.001627
59 -0.133486 -0.271201 0.771431 0.049543 -1.113924 -0.471701 -0.075084
69 -0.450778 -0.271201 1.100226 0.240005 -0.400258 -0.893468 -0.012757
13 -0.372445 -0.271201 -1.338342 -0.731351 -0.250908 -0.329589 0.080735
21 -0.490466 -0.223038 -0.872548 -1.012282 0.369177 0.132275 -0.451446
34 -0.464369 -0.656501 -1.379442 -0.331381 -0.420108 -0.242799 -0.247683
36 -0.385119 -0.030388 -1.434241 -0.098065 0.190525 -0.204226 -0.489801
39 -0.567507 -0.367526 -0.269756 -0.088542 -0.283047 -0.404197 -0.060701
67 -0.304660 -0.319364 -0.584852 -0.807536 0.354053 -0.890931 0.078337
70 -0.365525 -0.223038 -0.351955 -0.745636 -0.107230 0.055636 -0.111042
83 -0.139989 0.065937 -0.762950 -0.669451 0.486389 0.402288 -0.041523
14 -0.315123 0.017774 -0.242356 -0.431374 -0.080763 0.619009 -0.650414
20 -0.443316 -0.126713 -0.310856 -0.550412 -0.622393 -0.317915 -0.290833
24 0.014092 -0.800989 -0.091658 0.268574 -0.038226 0.004882 -0.700756
40 -0.183554 -0.415689 -0.584852 0.397136 -0.064694 0.397720 -0.377132
45 -0.417719 -0.656501 -0.119058 0.563790 -0.562842 0.910338 0.049571
61 -0.306035 -0.800989 0.031640 0.782821 0.034558 0.468269 0.114296
65 0.140534 -0.608339 -0.612252 -0.388520 -0.468316 0.399243 -0.417885
68 -0.713376 0.065937 0.196038 -0.002834 -0.522196 0.016048 -0.055906
74 -0.462868 -0.560176 -0.105358 0.025735 0.344601 0.278955 -0.082276
78 -0.431018 -0.608339 0.168638 0.597121 -0.678162 0.488063 -0.132617
2 0.246382 0.884700 -0.982147 -0.707543 1.135777 1.441228 -0.657606
54 0.791628 0.981025 0.648132 -0.364712 0.940110 1.300639 -0.719694
18 1.542984 -0.752826 0.374136 -1.755084 1.200054 0.153592 -0.718975
31 1.073236 -0.126713 0.538534 -1.255121 0.591312 0.403811 -0.671989
10 2.795688 -0.897314 2.127713 2.358894 2.356098 -0.758970 5.316241
62 0.916945 -0.993639 0.291937 -0.240912 3.882681 1.425494 5.364185
47 -0.431727 3.967103 -0.201257 1.549430 -0.964574 -0.246352 0.413947
79 -0.588726 3.389152 0.127539 2.635063 -1.607346 -2.105480 0.833458
X8 X9 Class
4 -0.163306 -0.827550 1
23 -0.201096 0.154444 1
28 0.554701 0.194090 1
53 -0.465625 -0.800103 1
59 -0.314465 -0.998331 1
69 -0.503415 0.401467 1
13 0.781440 0.471610 2
21 -0.314465 1.185232 2
34 -0.238886 0.496007 2
36 0.932599 1.532895 2
39 0.592491 0.755229 2
67 0.327962 0.230686 2
70 -0.427835 0.599696 2
83 0.365752 0.590547 2
14 0.403541 -0.059033 3
20 -0.352255 -0.080381 3
24 0.025643 0.105649 3
40 -0.465625 -0.559179 3
45 -0.012147 -0.357901 3
61 -0.654574 -0.116977 3
65 -0.730154 -1.041027 3
68 -1.070262 -0.424993 3
74 -1.448161 -0.306056 3
78 -0.352255 0.242884 3
2 -0.427835 -0.629321 4
54 -0.994683 -0.607974 4
18 2.066294 0.904663 5
31 0.781440 -0.138324 5
10 -0.767944 -1.181311 6
62 -0.049937 -0.873295 6
47 -0.049937 -1.477129 7
79 -0.654574 -2.733594 7
In [41]:
train_data.info()<class 'pandas.core.frame.DataFrame'> Index: 32 entries, 4 to 79 Data columns (total 10 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 X1 32 non-null float64 1 X2 32 non-null float64 2 X3 32 non-null float64 3 X4 32 non-null float64 4 X5 32 non-null float64 5 X6 32 non-null float64 6 X7 32 non-null float64 7 X8 32 non-null float64 8 X9 32 non-null float64 9 Class 32 non-null int32 dtypes: float64(9), int32(1) memory usage: 2.6 KB
In [42]:
def scatter_matrix(samples):
# является ли подклассом?
if isinstance(samples, pd.Series): # проверка если по каким то причинам наши значения признаков имеют тип данных series то их конвертирует в data frame
samples = samples.to_frame()
d = samples - samples.mean() # вычитает из значений признаков средние значения признаков
res = np.zeros((d.shape[1], d.shape[1])) #создаёт матрицу нулей размерностью 9 на 9
# приводит к виду int: 32, 24, ...
for _, row in d.iterrows(): # проходимся циклом по каждой строке матрицы D датафрейм(там где от значений - срзнач)
col = row.to_frame() # берём строчку из таблицы D( которая сейчас имеет вид seria и мы приводим её к виду dataframe)
res += col @ col.T # матрица из нулей 9х9 берём строчку из датафрейма и умножаем её на неё же только транспонированную и так проходим по всем строчкам
return res
def classes_scatter_matrix(samples, labels): # передаём samples (значение признаков в обучающей выборке табличкой) и labels (номера классов) shape если 0 то число строк если 1 то число колонок
A = np.zeros((samples.shape[1], samples.shape[1])) # zeros создаёт матрицу shape на shape то есть создаёт матрицу число признаков на число призноков (9 на 9) заполненую нулями
for cls in labels.unique(): # переменная cls счётчик по классам принимает значения уникальных классов то есть у нас от 1 до 7
A += scatter_matrix(samples[labels == cls]) # В матрицу А прибавляем соответствующие значения из матрицы полученые в результате работы skater matrix для текущего класса
return A
cov = pd.DataFrame(classes_scatter_matrix(train_data[FEATURES], train_data.Class) / (train_data.shape[0] - train_data.Class.unique().size), \
index=FEATURES, columns=FEATURES)
print('Ковариационная матрица')
print(cov)Ковариационная матрица
X1 X2 X3 X4 X5 X6 X7 \
X1 0.131594 -0.011451 0.043813 0.088946 -0.062342 -0.066520 -0.002256
X2 -0.011451 0.080149 0.010255 -0.041460 -0.008208 0.027059 -0.004493
X3 0.043813 0.010255 0.236552 0.125222 -0.042629 -0.105374 0.015503
X4 0.088946 -0.041460 0.125222 0.302019 -0.090391 -0.126444 0.039835
X5 -0.062342 -0.008208 -0.042629 -0.090391 0.196549 0.090677 -0.022302
X6 -0.066520 0.027059 -0.105374 -0.126444 0.090677 0.266100 -0.007099
X7 -0.002256 -0.004493 0.015503 0.039835 -0.022302 -0.007099 0.084685
X8 0.007213 0.003786 -0.085371 -0.054780 0.012072 0.045550 -0.008861
X9 -0.004419 -0.001432 -0.024718 -0.030794 0.043951 0.043405 -0.034830
X8 X9
X1 0.007213 -0.004419
X2 0.003786 -0.001432
X3 -0.085371 -0.024718
X4 -0.054780 -0.030794
X5 0.012072 0.043951
X6 0.045550 0.043405
X7 -0.008861 -0.034830
X8 0.269534 0.101226
X9 0.101226 0.232664
In [43]:
lda = LinearDiscriminantAnalysis()
lda.fit(train_data[FEATURES], train_data.Class)
means = pd.DataFrame(lda.means_, index=lda.classes_, columns=FEATURES) # type: ignore
print('Средние значения')
print(means)Средние значения
X1 X2 X3 X4 X5 X6 X7 \
1 -0.251479 -0.383580 0.967795 0.355869 -0.508490 -0.747804 0.332042
2 -0.386260 -0.253140 -0.874261 -0.560531 0.042356 -0.185193 -0.155391
3 -0.311838 -0.449402 -0.167007 0.126204 -0.265843 0.326461 -0.254395
4 0.519005 0.932863 -0.167007 -0.536128 1.037944 1.370933 -0.688650
5 1.308110 -0.439770 0.456335 -1.505103 0.895683 0.278701 -0.695482
6 1.856316 -0.945477 1.209825 1.058991 3.119390 0.333262 5.340213
7 -0.510227 3.678128 -0.036859 2.092247 -1.285960 -1.175916 0.623703
X8 X9
1 -0.182201 -0.312664
2 0.252382 0.732738
3 -0.465625 -0.259701
4 -0.711259 -0.618647
5 1.423867 0.383169
6 -0.408940 -1.027303
7 -0.352255 -2.105362
X1 X2 X3 X4 X5 X6 X7 \
1 -0.251479 -0.383580 0.967795 0.355869 -0.508490 -0.747804 0.332042
2 -0.386260 -0.253140 -0.874261 -0.560531 0.042356 -0.185193 -0.155391
3 -0.311838 -0.449402 -0.167007 0.126204 -0.265843 0.326461 -0.254395
4 0.519005 0.932863 -0.167007 -0.536128 1.037944 1.370933 -0.688650
5 1.308110 -0.439770 0.456335 -1.505103 0.895683 0.278701 -0.695482
6 1.856316 -0.945477 1.209825 1.058991 3.119390 0.333262 5.340213
7 -0.510227 3.678128 -0.036859 2.092247 -1.285960 -1.175916 0.623703
X8 X9
1 -0.182201 -0.312664
2 0.252382 0.732738
3 -0.465625 -0.259701
4 -0.711259 -0.618647
5 1.423867 0.383169
6 -0.408940 -1.027303
7 -0.352255 -2.105362
In [44]:
def find_mahl_sqr_dist(centers, samples, covr): # функция принимает в себя дважды средние значения и один раз матрицу ковариций
res = pd.DataFrame(index=samples.index, columns=centers.index) # создаём новый датафрейм вверху заголовки это номера классов и слева заголовки тоже номера классов
for i in centers.index: # двойной цикл идёт по одной и той же таблице
for j in samples.index:
res[i][j] = mahalanobis(centers.loc[i], samples.loc[j], np.linalg.inv(covr)) ** 2 # вычисляется растояние махаланобиса в квадрате. np.linalg.inv(covr)-возвращает матрицу обратную матрице ковариации centers.loc[i] и samples.loc[j] возвращают i и j строки таблицы means(ср знач) и это значение записывается в ячейку ij
return res
cen_dis = find_mahl_sqr_dist(means, means, cov)
print('Расстояние Махаланобиса (обучающая выборка)')
print(cen_dis)Расстояние Махаланобиса (обучающая выборка)
1 2 3 4 5 6 \
1 0.0 22.343817 15.22444 87.857078 76.915784 580.080057
2 22.343817 0.0 15.406909 72.484578 84.871051 620.59257
3 15.22444 15.406909 0.0 60.045571 82.766139 660.598405
4 87.857078 72.484578 60.045571 0.0 82.932649 618.932157
5 76.915784 84.871051 82.766139 82.932649 0.0 562.34511
6 580.080057 620.59257 660.598405 618.932157 562.34511 0.0
7 356.469575 335.176705 355.369339 267.77654 526.670574 975.229758
7
1 356.469575
2 335.176705
3 355.369339
4 267.77654
5 526.670574
6 975.229758
7 0.0
In [45]:
classes = np.unique(train_data.Class) # Получаем массив уникальных значений классов из обучающей выборки
# Создаём список групп - для каждого класса выбираем все строки с признаками FEATURES, которые принадлежат этому классу
# Результат: список из DataFrame'ов, каждый содержит объекты одного класса
groups = [train_data[FEATURES][train_data.Class == cls] for cls in classes]
n = [len(g) for g in groups] # Создаём список n - количество объектов в каждом классе
N = sum(n) # N - общее количество объектов в обучающей выборке (сумма всех элементов списка n)
p = train_data[FEATURES].shape[1] # количество признаков
# Вычисляем объединённую (pooled) ковариационную матрицу:
# rowvar=False означает, что переменные в столбцах, ddof=1 - поправка на смещение
S_pooled = sum((ni - 1) * np.cov(g, rowvar=False, ddof=1) for g, ni in zip(groups, n)) / (N - len(classes))
# Вычисляем обратную матрицу к объединённой ковариационной матрице
inv_S = np.linalg.inv(S_pooled)
# Создаём список средних векторов для каждого класса
means = [g.mean(axis=0) for g in groups]
# Вычисляем априорные вероятности классов (prior probabilities)
priors = np.array(n) / N
# Создаём пустые словари для хранения коэффициентов дискриминантных функций
coef_stat = {} # Словарь для хранения векторов коэффициентов a для каждого класса
const_stat = {} # Словарь для хранения константных членов c для каждого класса
# Проходим по всем классам и вычисляем коэффициенты дискриминантной функции для каждого
for cls, mu, p_j in zip(classes, means, priors):
# a - вектор коэффициентов при признаках X1...X9
a = inv_S @ mu # Матричное произведение
# c - константный член дискриминантной функции
# p_j - априорная вероятность класса
c = -0.5 * mu.T @ inv_S @ mu + np.log(p_j)
# Сохраняем вычисленные коэффициенты в словари по номеру класса
coef_stat[cls] = a
const_stat[cls] = c
# Создаём DataFrame для удобного отображения коэффициентов
# Строки: X1, X2, ..., X9; Столбцы: номера классов
df_stat = pd.DataFrame(coef_stat, index=[f"X{i+1}" for i in range(p)])
# Добавляем строку с константными членами в конец таблицы
df_stat.loc["Const"] = const_stat
print("Коэффициенты дискриминантных функций:")
display(df_stat)
Коэффициенты дискриминантных функций:
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|
| X1 | -5.234825 | -3.148253 | -4.217967 | 10.327559 | 16.656458 | 35.610185 | -11.002064 |
| X2 | -7.461316 | -2.651336 | -7.584118 | 13.105596 | -11.178379 | -4.201862 | 61.821871 |
| X3 | 5.873924 | -4.182788 | 0.356458 | -1.673889 | 9.525550 | 3.394379 | -16.353702 |
| X4 | -2.309080 | -1.313906 | 0.911015 | 2.145491 | -11.656789 | -10.600878 | 22.281651 |
| X5 | -3.567483 | -1.664356 | -4.297066 | 7.670232 | 4.942687 | 32.131610 | 5.652878 |
| X6 | -1.022878 | -3.339634 | 3.245656 | 5.188590 | 2.273415 | -1.875044 | -10.339378 |
| X7 | 2.264382 | 0.112277 | -5.471875 | -7.586796 | -3.817930 | 77.090589 | 0.590885 |
| X8 | 1.619403 | -1.400612 | -1.324989 | -2.957457 | 6.123916 | -2.939071 | 1.797481 |
| X9 | -0.671626 | 4.018564 | -1.121021 | -4.541955 | -3.230166 | 2.291908 | -7.498102 |
| Const | -7.817591 | -6.087361 | -5.803707 | -23.736379 | -34.669499 | -289.317890 | -153.089698 |
In [46]:
lda = LinearDiscriminantAnalysis()
lda.fit(train_data[FEATURES], train_data.Class)
means = pd.DataFrame(lda.means_, index=lda.classes_, columns=FEATURES) #type: ignore
def LDA_predict(lda, x):# принимает в себя результаты линейного дискр анализа и значения признаков всех объектов( исходная таблица только с значениями X1...X9)
return pd.DataFrame( # функция считает распределение по классам (классификация масива тестовых векторов Х )
lda.predict(x),
columns=["Class"],
index=x.index
)
lda_predict = LDA_predict(lda, data[FEATURES])
print('Распределение по классам')
print(lda_predict)
data_to_excel['Result Lda'] = lda_predictРаспределение по классам
Class
0 3
1 2
2 4
3 2
4 1
.. ...
79 7
80 1
81 5
82 5
83 2
[84 rows x 1 columns]
In [47]:
samp_dist = find_mahl_sqr_dist(means, data[FEATURES], cov)
print('Расстояние Махланобиса')
print(samp_dist)Расстояние Махланобиса
1 2 3 4 5 6 \
0 35.548618 16.784286 11.032019 67.652466 106.380096 685.635422
1 44.397398 15.17559 28.556746 72.867432 63.671602 621.00928
2 92.049277 64.122689 58.571719 3.93205 97.092489 629.959981
3 45.929629 37.680543 49.035861 60.870665 92.950853 658.258753
4 10.384607 36.685522 29.484569 81.018346 95.923694 485.761693
.. ... ... ... ... ... ...
79 361.914041 347.840907 367.588069 300.470231 556.010044 989.726755
80 11.323179 32.848072 20.11476 85.940301 75.565423 551.195134
81 529.682935 494.161709 479.158855 331.787176 276.596841 630.305397
82 369.683371 396.241564 330.687362 230.499407 183.163732 621.077741
83 26.559382 6.11464 19.063954 43.463375 66.032665 555.489917
7
0 362.911582
1 439.360623
2 247.054034
3 248.019142
4 299.084459
.. ...
79 6.031622
80 372.395174
81 977.516974
82 891.232175
83 309.482441
[84 rows x 7 columns]
In [48]:
def LDA_predict_probab(lda, x): # принимает в себя результаты линейного дискр анализа и значения признаков всех объектов( исходная таблица только с значениями X1...X9)
return pd.DataFrame( # возвращает апостериорные вероятности классификации в соответствии с каждым классом в массиве тестовых векторов
lda.predict_proba(x),
columns=lda.classes_,
index=x.index
)
lda_post_prob = LDA_predict_probab(lda, data[FEATURES])
print('Вероятности')
print(lda_post_prob)Вероятности
1 2 3 4 5 \
0 2.724508e-06 4.313697e-02 9.568603e-01 9.702972e-14 3.778389e-22
1 3.380276e-07 9.984489e-01 1.550797e-03 7.406988e-14 7.353441e-12
2 2.201448e-19 3.402740e-13 6.825218e-12 1.000000e+00 5.894780e-21
3 1.193238e-02 9.838574e-01 4.207907e-03 2.265694e-06 2.449497e-13
4 9.998788e-01 2.592478e-06 1.186566e-04 1.530751e-16 8.876653e-20
.. ... ... ... ... ...
79 1.578468e-77 2.393959e-74 1.541647e-78 1.157359e-64 3.747166e-120
80 9.798381e-01 2.767082e-05 2.013425e-02 2.047062e-17 3.664435e-15
81 3.312640e-55 2.282701e-47 5.166393e-44 1.036504e-12 1.000000e+00
82 9.438354e-41 2.151774e-46 4.619775e-32 5.262459e-11 1.000000e+00
83 2.720782e-05 9.980490e-01 1.923774e-03 1.936053e-09 2.432529e-14
6 7
0 6.217431e-148 7.450791e-78
1 6.952668e-133 1.934890e-93
2 1.147524e-136 1.609652e-53
3 4.305525e-136 5.205061e-47
4 1.976985e-104 6.799073e-64
.. ... ...
79 2.473522e-214 1.000000e+00
80 1.915701e-118 1.283056e-79
81 1.560137e-77 6.267969e-153
82 8.094335e-96 1.757483e-154
83 1.264055e-120 3.323616e-67
[84 rows x 7 columns]
In [49]:
def wilks_lambda(samples, labels): # samples - на каждой итерации мы передаём в функцию wilks lmbd dataframe сначала с одной колонкой Х1 постепенно увеличивая число колонок пока не дойдём до конца labels - колонка с номерами классов
if isinstance(samples, pd.Series):
samples = samples.to_frame()
# определитель матрицы рассеивания
dT = np.linalg.det(scatter_matrix(samples))
# определитель классовой матрицы рассеивания
dE = np.linalg.det(classes_scatter_matrix(samples, labels))
return dE / dT
def f_p_value(lmbd, n_obj, n_sign, n_cls): # sign это число признаков вошедших в модель lmbd - значение лямбды( число) n_obj - число объектов в обучающей выборке n_cls - число классов в обучающей выборке
num = (1-lmbd)*(n_obj - n_cls - n_sign)
den = lmbd * (n_cls - 1)
f_value = num / den
p = f.sf(f_value, n_cls-1, n_obj-n_cls-n_sign)
return f_value, p
def forward(samples, labels): # samples - значение признаков в обучаюзей выборке labels - колонка с номерами классов f_in точность( это F to inter в статистике)
st_columns = ["Wilk's lmbd", "Partial lmbd", "F to enter", "P value"] # создаётся список названий колонок таблицы
n_cls = labels.unique().size #число уникальных классов нашей обучающей выборки
n_obj = samples.shape[0] # число объектов в обучающей выборке
# хранение пременных вне и в модели(е)
out = {0: pd.DataFrame(columns=st_columns, index=samples.columns, dtype=float)} # создаётся словарик ключу(ключ показывает какой у нас шаг метода) 0 ставится в соответствие дата фрейм(пустой) с колонками из переменной st_colums и индексами(строки таблицы) Х1,,,Х9
into = {0: pd.DataFrame(columns=st_columns, dtype=float)} # создаётся словарик ключу(ключ показывает какой у нас шаг метода) 0 ставится в соответствие дата фрейм(пустой) с колонками из переменной st_colums
step = 0 # шаг нашего метода
while True:
model_lmbd = wilks_lambda(samples[into[step].index], labels)
# расчёт характеристик элементов вне модели
for el in out[step].index: # el переменная счётчик по списку индексов датафрейма внутри out (Х1 Х2,,, Х9)
lmbda = wilks_lambda(samples[into[step].index.tolist() + [el]], labels) # мы с датафрейма samples берём значение из колонок в квадратных скобках список колонок который есть в inta для текущего шага + текущая колонку из счётчика el и эти значения помещаются в функцию вилкс лямбда
partial_lmbd = lmbda / model_lmbd #
f_lmbd, p_value = f_p_value(partial_lmbd, n_obj, into[step].index.size, n_cls)
out[step].loc[el] = lmbda, partial_lmbd, f_lmbd, p_value # type: ignore
# расчёт характеристик элементов в моделе
for el in into[step].index:
lmbda = wilks_lambda(samples[into[step].index.drop(el)], labels)
partial_lmbd = model_lmbd / lmbda
f_lmbd, p_value = f_p_value(partial_lmbd, n_obj, into[step].index.size-1, n_cls)
into[step].loc[el] = lmbda, partial_lmbd, f_lmbd, p_value # type: ignore
if out[step].index.size == 0:
break
# добавление нового элемента
el_to_enter = out[step]["F to enter"].idxmax()
into[step+1] = pd.concat([into[step], out[step].loc[[el_to_enter]]])
out[step+1] = out[step].drop(index=el_to_enter)
step += 1
return into, out
into, out = forward(train_data[FEATURES], train_data.Class)
print("Forward stepwise")
for i, tab in into.items():
print("Step: ", i)
print(tab, end="\n\n")
forw_stepwise = into[CNT_ENTER].index.tolist() # смотрим каждый шаг выбираем последний шаг где p-value не превышает 0,05. смотрим сколько признаков на этом шаге вошло в модель
print(forw_stepwise)Forward stepwise
Step: 0
Empty DataFrame
Columns: [Wilk's lmbd, Partial lmbd, F to enter, P value]
Index: []
Step: 1
Wilk's lmbd Partial lmbd F to enter P value
X7 1.0 0.034339 117.172492 4.570418e-17
Step: 2
Wilk's lmbd Partial lmbd F to enter P value
X7 0.055327 0.034633 111.495293 2.539633e-16
X2 0.034339 0.055802 67.682511 7.462368e-14
Step: 3
Wilk's lmbd Partial lmbd F to enter P value
X7 0.009143 0.053446 67.889761 1.792580e-13
X2 0.007714 0.063345 56.682141 1.240902e-12
X5 0.001916 0.255016 11.198399 7.409067e-06
Step: 4
Wilk's lmbd Partial lmbd F to enter P value
X7 0.002715 0.055551 62.338586 1.092269e-12
X2 0.002441 0.061807 55.658064 3.490100e-12
X5 0.000531 0.283947 9.246521 4.132537e-05
X1 0.000489 0.308693 8.211369 9.737300e-05
Step: 5
Wilk's lmbd Partial lmbd F to enter P value
X7 0.000903 0.066279 49.307354 2.671740e-11
X2 0.000960 0.062367 52.619515 1.421293e-11
X5 0.000181 0.331033 7.072970 3.186012e-04
X1 0.000168 0.356066 6.329628 6.410529e-04
X3 0.000151 0.396708 5.322606 1.782148e-03
Step: 6
Wilk's lmbd Partial lmbd F to enter P value
X7 0.000346 0.079099 38.807979 5.450728e-10
X2 0.000470 0.058251 53.890403 2.662696e-11
X5 0.000064 0.427580 4.462485 5.054419e-03
X1 0.000078 0.349448 6.205518 8.358543e-04
X3 0.000077 0.354847 6.060375 9.604737e-04
X6 0.000060 0.457851 3.947056 9.143395e-03
Step: 7
Wilk's lmbd Partial lmbd F to enter P value
X7 0.000094 0.135646 20.178341 2.661233e-07
X2 0.000196 0.064997 45.553320 2.823152e-10
X5 0.000023 0.542947 2.665705 4.761796e-02
X1 0.000039 0.322062 6.665822 6.423386e-04
X3 0.000038 0.333369 6.332330 8.660969e-04
X6 0.000028 0.453227 3.820269 1.146526e-02
X4 0.000027 0.463805 3.660921 1.382036e-02
Step: 8
Wilk's lmbd Partial lmbd F to enter P value
X7 0.000061 0.135927 19.070677 6.712073e-07
X2 0.000090 0.091516 29.780979 2.085605e-08
X5 0.000015 0.548214 2.472314 6.373804e-02
X1 0.000026 0.321306 6.336900 1.015968e-03
X3 0.000022 0.379230 4.910770 3.881752e-03
X6 0.000016 0.501614 2.980693 3.350271e-02
X4 0.000016 0.501861 2.977754 3.362441e-02
X9 0.000013 0.648374 1.626957 1.970044e-01
Step: 9
Wilk's lmbd Partial lmbd F to enter P value
X7 0.000039 0.137183 17.820289 1.783586e-06
X2 0.000053 0.100352 25.400542 1.346467e-07
X5 0.000010 0.560292 2.223553 9.119062e-02
X1 0.000015 0.369418 4.836395 4.725286e-03
X3 0.000015 0.349955 5.262941 3.138798e-03
X6 0.000010 0.529227 2.520391 6.244848e-02
X4 0.000009 0.565830 2.174061 9.722558e-02
X9 0.000008 0.643040 1.572819 2.151830e-01
X8 0.000008 0.650606 1.521583 2.304240e-01
['X7', 'X2', 'X5', 'X1', 'X3', 'X6']
In [50]:
forw_stepwise_lda = LinearDiscriminantAnalysis().fit(train_data[forw_stepwise], train_data.Class)
print("Pi: ", forw_stepwise_lda.priors_)
forw_stepwise_pred = LDA_predict(forw_stepwise_lda, data[forw_stepwise])
print("Распределение")
print(forw_stepwise_pred.head())
data_to_excel["Result forward"] = forw_stepwise_predPi: [0.1875 0.25 0.3125 0.0625 0.0625 0.0625 0.0625] Распределение Class 0 3 1 2 2 4 3 2 4 1
In [51]:
forw_stepwise_lda = LinearDiscriminantAnalysis().fit(train_data[forw_stepwise], train_data.Class)
classes = np.unique(train_data.Class)
groups = [train_data[forw_stepwise][train_data.Class == cls] for cls in classes]
n = [len(g) for g in groups]
N = sum(n)
p = train_data[forw_stepwise].shape[1]
# Общая (pooled) ковариация как в Statistica
S_pooled = sum((ni - 1) * np.cov(g, rowvar=False, ddof=1) for g, ni in zip(groups, n)) / (N - len(classes))
inv_S = np.linalg.inv(S_pooled)
# Средние по классам и априорные вероятности
means = [g.mean(axis=0) for g in groups]
priors = np.array(n) / N
# Классификационные функции (Statistica)
coef_stat = {}
const_stat = {}
for cls, mu, p_j in zip(classes, means, priors):
a = inv_S @ mu
c = -0.5 * mu.T @ inv_S @ mu + np.log(p_j)
coef_stat[cls] = a
const_stat[cls] = c
df_stat = pd.DataFrame(coef_stat, index=[f"X{i+1}" for i in range(p)])
df_stat.loc["Const"] = const_stat
print("Функции Фишера ПДАсВ:")
display(df_stat)Функции Фишера ПДАсВ:
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|
| X1 | 1.736280 | -1.743108 | -4.758759 | -5.226456 | -6.519691 | 72.846915 | 10.519265 |
| X2 | -6.145526 | -2.381514 | -8.110648 | 11.978715 | -4.772944 | 0.412523 | 51.767905 |
| X3 | -3.172120 | -0.849988 | -4.681666 | 6.484264 | 6.900413 | 34.438056 | 0.270260 |
| X4 | -5.767018 | -3.575317 | -4.300118 | 9.991338 | 13.505700 | 31.253687 | -3.123795 |
| X5 | 4.403649 | -4.202633 | 1.213440 | 0.308135 | 2.874115 | 0.482478 | -8.547150 |
| X6 | -0.755884 | -2.768609 | 2.925547 | 4.204502 | 3.521690 | -0.577544 | -13.660256 |
| Const | -7.086063 | -4.589115 | -5.259904 | -18.973678 | -19.159540 | -280.003432 | -110.069841 |
In [52]:
def backward(samples, labels):
st_columns = ["Wilk's lmbd", "Partial lmbd", "F to remove", "P value"]
n_cls = labels.unique().size
n_obj = samples.shape[0]
# хранение пременных вне и в модели(е)
into = {0: pd.DataFrame(columns=st_columns, index=samples.columns, dtype=float)}
out = {0: pd.DataFrame(columns=st_columns, dtype=float)}
step = 0
while True:
# print(step)
model_lmbd = wilks_lambda(samples[into[step].index], labels)
# расчёт характеристик элементов вне модели
for el in out[step].index:
lmbda = wilks_lambda(samples[into[step].index.tolist() + [el]], labels)
partial_lmbd = lmbda / model_lmbd
f_lmbd, p_value = f_p_value(partial_lmbd, n_obj, into[step].index.size, n_cls)
out[step].loc[el] = lmbda, partial_lmbd, f_lmbd, p_value # type: ignore
# расчёт характеристик элементов в моделе
for el in into[step].index:
lmbda = wilks_lambda(samples[into[step].index.drop(el)], labels)
partial_lmbd = model_lmbd / lmbda
f_lmbd, p_value = f_p_value(partial_lmbd, n_obj, into[step].index.size-1, n_cls)
into[step].loc[el] = lmbda, partial_lmbd, f_lmbd, p_value # type: ignore
if into[step].index.size == 0:
break
# удаление элемента
el_to_remove = into[step]["F to remove"].idxmin()
out[step+1] = pd.concat([out[step], into[step].loc[[el_to_remove]]])
into[step+1] = into[step].drop(index=el_to_remove)
step += 1
return into, out
into, out = backward(train_data[FEATURES], train_data.Class)
print("Backward stepwise")
for i, tab in into.items():
print("Step: ", i)
print(tab, end="\n\n")
back_stepwise = into[len(into) - 1 - CNT_REMOVE].index.tolist()
print(back_stepwise)Backward stepwise
Step: 0
Wilk's lmbd Partial lmbd F to remove P value
X1 0.000015 0.369418 4.836395 4.725286e-03
X2 0.000053 0.100352 25.400542 1.346467e-07
X3 0.000015 0.349955 5.262941 3.138798e-03
X4 0.000009 0.565830 2.174061 9.722558e-02
X5 0.000010 0.560292 2.223553 9.119062e-02
X6 0.000010 0.529227 2.520391 6.244848e-02
X7 0.000039 0.137183 17.820289 1.783586e-06
X8 0.000008 0.650606 1.521583 2.304240e-01
X9 0.000008 0.643040 1.572819 2.151830e-01
Step: 1
Wilk's lmbd Partial lmbd F to remove P value
X1 0.000026 0.321306 6.336900 1.015968e-03
X2 0.000090 0.091516 29.780979 2.085605e-08
X3 0.000022 0.379230 4.910770 3.881752e-03
X4 0.000016 0.501861 2.977754 3.362441e-02
X5 0.000015 0.548214 2.472314 6.373804e-02
X6 0.000016 0.501614 2.980693 3.350271e-02
X7 0.000061 0.135927 19.070677 6.712073e-07
X9 0.000013 0.648374 1.626957 1.970044e-01
Step: 2
Wilk's lmbd Partial lmbd F to remove P value
X1 0.000039 0.322062 6.665822 6.423386e-04
X2 0.000196 0.064997 45.553320 2.823152e-10
X3 0.000038 0.333369 6.332330 8.660969e-04
X4 0.000027 0.463805 3.660921 1.382036e-02
X5 0.000023 0.542947 2.665705 4.761796e-02
X6 0.000028 0.453227 3.820269 1.146526e-02
X7 0.000094 0.135646 20.178341 2.661233e-07
Step: 3
Wilk's lmbd Partial lmbd F to remove P value
X1 0.000082 0.284888 8.367167 1.277550e-04
X2 0.000340 0.068884 45.057411 1.394937e-10
X3 0.000069 0.337744 6.536070 6.131245e-04
X4 0.000064 0.365254 5.792738 1.247026e-03
X6 0.000064 0.366742 5.755716 1.293539e-03
X7 0.000347 0.067374 46.141445 1.121008e-10
Step: 4
Wilk's lmbd Partial lmbd F to remove P value
X1 0.000232 0.275052 9.224878 5.243354e-05
X2 0.000908 0.070269 46.308408 4.898174e-11
X3 0.000187 0.341989 6.734261 4.357539e-04
X4 0.000181 0.353031 6.414144 5.907362e-04
X7 0.000947 0.067399 48.429295 3.179135e-11
Step: 5
Wilk's lmbd Partial lmbd F to remove P value
X1 0.000655 0.275907 9.622842 3.073058e-05
X2 0.003223 0.056083 61.712393 1.211773e-12
X3 0.000531 0.340281 7.108733 2.619443e-04
X7 0.003487 0.051840 67.063885 5.141883e-13
Step: 6
Wilk's lmbd Partial lmbd F to remove P value
X1 0.001916 0.277240 9.993420 1.834980e-05
X2 0.009351 0.056813 63.638909 3.595709e-13
X7 0.010934 0.048586 75.064819 6.044530e-14
Step: 7
Wilk's lmbd Partial lmbd F to remove P value
X2 0.034339 0.055802 67.682511 7.462368e-14
X7 0.055327 0.034633 111.495293 2.539633e-16
Step: 8
Wilk's lmbd Partial lmbd F to remove P value
X7 1.0 0.034339 117.172492 4.570418e-17
Step: 9
Empty DataFrame
Columns: [Wilk's lmbd, Partial lmbd, F to remove, P value]
Index: []
['X1', 'X2', 'X3', 'X4', 'X6', 'X7']
In [53]:
back_stepwise_lda = LinearDiscriminantAnalysis().fit(train_data[back_stepwise], train_data.Class)
print("Pi: ", back_stepwise_lda.priors_)
back_stepwise_pred = LDA_predict(back_stepwise_lda, data[back_stepwise])
print("Распределение")
print(back_stepwise_pred.head())
data_to_excel["Result backward"] = back_stepwise_predPi: [0.1875 0.25 0.3125 0.0625 0.0625 0.0625 0.0625] Распределение Class 0 3 1 3 2 4 3 2 4 1
In [54]:
back_stepwise_lda = LinearDiscriminantAnalysis().fit(train_data[back_stepwise], train_data.Class)
classes = np.unique(train_data.Class)
groups = [train_data[back_stepwise][train_data.Class == cls] for cls in classes]
n = [len(g) for g in groups]
N = sum(n)
p = train_data[back_stepwise].shape[1]
# Общая (pooled) ковариация как в Statistica
S_pooled = sum((ni - 1) * np.cov(g, rowvar=False, ddof=1) for g, ni in zip(groups, n)) / (N - len(classes))
inv_S = np.linalg.inv(S_pooled)
# Средние по классам и априорные вероятности
means = [g.mean(axis=0) for g in groups]
priors = np.array(n) / N # можно заменить на np.ones(len(classes))/len(classes), если Statistica = равные априоры
# Классификационные функции (Statistica)
coef_stat = {}
const_stat = {}
for cls, mu, p_j in zip(classes, means, priors):
a = inv_S @ mu
c = -0.5 * mu.T @ inv_S @ mu + np.log(p_j)
coef_stat[cls] = a
const_stat[cls] = c
df_stat = pd.DataFrame(coef_stat, index=[f"X{i+1}" for i in range(p)])
df_stat.loc["Const"] = const_stat
print("Функции Фишера ПДАсИ:")
display(df_stat)Функции Фишера ПДАсИ:
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
|---|---|---|---|---|---|---|---|
| X1 | -3.771640 | -2.857112 | -3.319760 | 7.084028 | 16.643087 | 25.314092 | -12.520287 |
| X2 | -6.063287 | -2.536821 | -6.090944 | 10.678588 | -12.177667 | -15.673274 | 60.676258 |
| X3 | 4.937764 | -3.905883 | 0.357235 | 0.196979 | 7.971863 | 7.845312 | -16.315901 |
| X4 | -1.851979 | -0.914431 | 1.743032 | 1.116366 | -13.010087 | -15.719599 | 21.180522 |
| X5 | -1.973156 | -3.164523 | 1.875504 | 6.247621 | 3.328162 | 6.663230 | -10.050651 |
| X6 | 3.300529 | -1.165738 | -4.143743 | -7.414054 | -3.475812 | 69.419338 | 2.432008 |
| Const | -6.656663 | -4.606433 | -3.962780 | -16.111413 | -29.617907 | -216.567018 | -146.680111 |
In [55]:
# Сохраняем в Excel без номеров строк (индексов)
data_to_excel.to_excel(OUTPUT_PATH, index=False)