mirror of
https://github.com/ada-dmitry/msmp-mephi_python.git
synced 2026-09-24 00:50:14 +00:00
66 KiB
66 KiB
In [1]:
# 1) Названия файлов с исходными данными и пустого Excel-файла для вывода программы
# Если файлы лежат не в одной директории с программой, необходимо указать полный путь к файлам
# В целом, если конечный файл не создать, то программа его создаст сама, но лучше явно ей указать его.
INPUT_PATH = 'nstand.xlsx'
OUTPUT_PATH = 'output.xlsx'
# 2) Количество признаков в методе с включением
CNT_ENTER = 6
# 3) Количество признаков в методе с исключением
CNT_REMOVE = 6
# 4) Обучающая выборка в формате словаря, где ключ - номер класса, значение - список номеров строк, принадлежащих этому классу
TRAIN_SAMPLE = {
1: [12, 44],
2: [5, 20, 36, 37, 62, 75],
3: [32, 82],
4: [15, 19, 43],
5: [3, 9, 60, 61, 76],
6: [18, 63, 78],
7: [2, 8, 21, 22, 24, 35, 71, 81]
}
In [2]:
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)]
data_to_excel = data[FEATURES]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 [3]:
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)
data_to_excel['Train sample'] = train_data.Class X1 X2 X3 X4 X5 X6 X7 \
12 -0.303659 1.558976 -1.845236 -0.902767 0.179182 -0.175803 -0.645620
44 -0.634709 2.522227 -0.831449 0.478082 -0.433342 0.395690 -0.609662
5 -0.379991 -0.512014 0.716631 0.459036 -0.938107 -0.019480 0.063954
20 -0.443316 -0.126713 -0.310856 -0.550412 -0.622393 -0.317915 -0.290833
36 -0.385119 -0.030388 -1.434241 -0.098065 0.190525 -0.204226 -0.489801
37 -0.553875 -0.512014 0.716631 0.359043 -0.273594 0.464208 -0.333983
62 0.916945 -0.993639 0.291937 -0.240912 3.882681 1.425494 5.364185
75 -0.594438 -0.415689 -0.995847 -0.121873 -0.583637 -0.611275 -0.099056
32 3.258474 1.270000 -1.310943 1.397061 0.234952 2.307096 -0.724968
82 3.735685 -1.041802 1.127626 0.068589 1.845662 2.972992 -0.720653
15 -0.462493 2.618552 1.566020 0.730444 -0.427670 -1.998389 0.867019
19 -1.018162 1.799788 2.333210 1.592284 -0.542991 0.296719 0.145459
43 -0.306744 0.162262 1.908516 1.592284 -1.274617 -1.128461 0.658462
3 -0.560670 0.354912 -0.105358 -0.221865 1.041252 -0.715829 -0.525759
9 -0.206566 -0.223038 0.374136 0.249528 0.545940 -1.281739 0.298881
60 -0.422597 -0.367526 -0.982147 0.401897 0.647082 -0.269191 -0.046318
61 -0.306035 -0.800989 0.031640 0.782821 0.034558 0.468269 0.114296
76 0.335596 -0.752826 -0.256056 -0.302812 0.441017 -0.593003 -0.698359
18 1.542984 -0.752826 0.374136 -1.755084 1.200054 0.153592 -0.718975
63 -0.695950 -0.463851 1.168725 -0.764682 0.039284 -0.549862 -0.293230
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
8 -0.305743 -0.271201 -0.968447 0.392374 -1.313372 1.133657 -0.463432
21 -0.490466 -0.223038 -0.872548 -1.012282 0.369177 0.132275 -0.451446
22 -0.265722 -0.271201 -0.119058 -1.136083 -0.975917 -0.112360 -0.396310
24 0.014092 -0.800989 -0.091658 0.268574 -0.038226 0.004882 -0.700756
35 -0.004418 -0.126713 0.182338 -0.578982 0.809665 -0.166160 -0.458638
71 -0.351935 0.162262 -0.365655 -1.388445 1.284182 -0.347353 -0.667195
81 4.458442 -0.945477 -2.242531 -2.355039 2.465747 2.111184 -0.726406
X8 X9 Class
12 -0.654574 -0.068182 1
44 -0.919103 1.795167 1
5 0.781440 -1.397838 2
20 -0.352255 -0.080381 2
36 0.932599 1.532895 2
37 0.290172 -0.830599 2
62 -0.049937 -0.873295 2
75 0.479121 -0.388397 2
32 -0.125516 -2.044368 3
82 -0.843523 -1.483229 3
15 2.066294 -0.909891 4
19 -1.599320 -0.376199 4
43 -0.730154 -0.385348 4
3 0.894809 -0.162722 5
9 0.214592 0.355722 5
60 0.025643 1.139487 5
61 -0.654574 -0.116977 5
76 -0.919103 0.087351 5
18 2.066294 0.904663 6
63 -0.692364 0.316077 6
78 -0.352255 0.242884 6
2 -0.427835 -0.629321 7
8 -0.994683 0.431964 7
21 -0.314465 1.185232 7
22 1.952925 1.215729 7
24 0.025643 0.105649 7
35 -0.465625 -0.159672 7
71 -0.276676 -1.873587 7
81 -0.163306 -2.239547 7
In [4]:
train_data.info()<class 'pandas.core.frame.DataFrame'> Index: 29 entries, 12 to 81 Data columns (total 10 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 X1 29 non-null float64 1 X2 29 non-null float64 2 X3 29 non-null float64 3 X4 29 non-null float64 4 X5 29 non-null float64 5 X6 29 non-null float64 6 X7 29 non-null float64 7 X8 29 non-null float64 8 X9 29 non-null float64 9 Class 29 non-null int32 dtypes: float64(9), int32(1) memory usage: 3.4 KB
In [5]:
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 1.121103 -0.259221 -0.306454 -0.461469 0.767943 0.419094 0.225543
X2 -0.259221 0.454943 -0.123711 0.057824 -0.122517 -0.139440 -0.144377
X3 -0.306454 -0.123711 0.617679 0.113714 -0.034254 -0.068918 0.148946
X4 -0.461469 0.057824 0.113714 0.554228 -0.500424 0.028564 0.011445
X5 0.767943 -0.122517 -0.034254 -0.500424 1.423319 0.378642 0.826053
X6 0.419094 -0.139440 -0.068918 0.028564 0.378642 0.625262 0.278222
X7 0.225543 -0.144377 0.148946 0.011445 0.826053 0.278222 1.253289
X8 0.114257 0.186514 -0.100983 -0.245522 0.004822 -0.324924 -0.067632
X9 -0.445611 0.039838 0.016149 0.194358 -0.428706 -0.164860 -0.129618
X8 X9
X1 0.114257 -0.445611
X2 0.186514 0.039838
X3 -0.100983 0.016149
X4 -0.245522 0.194358
X5 0.004822 -0.428706
X6 -0.324924 -0.164860
X7 -0.067632 -0.129618
X8 0.944867 0.130613
X9 0.130613 0.920195
In [6]:
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.469184 2.040601 -1.338342 -0.212342 -0.127080 0.109943 -0.627641
2 -0.239966 -0.431743 -0.169291 -0.032197 0.275912 0.122801 0.702411
3 3.497080 0.114099 -0.091658 0.732825 1.040307 2.640044 -0.722810
4 -0.595800 1.526867 1.935916 1.305004 -0.748426 -0.943377 0.556980
5 -0.232055 -0.357894 -0.187557 0.181914 0.541970 -0.478299 -0.171452
6 0.138672 -0.608339 0.570500 -0.640882 0.187059 0.030598 -0.381607
7 0.412579 -0.198957 -0.682463 -0.814678 0.467129 0.524669 -0.565224
X8 X9
1 -0.786838 0.863492
2 0.346857 -0.339603
3 -0.484520 -1.763798
4 -0.087726 -0.557146
5 -0.087726 0.260573
6 0.340558 0.487875
7 -0.083003 -0.245444
In [7]:
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 7
1 0.0 34.569829 57.501138 42.238035 32.629329 28.653531 23.602533
2 34.569829 0.0 58.289833 62.028535 8.254872 6.625923 4.538061
3 57.501138 58.289833 0.0 40.907638 49.059555 48.298045 53.450568
4 42.238035 62.028535 40.907638 0.0 51.948478 50.634712 65.789422
5 32.629329 8.254872 49.059555 51.948478 0.0 9.267653 10.766549
6 28.653531 6.625923 48.298045 50.634712 9.267653 0.0 4.597995
7 23.602533 4.538061 53.450568 65.789422 10.766549 4.597995 0.0
/tmp/ipykernel_21743/2884540641.py:5: FutureWarning: ChainedAssignmentError: behaviour will change in pandas 3.0! You are setting values through chained assignment. Currently this works in certain cases, but when using Copy-on-Write (which will become the default behaviour in pandas 3.0) this will never work to update the original DataFrame or Series, because the intermediate object on which we are setting values will behave as a copy. A typical example is when you are setting values in a column of a DataFrame, like: df["col"][row_indexer] = value Use `df.loc[row_indexer, "col"] = values` instead, to perform the assignment in a single step and ensure this keeps updating the original `df`. See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy 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
In [8]:
classes = np.unique(train_data.Class)
groups = [train_data[FEATURES][train_data.Class == cls] for cls in classes]
n = [len(g) for g in groups]
N = sum(n)
p = train_data[FEATURES].shape[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]
priors = np.array(n) / N
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 | 2.362506 | -4.016844 | 8.530316 | 10.029015 | -1.338690 | -1.197691 | -3.363942 |
| X2 | 7.374780 | -3.632032 | 4.772253 | 9.883164 | -2.995975 | -1.682457 | -2.044837 |
| X3 | -0.059094 | -2.287065 | 4.369934 | 7.442596 | -1.640451 | 0.784721 | -1.990840 |
| X4 | -3.304585 | -0.642709 | 7.999418 | 4.324471 | 3.919683 | -1.867713 | -2.879039 |
| X5 | -1.829592 | 0.586150 | 1.267570 | -2.298520 | 4.057314 | 0.180565 | 0.335214 |
| X6 | -0.476711 | 2.062190 | 0.326981 | -6.260008 | -1.995700 | 1.397920 | 2.894594 |
| X7 | 1.299385 | 0.246256 | -3.019053 | 1.705241 | -2.270157 | -0.690854 | -0.711296 |
| X8 | -3.811598 | 2.103650 | -0.117392 | -3.540090 | 0.621673 | 0.793223 | 0.801712 |
| X9 | 2.248435 | -1.602480 | 0.481743 | 1.329583 | 0.090191 | 0.528417 | -0.703412 |
| Const | -12.187394 | -3.976382 | -22.378461 | -20.924723 | -4.715498 | -3.953788 | -3.741128 |
In [9]:
lda = LinearDiscriminantAnalysis()
lda.fit(train_data[FEATURES], train_data.Class)
means = pd.DataFrame(lda.means_, index=lda.classes_, columns=FEATURES)
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 7
1 7
2 7
3 5
4 2
.. ...
79 4
80 2
81 7
82 3
83 7
[84 rows x 1 columns]
In [10]:
samp_dist = find_mahl_sqr_dist(means, data[FEATURES], cov)
print('Расстояние Махланобиса')
print(samp_dist)Расстояние Махланобиса
1 2 3 4 5 6 \
0 35.668975 18.086042 97.0607 104.7747 39.287311 22.294804
1 28.121316 11.620122 63.118854 75.227086 9.049758 6.399131
2 20.971141 14.779434 53.520726 73.916814 22.538726 15.973884
3 37.74131 10.545044 55.809859 60.14781 3.88976 12.941525
4 30.772645 16.642726 37.740081 18.002732 18.136879 16.294222
.. ... ... ... ... ... ...
79 100.840777 147.154814 101.139297 37.825807 128.480876 148.89583
80 49.861969 2.600999 67.317051 76.328106 11.197973 11.170235
81 43.418763 28.329574 45.341597 90.95962 34.612055 26.013846
82 71.02434 59.211977 6.918547 53.727969 49.766347 42.226777
83 23.32652 4.672571 63.038113 71.118769 12.311471 6.113339
7
0 11.271232
1 6.021595
2 6.10029
3 12.655462
4 23.337915
.. ...
79 156.485769
80 8.699467
81 16.342598
82 52.451872
83 1.594039
[84 rows x 7 columns]
/tmp/ipykernel_21743/2884540641.py:5: FutureWarning: ChainedAssignmentError: behaviour will change in pandas 3.0! You are setting values through chained assignment. Currently this works in certain cases, but when using Copy-on-Write (which will become the default behaviour in pandas 3.0) this will never work to update the original DataFrame or Series, because the intermediate object on which we are setting values will behave as a copy. A typical example is when you are setting values in a column of a DataFrame, like: df["col"][row_indexer] = value Use `df.loc[row_indexer, "col"] = values` instead, to perform the assignment in a single step and ensure this keeps updating the original `df`. See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy 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
In [11]:
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 1.226696e-06 2.420715e-02 5.723915e-20 1.814306e-21 5.023024e-07
1 2.659498e-06 3.055692e-02 6.686223e-14 2.355065e-16 9.206105e-02
2 1.456320e-04 9.658066e-03 1.245100e-11 6.955707e-16 1.662647e-04
3 1.667614e-08 4.025419e-02 1.988658e-12 3.409394e-13 9.349886e-01
4 1.243606e-04 4.365931e-01 3.817010e-06 1.105924e-01 1.723632e-01
.. ... ... ... ... ...
79 1.381599e-14 3.635280e-24 1.190033e-14 1.000000e+00 3.438354e-20
80 1.683838e-11 9.247209e-01 2.728836e-15 4.522129e-17 1.047177e-02
81 3.283185e-07 1.862054e-03 1.255328e-07 2.338948e-17 6.707915e-05
82 1.201168e-14 1.323570e-11 1.000000e+00 1.026962e-10 1.240676e-09
83 3.967599e-06 1.337527e-01 9.446368e-15 2.492665e-16 2.445430e-03
6 7
0 1.475681e-03 9.743154e-01
1 2.078756e-01 6.695038e-01
2 2.657593e-03 9.873724e-01
3 6.072837e-03 1.868436e-02
4 2.598505e-01 2.047264e-02
.. ... ...
79 7.611163e-25 4.563366e-26
80 6.370811e-03 5.843647e-02
81 2.963578e-03 9.951068e-01
82 3.228541e-08 5.183510e-10
83 3.253971e-02 8.312582e-01
[84 rows x 7 columns]
In [12]:
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
X2 1.0 0.351167 6.774712 0.00036
Step: 2
Wilk's lmbd Partial lmbd F to enter P value
X2 0.409057 0.332926 7.012841 0.000337
X3 0.351167 0.387810 5.525047 0.001440
Step: 3
Wilk's lmbd Partial lmbd F to enter P value
X2 0.180700 0.267016 9.150317 0.000070
X3 0.155335 0.310617 7.398007 0.000285
X1 0.136186 0.354293 6.075088 0.000947
Step: 4
Wilk's lmbd Partial lmbd F to enter P value
X2 0.076802 0.290539 7.732610 0.000261
X3 0.057736 0.386482 5.026892 0.003062
X1 0.090944 0.245360 9.739531 0.000058
X4 0.048250 0.462470 3.680624 0.013503
Step: 5
Wilk's lmbd Partial lmbd F to enter P value
X2 0.043929 0.298455 7.051774 0.000554
X3 0.031571 0.415276 4.224106 0.007952
X1 0.039608 0.331012 6.063121 0.001296
X4 0.033848 0.387342 4.745103 0.004594
X6 0.022314 0.587557 2.105884 0.103208
Step: 6
Wilk's lmbd Partial lmbd F to enter P value
X2 0.032330 0.287177 7.032827 0.000684
X3 0.023916 0.388207 4.465174 0.006850
X1 0.027402 0.338816 5.529122 0.002453
X4 0.024427 0.380082 4.621208 0.005849
X6 0.017569 0.528435 2.528407 0.061822
X5 0.013111 0.708146 1.167724 0.368249
Step: 7
Wilk's lmbd Partial lmbd F to enter P value
X2 0.020494 0.280282 6.847552 0.000968
X3 0.014397 0.398984 4.016983 0.012080
X1 0.015287 0.375736 4.430515 0.007954
X4 0.015947 0.360189 4.736860 0.005908
X6 0.010992 0.522552 2.436490 0.072419
X5 0.009472 0.606445 1.730548 0.177905
X7 0.009284 0.618680 1.643587 0.199304
Step: 8
Wilk's lmbd Partial lmbd F to enter P value
X2 0.017474 0.229532 8.391751 0.000413
X3 0.010229 0.392088 3.876114 0.015467
X1 0.010836 0.370138 4.254231 0.010637
X4 0.010616 0.377814 4.117010 0.012162
X6 0.007683 0.522050 2.288811 0.090604
X5 0.006791 0.590595 1.733022 0.181304
X7 0.006568 0.610671 1.593858 0.216589
X8 0.005744 0.698248 1.080390 0.416908
Step: 9
Wilk's lmbd Partial lmbd F to enter P value
X2 0.013505 0.210590 8.746630 0.000440
X3 0.007464 0.381020 3.790575 0.018724
X1 0.007812 0.364091 4.075329 0.014227
X4 0.007528 0.377789 3.842949 0.017789
X6 0.005667 0.501852 2.316111 0.091680
X5 0.004818 0.590368 1.619003 0.214105
X7 0.004619 0.615695 1.456424 0.262463
X8 0.004523 0.628867 1.377043 0.289960
X9 0.004011 0.709119 0.957134 0.487364
['X2', 'X3', 'X1', 'X4', 'X6', 'X5']
In [13]:
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.06896552 0.20689655 0.06896552 0.10344828 0.17241379 0.10344828 0.27586207] Распределение Class 0 7 1 7 2 7 3 5 4 5
In [14]:
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 | 4.306452 | -2.095721 | 5.390990 | 7.118680 | -2.067911 | -1.219684 | -1.212971 |
| X2 | -1.263980 | -1.504646 | 4.167192 | 6.545363 | -1.605919 | 0.747939 | -1.679540 |
| X3 | -1.026588 | -2.062817 | 8.561871 | 7.236414 | -0.792322 | -0.986318 | -2.473676 |
| X4 | -1.885891 | -0.612939 | 5.537484 | 6.120972 | 1.984119 | -2.603297 | -3.506617 |
| X5 | 2.008894 | 0.529407 | 0.865689 | -4.072759 | -2.078219 | 1.006678 | 2.418692 |
| X6 | -0.392602 | 0.733885 | -1.607645 | -0.424377 | 1.842079 | -0.606497 | -0.358303 |
| Const | -8.490278 | -2.546424 | -20.097006 | -17.957060 | -3.547085 | -3.577515 | -3.450540 |
In [15]:
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.007812 0.364091 4.075329 0.014227
X2 0.013505 0.210590 8.746630 0.000440
X3 0.007464 0.381020 3.790575 0.018724
X4 0.007528 0.377789 3.842949 0.017789
X5 0.004818 0.590368 1.619003 0.214105
X6 0.005667 0.501852 2.316111 0.091680
X7 0.004619 0.615695 1.456424 0.262463
X8 0.004523 0.628867 1.377043 0.289960
X9 0.004011 0.709119 0.957134 0.487364
Step: 1
Wilk's lmbd Partial lmbd F to remove P value
X1 0.010836 0.370138 4.254231 0.010637
X2 0.017474 0.229532 8.391751 0.000413
X3 0.010229 0.392088 3.876114 0.015467
X4 0.010616 0.377814 4.117010 0.012162
X5 0.006791 0.590595 1.733022 0.181304
X6 0.007683 0.522050 2.288811 0.090604
X7 0.006568 0.610671 1.593858 0.216589
X8 0.005744 0.698248 1.080390 0.416908
Step: 2
Wilk's lmbd Partial lmbd F to remove P value
X1 0.015287 0.375736 4.430515 0.007954
X2 0.020494 0.280282 6.847552 0.000968
X3 0.014397 0.398984 4.016983 0.012080
X4 0.015947 0.360189 4.736860 0.005908
X5 0.009472 0.606445 1.730548 0.177905
X6 0.010992 0.522552 2.436490 0.072419
X7 0.009284 0.618680 1.643587 0.199304
Step: 3
Wilk's lmbd Partial lmbd F to remove P value
X1 0.027402 0.338816 5.529122 0.002453
X2 0.032330 0.287177 7.032827 0.000684
X3 0.023916 0.388207 4.465174 0.006850
X4 0.024427 0.380082 4.621208 0.005849
X5 0.013111 0.708146 1.167724 0.368249
X6 0.017569 0.528435 2.528407 0.061822
Step: 4
Wilk's lmbd Partial lmbd F to remove P value
X1 0.039608 0.331012 6.063121 0.001296
X2 0.043929 0.298455 7.051774 0.000554
X3 0.031571 0.415276 4.224106 0.007952
X4 0.033848 0.387342 4.745103 0.004594
X6 0.022314 0.587557 2.105884 0.103208
Step: 5
Wilk's lmbd Partial lmbd F to remove P value
X1 0.090944 0.245360 9.739531 0.000058
X2 0.076802 0.290539 7.732610 0.000261
X3 0.057736 0.386482 5.026892 0.003062
X4 0.048250 0.462470 3.680624 0.013503
Step: 6
Wilk's lmbd Partial lmbd F to remove P value
X1 0.136186 0.354293 6.075088 0.000947
X2 0.180700 0.267016 9.150317 0.000070
X3 0.155335 0.310617 7.398007 0.000285
Step: 7
Wilk's lmbd Partial lmbd F to remove P value
X2 0.409057 0.332926 7.012841 0.000337
X3 0.351167 0.387810 5.525047 0.001440
Step: 8
Wilk's lmbd Partial lmbd F to remove P value
X2 1.0 0.351167 6.774712 0.00036
Step: 9
Empty DataFrame
Columns: [Wilk's lmbd, Partial lmbd, F to remove, P value]
Index: []
['X1', 'X2', 'X3', 'X4', 'X5', 'X6']
In [16]:
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.06896552 0.20689655 0.06896552 0.10344828 0.17241379 0.10344828 0.27586207] Распределение Class 0 7 1 7 2 7 3 5 4 5
In [17]:
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 | -1.026588 | -2.062817 | 8.561871 | 7.236414 | -0.792322 | -0.986318 | -2.473676 |
| X2 | 4.306452 | -2.095721 | 5.390990 | 7.118680 | -2.067911 | -1.219684 | -1.212971 |
| X3 | -1.263980 | -1.504646 | 4.167192 | 6.545363 | -1.605919 | 0.747939 | -1.679540 |
| X4 | -1.885891 | -0.612939 | 5.537484 | 6.120972 | 1.984119 | -2.603297 | -3.506617 |
| X5 | -0.392602 | 0.733885 | -1.607645 | -0.424377 | 1.842079 | -0.606497 | -0.358303 |
| X6 | 2.008894 | 0.529407 | 0.865689 | -4.072759 | -2.078219 | 1.006678 | 2.418692 |
| Const | -8.490278 | -2.546424 | -20.097006 | -17.957060 | -3.547085 | -3.577515 | -3.450540 |
In [18]:
data_to_excel.to_excel(OUTPUT_PATH)