Стъпка по стъпка преглед на проекта
В тази публикация ще преминем през процеса на изграждане на регресионен модел, използвайки набора от данни за жилищното строителство в Калифорния, получен от преброяването на населението в САЩ през 1990 г. Целта е да се предскаже средната стойност на дома на даден район възможно най-точно, като се имат предвид около 20 000 точки от данни в 8 функции. Кодът по-долу импортира необходимите пакети и използвания набор от данни, както и отпечатва описание на данните, така че да добием добра представа с какво работим.
import numpy as np
import pandas as pd
from scipy.stats import t
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
from sklearn.datasets import fetch_california_housing
cal_data = fetch_california_housing(as_frame=True)
print(cal_data.DESCR)
# output ->
.. _california_housing_dataset:
California Housing dataset
--------------------------
**Data Set Characteristics:**
:Number of Instances: 20640
:Number of Attributes: 8 numeric, predictive attributes and the target
:Attribute Information:
- MedInc median income in block group
- HouseAge median house age in block group
- AveRooms average number of rooms per household
- AveBedrms average number of bedrooms per household
- Population block group population
- AveOccup average number of household members
- Latitude block group latitude
- Longitude block group longitude
:Missing Attribute Values: None
This dataset was obtained from the StatLib repository.
https://www.dcc.fc.up.pt/~ltorgo/Regression/cal_housing.html
The target variable is the median house value for California districts,
expressed in hundreds of thousands of dollars ($100,000).
This dataset was derived from the 1990 U.S. census, using one row per census
block group. A block group is the smallest geographical unit for which the U.S.
Census Bureau publishes sample data (a block group typically has a population
of 600 to 3,000 people).
An household is a group of people residing within a home. Since the average
number of rooms and bedrooms in this dataset are provided per household, these
columns may take surpinsingly large values for block groups with few households
and many empty houses, such as vacation resorts.
It can be downloaded/loaded using the
:func:`sklearn.datasets.fetch_california_housing` function.
.. topic:: References
- Pace, R. Kelley and Ronald Barry, Sparse Spatial Autoregressions,
Statistics and Probability Letters, 33 (1997) 291-297
За да предскажем точно дадена цел, трябва да идентифицираме характеристиките, които я влияят (променят) най-много. Това е основата на регресионен модел, определящ зависима (целева) променлива с помощта на независими (характерни) променливи. Ако характеристика влияе върху целта, това означава, че целевата стойност се променя (увеличава/намалява) в известна пропорция на промените (увеличава/намалява) в стойността на характеристиката. Тази връзка може да се види чрез начертаване на всички точки от данни върху 2D (X-Y) равнина, като X е стойността на характеристиката, а Y е целевата стойност. Начертаването на целта (средна стойност на къщата) върху всяка променлива на характеристиката създава следните диаграми:
Ковариация
Ако две променливи са идеално свързани, това означава, че когато едната се променя, другата ще се променя в точна пропорция всеки път, образувайки права линия, когато точките от данни са нанесени (перфектно линейни). Тъй като точките от данни се разпространяват (отдалечават се от) линията, това означава, че връзката между двете става по-малко линейна (едната не се променя правопропорционално на другата). Това разпространение се нарича ковариация и може да бъде измерено, за да се определи колко свързани са две променливи. Връзката се измерва чрез намиране на средното разстояние на точката от данни на всяка променлива от техните средни стойности, умножени заедно.
изваждаме 1 от броя на извадките, тъй като нашите данни представляват само извадка, събрана от Бюрото за преброяване на населението на САЩ от точките с данни за общото население
Повечето от графиките по-горе изглежда имат висока ковариация и нямат линейна връзка с целта, с изключение на средния доход.
Кодът по-долу изчислява ковариацията между всяка променлива на характеристиката и целта и след това показва резултата в рамка с данни на Pandas.
# calculate covariance using means of each variable
def covariance(x, y):
return np.sum(((x - x.mean())*(y - y.mean()))/(len(x)-1))
covs = [covariance(cal_data['data'][f],cal_data['target']) for f in cal_data['data'].columns]
metrics = pd.DataFrame(index=cal_data['data'].columns)
metrics['covariance'] = covs
metrics
Корелация
Данните с висока ковариация означават, че има малка връзка между променливите. Но каква стойност се счита за висока? Проблемът с ковариацията е, че тя не зависи от мащаба, което означава, че ковариацията между две променливи не може да се сравни с ковариацията между две други променливи, тъй като стойностите са в 2 различни скали. За да поставите стойностите в една и съща скала, така че да могат лесно да се сравняват, стойностите трябва да бъдат разделени на произведението на стандартните отклонения на всяка променлива. В математиката тази метрика на връзката между променливите се нарича корелация, а силата на връзката между тях е корелационният коефициент или корелацията на Пиърсън.
Кодът по-долу дефинира корелационна функция, която връща корелацията между входните променливи x и y, като използва ковариацията и стандартните отклонения на всяка от тях. След това тези данни се добавят към рамката с данни на показателите и се показват за сравнение.
def correlation(x,y):
cov = covariance(x,y)
st_dev_x = standard_deviation(x)
st_dev_y = standard_deviation(y)
corr = cov / (st_dev_x * st_dev_y)
return corr
def standard_deviation(x):
return np.sqrt(sum(np.square(x-x.mean())) / (len(x)-1))
corrs = [correlation(cal_data['data'][f],cal_data['target']) for f in cal_data['data'].columns]
metrics['correlation'] = corrs
metrics
Забележете как ковариацията между съвкупността и средната стойност на къщата е около -32 (най-високата за всички променливи на характеристиките), но корелацията е около -0,025 (почти най-ниската за всички променливи на характеристиките). Това демонстрира необходимостта от използване на корелация за сравняване на връзката на всяка променлива на характеристиката с целта, тъй като ковариацията не зависи от мащаба (всяка стойност е в различен мащаб).
Корелацията между средната стойност на къщата и средния доход е силно положителна (около 0,69), което може да се очаква. Обикновено с нарастването на доходите ви нараства и стойността на нещата, които притежавате. Съществува и известна положителна корелация между възрастта на къщата и броя на стаите в къщата със средната стойност на къщата. Отново, това има смисъл. Колкото по-нова е една къща, толкова по-висока е стойността й. Колкото повече стаи има една къща, толкова по-голяма е тя и следователно толкова по-ценна е. Но как да знаем, че тази корелация не е просто случайна? Случайността може да е отговорна за корелацията. Преди да продължим напред, трябва да определим, че корелацията не е случайна и следователно „значима“ и трябва да се използва в регресионния модел.
Статистическа значимост
Статистическата значимост е методът, чрез който се определя, че даден показател не е случаен. За да сме сигурни, че корелацията не се дължи на случайност, шансът, дължащ се на случайност, трябва да е много нисък. Шансът, дължащ се на обикновено използвания праг на случайност, е 5% (0,05), което означава, че искаме да сме 95% сигурни, че корелацията е значителна (не случайно). За да проверим дали корелацията е значима, проверяваме дали не е (това се нарича нулева хипотеза). Ако шансът за липса на корелация падне под 5% (0,05), можем да бъдем 95% или повече сигурни, че корелацията не е случайна. Това се постига чрез съпоставяне на корелацията с вероятностно разпределение (вероятностно разпределение = вероятност за всички стойности, на които една променлива може да бъде равна), за да се види дали не попада извън 95% от възможностите за случайно възникване (измерено в стандартни отклонения от липса на корелация ).
Ако корелацията се дължи на случайност, това означава, че няма корелация между двата показателя и следователно коефициентът на корелация е 0. Това, което изчисляваме, е вероятността действителната корелация да е 0 (нулева хипотеза)
Тъй като тестваме за вероятността корелацията да е 0, трябва да се вземе предвид спектърът от стойности над и под 0. Това се нарича тест с „две опашки“ (с една опашка се тества вероятността дадена стойност да е само над или под прага, а не едновременно над И под). Да бъдеш 95% сигурен, че корелацията не се дължи на случайност, е същото като да си 5% сигурен, че е така. Тъй като трябва да отчетем вероятностите както над, така и под 0, трябва да разделим 5% на два диапазона на вероятност от 2,5% и да ги поставим в края на двете страни на вероятностното разпределение. Вероятността в ниския диапазон (под 0) покрива диапазона 0%-2,5% (2,5 персентил), а високият диапазон покрива диапазона 97,5%-100% (1-97,5 персентил). Много по-лесно е за разбиране, когато се визуализира, така че по-долу е t-разпределение с 95% от вероятностите защриховани, оставяйки 2,5% незащриховани от всяка страна на разпределението. Ако вероятността се окаже извън защрихованата област, можем да сме 5% (или по-малко) сигурни, че корелацията се дължи на случайност, което означава, че можем да бъдем 95% (или повече) сигурни, че корелацията не се дължи на случайност.
Стойностите по оста x се наричат t-резултати и както споменах по-горе, представляват прогнозни стандартни отклонения от нулевата хипотеза (корелация = 0). T-резултатите при 2,5 и 97,5 персентил могат да бъдат определени с помощта на функцията за процентни точки. Тези t-резултати се наричат „критични стойности“.
Няма да навлизам в математиката зад изчисляването на тези стойности. Те са важни ценности, които трябва да знаете и разбирате, но как са извлечени е по-малко важно и би отнело цяла статия сама по себе си.
Функцията за процентна точка (PPF) изчислява t-резултат от процентил. По-долу използвам функцията PPF за t-разпределения от библиотеката scipy, за да определя критичните стойности при 2,5 и 97,5 процентила.
# calculate points on distribution where probability is less than 2.5% and 97.5% # to find middle 95% of probabilities (critical values) using the percent point # function cv_h = t(cal_data['target'].size - len(cal_data['data'].columns)-1).ppf(0.975) cv_l = t(cal_data['target'].size - len(cal_data['data'].columns)-1).ppf(0.025) # > critical values = -1.9600789769323244, 1.960078976932324
Логично е критичните стойности, определящи нашите долни и горни квантили (2,5 и 97,5), да са противоположни, тъй като вероятностното разпределение е симетрична крива и квантилите са симетрично противоположни.
За да определим къде попадат нашите корелационни стойности върху това разпределение, трябва да изчислим техните t-резултати. T-резултатът за корелация се изчислява по формулата по-долу:
Нека да проверим всеки корелационен резултат между променливите на характеристиките и целта (средна стойност на къщата), за да определим колко сигурни можем да бъдем, че не се дължат на случайност.
# map the correlations between each feature and the target to the distribution
# by calculating the correlations t-score
def t_score(r,n,p):
return r * np.sqrt((n-p) / ((1-r**2)+1e-5))
t_scores = [t_score(r,cal_data['target'].size, len(cal_data['data'].columns)+1) for r in metrics['correlation'].values]
metrics['t-scores'] = t_scores
metrics
Всички тези стойности са извън обхвата на критичната стойност, някои с голяма разлика. Можем да използваме „Функцията на кумулативното разпределение“ (CDF), за да определим процентила на t-резултата и от това да извлечем неговата p-стойност. Ако t-резултатът е по-голям от 0, трябва да извадим персентила от 1, за да изчислим правилно р-стойността (97,5 персентил оставя 2,5% от дясната опашка от 1). Тъй като провеждаме тест с две опашки (отчитайки вероятностните диапазони от двете страни на разпределението), стойността трябва да бъде умножена по 2 за точна p-стойност.
Изпълнението на кода по-долу изчислява p-стойността за всяка корелация.
def p_value(t_score,n,f):
if t_score > 0:
return 1.0 - t(n-f).cdf(t_score)
else:
return t(n-f).cdf(t_score)
p_values = [p_value(t_score, cal_data['target'].size, len(cal_data['data'].columns)+1)*2 for t_score in t_scores]
metrics['p_value'] = p_values
metrics
Колкото по-силна е корелацията между характеристиката и целта (положителна или отрицателна), толкова по-далеч от нулата е t-резултатът и толкова по-вероятно е да бъде извън обхвата на критичната стойност. Колкото по-навътре е границата на критичната стойност, толкова по-ниска е p-стойността (корелацията на вероятността е равна на 0/поради случайност). Всички корелации имат р-стойност много под 0,05 (5%), което означава, че корелациите със сигурност не са случайни.
Сега, когато знаем, че всички наши корелации са значими, можем да се доверим на тези числа и да ги използваме, за да определим кои функции ще използваме в регресионния модел.
Колинеарност
По-горе обсъдихме как средният доход, възрастта на къщата и средните стаи изглежда са силно свързани със средната стойност на къщата. Единствената друга характеристика, която изглежда има силна корелация с целта, е географската ширина на блоковата група, която има отрицателна корелация от около -0,14. Може би се питате защо географската ширина показва много по-силна (тройна) корелация от географската дължина. Въпреки че е възможно географската ширина да има много по-силно влияние върху средната стойност на къщата, отколкото дължината, това не е вероятно.
Не забравяйте, че регресията е мярка за влиянието, което независимите променливи (характеристики) имат върху зависимата променлива (целта), което означава, че ако стойността на една независима променлива (характеристика) се промени и всички други независими променливи (характеристики) останат същите, как които влияят на зависимата променлива (цел). За да бъде една променлива независима, тя не може да бъде повлияна от друга променлива. Ако една променлива е повлияна от друга, това означава, че когато едната се променя, се променя и другата, и следователно мярката за нейното влияние върху целта не може да бъде правилно изчислена. Регресионният модел все пак може да намери коефициенти, които изглежда добре предсказват целта, но те няма да представят истинското влияние върху целта и ще направят модела много чувствителен към промени (като добавяне или премахване на променливи на характеристики). Точно както използвахме корелация, за да определим влиянието на променливите на характеристиките върху целта, така може да се използва и за определяне на влиянието една върху друга. Използвайки това, което се нарича корелационна матрица, можем да определим дали две променливи на характеристиките са свързани (влияят) една на друга и следователно не са независими една от друга. Библиотеката Pandas прави това много лесно с метода corr(). Всичко, което трябва да направим, е да предадем нашата рамка от данни на всички променливи на характеристиките и целевата променлива (cal_data[‘frame’]) към функцията corr() и да покажем резултатите с помощта на топлинна карта (използвайки библиотеката Seaborn).
corr = cal_data['frame'].corr(method='pearson') sns.heatmap(data=corr.round(2), cmap='coolwarm', annot=True)
Можем да видим, че променливите за географска ширина и дължина са изключително отрицателно корелирани (-.92), което означава, че когато едната нараства, другата намалява с почти същата скорост, което ги прави зависими една от друга. Средните стаи и средните спални също имат силна положителна корелация (.85), което има смисъл. По-голямата къща обикновено има повече спални. Тази зависимост (корелация) между променливите на характеристиките се нарича колинеарност и в идеалния случай трябва да бъде намалена възможно най-много. Важно е да се разбере, че колинеарността не може да бъде елиминирана напълно от регресионния модел, тъй като променливите на характеристиките, които са свързани (корелирани) с целта, неизбежно ще споделят някаква връзка (корелация) една с друга (например възрастта на домовете в повече богатите квартали обикновено са по-ниски от тези в кварталите с по-ниски доходи).
Корелацията между две променливи на характеристики е добра индикация, че съществува колинеарност, но какво ще стане, ако променлива на характеристики е свързана с комбинацията от две или повече променливи на характеристики? Този тип корелация се нарича мултиколинеарност и не може да бъде идентифицирана чрез анализиране на нашата корелационна матрица по-горе. Има обаче метрика, която можем да използваме, която може да измери доколко дадена променлива е повлияна от другите променливи, присъстващи в регресията, наречен фактор на инфлация на вариацията (VIF). Това се постига чрез провеждане на регресия с променливата на функцията, зададена като цел (зависима променлива) и използване на останалите характеристики като независими променливи. След завършване на регресията, ние изчисляваме „коефициента на определяне“, който измерва каква част от вариацията на целевата променлива е предвидена от променливите на характеристиките (каква част от целта се обяснява с характеристиките). Оттам може да се изчисли VIF и да се определи дали съществува мултиколинеарност.
Има няколко начина за изчисляване на коефициента на детерминация. Един от начините е чрез повдигане на корелацията на квадрат, поради което коефициентът на детерминация често се нарича „r-квадрат“ (r се използва като променлива, представляваща корелация). Вече изчислихме корелацията между всички променливи по-горе. Поставянето на квадрат на стойностите и начертаването на топлинна карта на резултатите ще подчертае двойки функции, където вариантите съвпадат много (те се обясняват добре).
coeff_det = corr**2 sns.heatmap(data=coeff_det.round(2), cmap='coolwarm', annot=True)
Това всъщност не ни казва нищо ново. Корелационната матрица подчерта същите данни. Ако две характеристики са корелирани, променливостта на едната ще съответства много на другата.
Това засяга само колинеарността (между две променливи). Ако искаме да търсим възможността за мултиколинеарност между една променлива и останалите, трябва да се опитаме да предскажем въпросната променлива на характеристиката, като използваме само другите променливи на характеристиката (и член на отклонение). Точно това прави регресията. След като регресията се сближи, ние измерваме приликата между истинските данни (Y) и прогнозите чрез изчисляване на стойностите на r-квадрат и VIF. Ако дадена променлива на характеристиката може да бъде предвидена от другите, това означава, че голяма част от вариацията на данните в променливата на характеристиката вече присъства в (обяснена от) вариацията на другите и нейната стойност на r-квадрат ще бъде близка до 1. VIF разделя 1 на това, което не е обяснено от другите променливи на характеристиките (1-r-квадрат), което създава високи стойности за променливите на характеристиките, където е налице мултиколинеарност. Стойностите, по-големи от 10, се считат за значими и прагът, при който много статистици и специалисти в областта на науката за данни ще премахнат променливата от регресията. Ще преминем през няколко теста и ще сравним данните между тях, за да определим кой регресионен модел/и страдат от най-малко количество колинеарност (всички променливи на характеристиките са предимно независими) и прогнозира целта най-точно (отчита/обяснява повечето от променливостта и има нисък MSE).
Кодът по-долу очертава класа StandardScaler, който използвахме в моята предишна публикация. Добавих функции за изпълнение на регресия и изчисляване на стойностите на r-квадрат и VIF, които след това се използват във функцията collinearity_test, която ще приеме база данни с променливи на характеристиките и ще върне стойностите на r-квадрат и VIF, след като се изпълни регресия на всяка .
Функцията r-квадрат по-долу използва различен метод за изчисляване на коефициента на определяне.
class StandardScaler():
def __init__(self):
self.mean = None
self.standard_deviation = None
def fit(self, data):
self.mean = np.mean(data, axis=0)
self.standard_deviation = np.sqrt(np.mean(np.square(data-self.mean), axis=0))
def transform(self, data):
return (data-self.mean)/self.standard_deviation
def Regressor(X, Y, iterations, learning_rate):
# set Theta
theta = np.random.randn(X.shape[1],1)
mses = []
for i in range(iterations):
Yh = X.dot(theta)
residual = Y-Yh
mse = np.sum((np.square(residual)))/len(Y)
mses.append(mse)
Theta -= (learning_rate/len(Yh))*(-2*X.T.dot(residual))
return mses, theta
def r_squared(Y,predictions):
# calculate residual sum of squares
residuals = Y-predictions
SSR = np.sum(residuals**2)
# calculate total sum of squares
SST=np.sum(np.square(Y-Y.mean()))
# calculate coefficient of determination
r_sq = 1 - (SSR / SST)
return r_sq
def VIF(r_sq):
return 1/(1-r_sq)
# run a regression on an all features to calculate coefficient of determination
# and value inflation factor
def collinearity_test(data, iterations, learning_rate):
# create arrays to store R2 and VIF of each feature
r2s = []
VIFs = []
for col in data.columns:
# set X and Y data
X = data.drop(columns=col).to_numpy()
Y = data[col].to_numpy().reshape((-1,1))
# scale X
scaler = StandardScaler()
scaler.fit(X)
X = scaler.transform(X)
# add column of ones to X
X = np.column_stack((np.ones((X.shape[0],1)),X))
_, theta = Regressor(X, Y, iterations, learning_rate)
# calculate r-squared
r_sq = r_squared(Y,X.dot(theta))
r2s.append(r_sq)
# calculate VIF
vif = VIF(r_sq)
VIFs.append(vif)
return r2s, VIFs
Нека да пуснем променливите на функцията чрез функцията collinearity_test и да добавим резултатите към нашата база данни с показатели
r2s, VIFs = collinearity_test(cal_data['data'], 1000, 0.01) # create dataframe to display VIF and R2 per feature df = pd.DataFrame(index=cal_data['data'].columns) # df['Feature']=cal_data['data'].columns metrics['R2']=r2s metrics['VIF'] = VIFs metrics
Средният доход се обяснява с почти 60% от другите променливи, което е доста високо. VIF е нисък, така че засега ще го запазим там и ще видим как ще реагира на промените в модела. Средните стаи и средните спални вероятно представляват голяма част от обяснението едно в друго, особено след като те са силно свързани (0,85). VIFs също са малко по-високи. Тъй като средните стаи имат по-висок r-квадрат и VIF от средните спални, нека премахнем средните стаи от набора от данни и да видим дали това намалява стойностите на r-квадрат и VIF на средните спални. R-квадрат на географската ширина и дължина и стойностите на VIF също са високи, което потвърждава нашето подозрение, че те обясняват голяма част една от друга. Следвайки логиката в предишната стъпка, ще премахнем географската ширина от набора от данни и ще видим как реагира дължината.
Разделяне/кръстосано валидиране на тест за обучение
Първо трябва да установим базова линия, за да имаме с какво да сравним промените си. Базовата линия, която ще използваме, е модел, който използва всички променливи на характеристиките, за да предвиди целта. Преди да пуснете данните през модела, данните трябва да бъдат разделени на набори от данни за обучение и тестване. За да получите точен показател за прогнозиране, прогнозите трябва да се правят с помощта на нови данни (данни, които моделът не е виждал преди), за които истинската цел е известна. Това предпазва модела от това да се нагоди толкова добре към данните, че да не е в състояние да предвиди точно нови данни (известно като пренастройване).
Методът, по който става това, се нарича кръстосано валидиране. Накратко, кръстосаното валидиране се извършва чрез разделяне на данните в 2 набора от данни, единият се използва за обучение на модела (данни за обучение) и един, използван за тестване/валидиране на модела (данни за валидиране). Освен това можем да разделим данните за обучение на избран брой равни части, наречени „гънки“. Едно сгъване се съхранява като данни от теста, а останалите се подават към регресията като данни за обучение. Итерацията на регресия продължава, докато всяко сгъване бъде циклично като тестови данни. MSE за итерацията е средната MSE за всяко тестово сгъване. Този тип кръстосано валидиране се нарича k-кратно кръстосано валидиране и се счита за златен стандарт за тестване на модели и облекчаване на пренастройването. Процесът на k-кратно кръстосано валидиране е илюстриран по-долу за яснота.
Ето функцията cross_validate, която написах, за да капсулирам процеса по-горе.
def cross_validate(X, Y, iterations, learning_rate):
#scale data
scaler = StandardScaler()
scaler.fit(X)
X = scaler.transform(X)
# add column of ones to X
X = np.column_stack((np.ones((X.shape[0],1)),X))
# initialize theta term
theta = np.random.randn(X.shape[1],1)
# split data into train and validation set
## shuffle data
p = np.random.permutation(len(X))
X = X[p]
Y = Y[p]
## split
test_percent = 10
test_size = len(X) // test_percent
x_train = X[:-test_size]
y_train = Y[:-test_size]
x_test = X[-test_size:]
y_test = Y[-test_size:]
# split train into k folds
folds = 10
fold_size = len(x_train) // folds
# initialize array to store mses
mses = []
for i in range(iterations):
# reset fold mses array to be empty
fold_mses = []
for k in range(folds):
# get range of indices to use as test fold
test_fold = np.arange(k*fold_size,(k+1)*fold_size)
# set train and test data according to current test fold
x_train_fold = np.delete(x_train, test_fold, axis=0)
y_train_fold = np.delete(y_train, test_fold, axis=0)
x_test_fold = np.take(x_train, test_fold, axis=0)
y_test_fold = np.take(y_train, test_fold, axis=0)
# run regression on training folds
Yh = x_train_fold.dot(theta)
residual = y_train_fold - Yh
train_mse = np.sum((np.square(residual)))/len(y_train_fold)
theta -= (learning_rate/len(Yh))*(-2*x_train_fold.T.dot(residual))
# test on test fold
Yh = x_test_fold.dot(theta)
residual = y_test_fold - Yh
test_mse = np.sum((np.square(residual)))/len(y_test_fold)
fold_mses.append([train_mse, test_mse])
mses.append(np.mean(fold_mses, axis=0))
# if i % 10 == 0:
# print(f'Iteration {i}: {mses[-1]}')
# test validation data to determine true MSE
Yh = x_test.dot(theta)
residual = y_test - Yh
mse = np.sum((np.square(residual)))/len(y_test)
# print(f'Final MSE on validation data: {mse}')
return mse, mses, {'x_train':x_train, 'y_train':y_train, 'x_test':x_test, 'y_test':y_test, 'theta':theta}
Избор на функция
Сега можем да установим базова метрика чрез кръстосано валидиране на регресия, използвайки всяка променлива на функцията.
mse, mses, data = cross_validate(cal_data['data'].to_numpy(), cal_data['target'].to_numpy().reshape((-1,1)), 1000, 0.01)
# print metrics
print(f'mse={mse}')
rs = r_squared(data['y_test'], data['x_test'].dot(data['theta']))
print(f'r-squared={rs}')
mses = np.array(mses)
plt.plot(mses[:,0], label='Train')
plt.plot(mses[:,1], label='Test')
plt.scatter(x=1000, y=mse, label='Validation MSE')
plt.title("MSE's")
plt.xlabel('Iteration')
plt.ylabel('MSE')
plt.legend()
metrics['coefficients']=data['theta'][1:]
metrics
# -> mse=0.5471310653792788
# -> r-squared=0.5838323213704486
Средната квадратна грешка както за влака, така и за тестовите комплекти се сближава бързо. Линиите на конвергенция за двата набора от данни никога не се разделят, което означава, че моделът не показва признаци на пренастройване на данните за обучение. Истинският тест за точност е върху данните за валидиране, които моделът все още не е видял. MSE на набора за валидиране съвпада идеално с MSE както на влака, така и на тестовия набор в края на обучението на модела (итерация 1000), доказвайки, че моделът може да предвиди точно данни, които никога преди не е виждал. Стойността на r-квадрат показва, че моделът е в състояние да обясни 58,3% от променливостта на показателя за средна стойност на дома. Нека видим как се променят тези показатели и дали можем да подобрим някои, ако променим регресионния модел.
Преди да продължа, ще модифицирам кода, за да капсулирам целия процес до момента във функция, която ще нарека „оценява“.
def evaluate(X,Y, iterations, learning_rate):
# calculate covariance between each feature and the target
covs = [covariance(X[f],Y) for f in X.columns]
# calculate correlation between each feature and the target
corrs = [correlation(X[f],Y) for f in X.columns]
# calculate t-scores of each feature
t_scores = [t_score(r,X.shape[0], X.shape[1]) for r in corrs]
# calculate p-value of each feature
p_values = [p_value(t_score, X.shape[0], X.shape[1])*2 for t_score in t_scores]
# run collinearity test on features to get r-squared and value inflation
r2s, VIFs = collinearity_test(X, iterations)
# run cross validation of regression model
val_mse, mses, data = cross_validate(X.to_numpy(), Y.to_numpy().reshape((-1,1)), iterations, learning_rate)
# calculate r-squared of the target
rs = r_squared(data['y_test'], data['x_test'].dot(data['theta']))
# build database of metrics
metrics = pd.DataFrame(index=X.columns)
metrics['covariance'] = covs
metrics['correlation'] = corrs
metrics['t-score'] = t_scores
metrics['p-value'] = p_values
metrics['r-squared'] = r2s
metrics['VIF'] = VIFs
metrics['coefficients']=data['theta'][1:]
# return dictionary of metrics
return {'validation MSE': val_mse,
'MSEs': mses,
'r-squared': rs,
'metrics': metrics,
'training data': data}
Сега ще изпълним регресия на набор от данни, който не включва средните стаи и променливите за географската ширина, за да видим как това се отразява на модела.
data = evaluate(cal_data['data'].drop(columns=['AveRooms','Latitude']), cal_data['target'], 1000, 0.01) print(data['validation MSE']) print(data['r-squared']) data['metrics'] # -> 0.5928063257288153 # -> 0.511122753967993
MSE се увеличи малко и по същество няма колинеарност, но на цената на около 7% от обяснимостта на целта (r-квадрат). Стойността на r-квадрат на средния доход спадна до 2% от 60%, което означава, че голяма част от неговата променливост се съдържа в характеристиките на средните стаи. Това има смисъл. Кварталите с по-високи средни доходи обикновено имат по-големи къщи, а по-големите къщи имат повече стаи. Какво ще кажете да замените средните спални със средни стаи? Може би съдържа някаква променливост в целта, която е била загубена при премахването й.
data = evaluate(cal_data['data'].drop(columns=['AveBedrms', 'Latitude']), cal_data['target'], 1000, 0.01) print(data['validation MSE']) print(data['r-squared']) data['metrics'] # -> 0.6348141560950129 # -> 0.5205750454284308
По принцип получаваме същия резултат. Стойността на r-квадрат на средния доход се увеличи до 11%, което ни казва, че повече от неговата информация (променливост) присъства в характеристиките на средните стаи, което изглежда логично (Кварталът с по-висок среден доход не означава непременно, че къщите ще имат повече спални, но най-вероятно ще има повече стаи общо).
Заключение
Какво ще вземем за вкъщи? Намаляването на колинеарността на модел, така че коефициентите на променливата на характеристиките да представят точно тяхната действителна тежест (влияние) може и често намалява способността на модела да обяснява (отчита променливостта на) целта. Ако целта на модела е да предвиди целта възможно най-точно, можете да изберете да запазите всички променливи на характеристиките, за да обхванете най-голямата променливост на целта. Ако искате да знаете какви са основните движещи сили и доколко те влияят на целта, намаляването на колинеарността на модела, за да намерите тези основни характеристики, би било по-важно от постигането на най-висока точност.
Какво ще стане, ако искаме да предвидим категория вместо число или да предвидим вероятността за един или множество резултати? Тези въпроси са това, за което се използва логистичната регресия и ще бъдат тема на следващата публикация.
