В Data Science, в отличие от физики, мы часто пренебрегаем согласованностью размеров. Однако применение анализа размеров может помочь в построении более совершенных моделей.
Что мы сделаем в этой статье:
1. Изучите основы размерного анализа.
2. Научитесь применять размерный анализ к задачам регрессии.
3. Разработайте размерно-надежную и объяснимую регрессию. модель со средним показателем R² 0,93. Автор книги, в которой был найден этот набор данных, получил R² 0,89 при прогнозировании журнала используемой цели.
Как следовать этой статье:
Хотя в этой статье есть код, он не важен и не является сложным. Я бы рекомендовал сначала прочитать, игнорируя код, а затем вернуться, чтобы понять код, если хотите.
Введение
Часто все, о чем мы заботимся, - это построение модели, которая может делать хорошие прогнозы. Однако иногда мы также заботимся об объяснимости модели, и многие популярные алгоритмы машинного обучения не особо объяснимы. Когда дело доходит до физических данных (килограммы, метры, кулоны и т. Д.), Мы также можем заботиться о размерной точности. К счастью, эти измерения имеют значение; они не произвольны. Поскольку они имеют значение, отношения между переменными с такими единицами должны соответствовать правилам. Мы можем использовать эти правила для построения размерно обоснованной и объяснимой модели.
Размерный анализ
Следующий подход во многом зависит от понимания размерного анализа. Общая идея размерного анализа заключается в том, что если мы выполняем операции с величинами, имеющими единицы измерения, операции выполняются также и на самих единицах. Затем мы должны обеспечить равенство размеров, а также количественное равенство. Например, давайте посмотрим на формулу объемной плотности цилиндра.
ρ= m/π r²h
* Для следующего анализа размеров мы будем использовать стандартные единицы измерения. Радиус и высота будут считаться единицами длины и называться L. Массе будет присвоено общее обозначение M. π - безразмерная константа.
Отсюда мы можем получить единицы плотности (ρ), заданные уравнением. Чтобы избежать путаницы, я буду обозначать «имеет единицы» как:.
ρ : M/(L*L*L):M/L³
Мы также можем использовать эту идею «проверять» отношения. Чтобы данное уравнение было размерно обоснованным, обе стороны уравнения должны быть согласованы в размерном отношении.
* Примечание. При работе с данными, которые включают фактические единицы измерения, имеет смысл работать с одной и той же единицей для всех единиц такого типа (например, использовать мм для всех единиц длины). На самом деле это не улучшит предсказательную способность модели, поскольку преобразование единиц измерения будет запекаться в подгоночные коэффициенты, но сделает вещи «чище».
Гадание на отношения
Внимание: мы не используем этот процесс в примере с реальным набором данных. Тем не менее, это полезно знать, и в том же духе, что и остальная работа здесь, поэтому я решил, что добавлю это.
Эта идея также может помочь нам угадать отношения. Если это просто, мы обычно можем быстро угадать взаимосвязь и подтвердить согласованность размеров. Для более сложных отношений помогает метод размерного анализа Рэлея. Я обрисую идею, используя терминологию науки о данных «цель» и «функции».
Пусть F₁, F₂,… Fₘ будет набором прогнозных функций m. Пусть T будет целью. * Цель - найти:
T = CFⁿ¹₁F₂ⁿ² ***Fₘⁿᵐ
Учитывая, что нам известны размеры каждого объекта Fᵢ и размеры цели T, мы можем решить эту взаимосвязь размерно, подставив единицы и разбив проблему на в l линейные уравнения, где l определяется как количество уникальных базовых единиц в F₁, F₂,… Fₘ. В качестве примера я буду использовать объемную плотность цилиндра.
Примечание. Мы также можем рассмотреть условия суммирования. Если каждый отдельный член соответствует размерам в левой части уравнения. Фактически мы сделаем это позже в основном примере.
Пример
В качестве примера я буду использовать объемную плотность цилиндра. Это просто, и мы знаем решение, поэтому имейте в виду, что для решения более сложных задач потребуется больше работы и немного линейной алгебры.
Целевой показатель - плотность ρ, а размеры - м, r и h. Мы будем использовать те же базовые единицы, что и раньше, для каждого из этих значений. Итак, теперь у нас есть это:
ρ = Cmᵃrᵇhᶜ
Используя базовые единицы каждого из них, мы можем составить одно уравнение, представляющее согласованность размеров. Лучше не использовать дроби и вместо этого использовать отрицательные показатели, чтобы уменьшить количество ошибок на следующем этапе. Мы также можем пока игнорировать константу C, поскольку мы предполагаем безразмерную константу.
M¹L⁻³=MᵃLᵇLᶜ
Теперь мы разбиваем это на несколько линейных уравнений.
1 = a
-3 = b+c
Первое уравнение за нас уже решено, а второе имеет бесконечное количество решений. Так будет не всегда. Часто будет одно решение, а в других случаях решение не существует (что означает, что вам нужно взглянуть на предположения о единицах и спросить себя, существует ли такая взаимосвязь). Для более сложных отношений линейная алгебра ускорит процесс.
Но достаточно о других возможностях; давайте продолжим смотреть на это конкретное уравнение. Мы знаем, что возможно бесконечное количество решений, поэтому мы предполагаем и проверяем вероятные решения на основе знаний и опыта в предметной области. В этом примере мы понимаем, что уравнение представляет собой «объемную» часть плотности, и по опыту знаем, что V ∼ Ah, поэтому мы предполагаем, что b = -2 и c = -1.
Теперь мы снова подключаемся к этому решению, чтобы получить:
ρ = C m/(r² h)
Так в чем дело с C? Мы знаем, что оно должно быть 1 / π, но, используя этот метод, мы получаем точность только до константы. Помните, что наш вариант использования относится к области науки о данных, поэтому мы можем запустить регрессию, чтобы определить C. По сути, все, что мы сделали до сих пор, - это извлекли функцию m / (r² h), для которой мы можем запустить линейную регрессию. Далее я рассмотрю небольшой пример с реальными данными, которые не так просты, и сделаю именно это.
Мескитовый
Набор данных, который мы будем использовать в наборе данных мескитового дерева. Набор данных используется в Анализе данных с использованием регрессии и многоуровневых / иерархических моделей Гельмана и Хилла. Книжные ресурсы можно найти здесь. Конкретную главу можно найти здесь. Нашей целью будет LeafWeight. Вместо этого книга нацелена на журнал LeafWeight, который является распространенной техникой, когда продукты функций предсказывают цель. Подойдем к проблеме иначе, не требуя ведения журнала.
Сводка основных данных
Во-первых, мы начнем с вывода базовой информации о пандах.
mesquite = pd.read_csv("./mesquite.tsv", sep='\s+', header=0)
mesquite_original = mesquite.copy()
mesquite.info()
Output:
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 46 entries, 0 to 45
Data columns (total 8 columns):
Obs 46 non-null int64
Group 46 non-null object
Diam1 46 non-null float64
Diam2 46 non-null float64
TotHt 46 non-null float64
CanHt 46 non-null float64
Dens 46 non-null int64
LeafWt 46 non-null float64
dtypes: float64(5), int64(2), object(1)
memory usage: 3.0+ KB
Чтобы использовать размерность, нам необходимо понимать значение каждой функции и их размеров (это всегда нужно делать, но это особенно важно при использовании этого метода). В этом разделе мне хватит EDA, чтобы ответить на мои вопросы. Имейте в виду, что обычно требуется больше EDA.
- Diam1: диаметр кроны (покрытой листвой площади куста) в метрах, измеренный по более длинной оси куста.
- Diam2: диаметр купола, измеренный по более короткой оси.
- высота навеса: высота навеса
- общая высота: общая высота куста
- Плотность: удельная плотность растения (количество первичных стеблей на единицу растения)
- group: группа измерений (0 для первой группы, 1 для второй группы)
Давайте переименуем некоторые из этих столбцов, чтобы с ними было легче работать.
mesquite = mesquite.rename(columns =
{
'Obs':'Observation',
'TotHt':'TotalHeight',
'CanHt':'CanopyHeight',
'Dens':'Density',
'LeafWt':'LeafWeight'
})
Исследовательский анализ данных
Общая высота и высота навеса
CanopyHeight будет представлен в уравнениях как hc, а TotalHeight - как hT. Оба имеют единицы измерения L (длина).
В книгах даются следующие определения этих двух функций:
- высота навеса: высота навеса
- общая высота: общая высота куста
Однако я все еще не уверен, что понимаю разницу. Высота купола - это высота начала купола? Надеюсь, диаграмма рассеяния поможет принять решение. Я добавлю красную линию, чтобы отметить равенство между осями.
figure = plt.figure(figsize=(12,8))
axes = figure.add_subplot(1, 1, 1)
axes.scatter(mesquite.CanopyHeight, mesquite.TotalHeight, color='dimgray', alpha=0.7)
axes.set_xlabel('CanopyHeight')
axes.set_ylabel('TotalHeight')
axes.set_title('TotalHeight vs. CanopyHeight')
equal = np.linspace(0, max(max(mesquite.CanopyHeight), max(mesquite.TotalHeight)), 100)
axes.plot(equal, equal, '-', color='firebrick')
Красная линия выше представляет линию y = x. Поскольку все точки данных находятся над этой линией, мы можем сделать вывод, что TotalHeight больше или равно CanopyHeight, о чем мы уже знали. На данный момент у нас есть две гипотезы о том, что представляет собой CanopyHeight:
- Расстояние от земли до низа навеса
- Расстояние от низа козырька до верха козырька
Мы знаем, что расстояние от низа до верха купола должно быть положительным. Итак, предполагая первую интерпретацию CanopyHeight, мы знаем, что TotalHeight -CanopyHeight должно быть положительным. Давай проверим это.
total = sum((mesquite.TotalHeight-mesquite.CanopyHeight)>0)/mesquite.shape[0]
print("{0:.2f} of observations have a TotalHeight > CanopyHeight".format(total))
Output:
0.96 of observations have a TotalHeight > CanopyHeight
Поскольку некоторые наблюдения не соответствуют заявленному требованию, мы можем сделать вывод, что вторая интерпретация верна.
Плотность
Найденные примечания содержат следующее определение:
* плотность: плотность единицы растения (количество первичных стеблей на единицу растения).
Тем не менее, это не совсем понятно; Для интерпретации этого определения может потребоваться знание предметной области. Считается ли единичное наблюдение единым целым? Что такое первичный стебель? Быстрый поиск в Google не дал дополнительной информации. Давайте посмотрим на подсчеты, чтобы понять, сможем ли мы лучше понять.
mesquite['Density'].value_counts().sort_index() Output: 1 33 2 7 3 3 5 1 7 1 9 1 Name: Density, dtype: int64
«Плотность» - это строго положительное целочисленное свойство, при этом у большинства кустов «Плотность» равна 1. При просмотре изображений кустов мескита я заметил, что у многих кустов есть только один стебель от земли, а у некоторых есть несколько дополнительных стеблей. Таким образом, мы сделаем вывод, что «Плотность» - это количество стеблей, исходящих из земли. В этом случае в параметре Density должно быть указано «количество стержней» N, и оно будет представлено как p. Для практических целей мы можем рассматривать это как безразмерную константу. Вы можете думать об этом как о множителе, не имеющем общей базовой единицы. Что касается влияния на общий вес листа, однако, я не уверен, что это даст важную информацию. Возможно, нам придется проявить творческий подход к тому, как добавить эту функцию.
Построение модели
При построении этой модели мы пока будем игнорировать `Group`. Судя по описанию, нет четкой логической причины, по которой это существенно поможет, и уж точно не поможет при построении нашей размерно обоснованной модели. Это только отвлечет.
Поскольку влияние «Плотности» на цель неочевидно, мы сравним результаты оценки без ее включения с результатами оценки с ее включением.
Первая теоретическая модель
Давайте быстро определим еще несколько переменных уравнения. Целевой объект LeafWeight имеет базовые единицы массы M и будет представлен в уравнениях как WL. Два диаметра будут представлены как d₁ и d ₂. Теперь цель состоит в том, чтобы найти пространственно обоснованные отношения следующим образом:
WL = F(d₁, d₂, hT, hc, p)
Идея состоит в том, чтобы найти некоторую комбинацию характеристик, которая представляет собой объем, охватывающий густоту листьев. Коэффициент, извлеченный в регрессии, будет включать объемную плотность листьев в этом объеме.
Начнем с того, что V ∼ Ah. Мы представим эту область как площадь эллипса, Aₑₗₗᵢₚₛₑ ∼ d₁d₂. Мы также могли бы добавить безразмерный элемент, `shape`, (s = d₁ / d₂), чтобы представить асимметрию куста. Мы добавим это в более поздние модели. В этой попытке мы также включаем p, что может помочь, а может и не помочь. Итак, теперь у нас есть модель.
WL = α₀ + α pd₁d₂hc
α₀ должен иметь M единиц, а α₁ должен иметь единицы M / L³ для согласованности размеров. Теперь это наша стартовая модель для регрессии. Мы разработаем новую функцию для pd₁d₂hc и запустим регрессию с этой функцией в качестве единственной функции прогнозирования.
Первая регрессия
mesquite['model_1'] = mesquite['Density']*mesquite['Diam1']*mesquite['Diam2']*mesquite['CanopyHeight']
Теперь, когда мы создали новую функцию, запустим линейную регрессию. Вместо того, чтобы запускать линейную регрессию один раз, мы запустим загрузочную регрессию. По сути, вместо того, чтобы подгонять регрессию один раз, модель подбирается на выборках бутстрапа. Количество подгоняемых образцов указывается в параметре «samples». Затем функция выводит 95% достоверный интервал для статистики оценки: среднеквадратическая ошибка (сигма), средняя абсолютная ошибка и R². Также возвращаются результаты для каждого образца, чтобы можно было проводить сравнения с другими моделями. Использование этого метода позволяет нам видеть распределение ожидаемой производительности вместо отдельных значений. Функция ниже выполняет этот процесс.
def bootstrap_regression(data, formula, estimator, samples=100):
bootstrap_results = {}
sigmas = []
r2s = []
maes = []
n = data.shape[ 0]
bootstrap_results[ "n"] = n
for i in range(samples):
sampling_indices = pd.Index([i for i in [np.random.randint(0, n - 1) for _ in range(0, n)]])
sample = data.loc[sampling_indices]
y_sample, X_sample = dmatrices(formula, sample)
estimator.fit(X_sample, y_sample)
y_predict = estimator.predict(X_sample)
sigma = np.sqrt(mean_squared_error(y_sample, y_predict))
r2 = r2_score(y_sample, y_predict)
mae = mean_absolute_error(y_sample, y_predict)
sigmas.append(sigma)
r2s.append(r2)
maes.append(mae)
# Get a set of residuals
y, X = dmatrices(formula, data)
estimator.fit(X, y)
y_predict = estimator.predict(X)
bootstrap_results["residuals"] = y-y_predict
bootstrap_results["resampled_mae"] = pd.Series(maes)
bootstrap_results["resampled_sigma"] = pd.Series(sigmas)
bootstrap_results["resampled_r^2"] = pd.Series(r2s)
# Find credible intervals
mae_cred = np.array2string(stats.mstats.mquantiles(maes, [0.25, .975]), precision=3, separator=',')
sigma_cred = np.array2string(stats.mstats.mquantiles(sigmas, [0.25, .975]), precision=3, separator=',')
r2_cred = np.array2string(stats.mstats.mquantiles(r2s, [0.25, .975]), precision=3, separator=',')
# Print results
print("MAE Resample Mean: {0:.3f}".format(np.mean(maes)))
print("95% Credible Interval for MAE: ", mae_cred)
print("Sigma Resample Mean: {0:.3f}".format(np.mean(sigmas)))
print("95% Credible Interval for sigma: ", sigma_cred)
print("R2 Resample Mean: {0:.3f}".format(np.mean(r2s)))
print("95% Credible Interval for R2: ", r2_cred)
return bootstrap_results
Теперь мы можем использовать функцию:
formula = 'LeafWeight ~ model_1' score_1 = bootstrap_regression(mesquite, formula, LinearRegression(), samples = 500) Output: MAE Resample Mean: 229.580 95% Credible Interval for MAE: [207.448,299.46 ] Sigma Resample Mean: 288.533 95% Credible Interval for sigma: [262.388,361.327] R2 Resample Mean: 0.676 95% Credible Interval for R2: [0.529,0.924]
Учитывая, что мы использовали только одну прогностическую функцию и два подгоночных параметра, 95% доверительный интервал для R² [0,552, 0,923] является довольно хорошим. Однако напомним, что мы не были уверены, что `Density` (p) оказывает значительное влияние на нашу цель. Давайте создадим еще одну функцию, которая не включает `Density` (p), и повторно запустим регрессию начальной загрузки.
mesquite['model_2'] = mesquite['Diam1']*mesquite['Diam2']*mesquite['CanopyHeight'] formula = 'LeafWeight ~ model_2' score_2 = bootstrap_regression(mesquite, formula, LinearRegression(), samples = 500) Output: MAE Resample Mean: 143.915 95% Credible Interval for MAE: [129.532,183.581] Sigma Resample Mean: 191.788 95% Credible Interval for sigma: [170.946,248.343] R2 Resample Mean: 0.845 95% Credible Interval for R2: [0.753,0.972]
На первый взгляд, это значительное улучшение. К счастью, мы можем измерить, насколько * значимо *, поскольку у нас есть * распределение * оценок для обеих моделей. Ниже мы определим вероятность того, что модель, исключающая плотность, работает лучше, чем модель, включающая плотность, на основе выборок начальной загрузки. Мы будем называть модель, включающую «Density», как model_1, а модель без «Density» - как model_2.
difference = {}
difference['mae'] = score_2['resampled_mae'] - score_1['resampled_mae']
difference['sigma'] = score_2['resampled_sigma'] - score_1['resampled_sigma']
difference['R2'] = score_2['resampled_r^2'] - score_1['resampled_r^2']
print("P(model_2 mae < model_1 mae) = {0:.2f}".format(np.mean(difference['mae']<0)))
print("P(model_2 sigma < model_1 sigma) = {0:.2f}".format(np.mean(difference['sigma']<0)))
print("P(model_2 R2 > model_1 R2) = {0:.2f}".format(np.mean(difference['R2']>0)))
Output:
P(model_2 mae < model_1 mae) = 0.98
P(model_2 sigma < model_1 sigma) = 0.97
P(model_2 R2 > model_1 R2) = 0.77
Основываясь на этих числах, весьма вероятно, что исключение «Плотности» из нашей спроектированной функции прогнозирования дает лучшую модель. Что ж, это утверждение не совсем верно. Фактически, разница между функциями p⁰d₁d₂hc и p¹d₁d₂hc. Напомним, что в нашей модели p можно считать безразмерным. Следовательно, мы можем возвести его в любую степень и оставаться неизменными по размерам. Пока мы проверили только две действующие силы. Мы можем рассматривать эту мощность как параметр модели и использовать график проверки, чтобы найти оптимальную мощность. Поскольку это нетипичный параметр, нам придется использовать настраиваемую функцию графика проверки. Приведенную ниже функцию можно легко обобщить с помощью нескольких правок, если вы захотите это сделать.
def power_validation_plot(data, features_to_mult, feature_to_optimize, powers, fit_intercept=True,
repetitions=20, folds=5):
scores = {}
scores['p'] = []
scores['train'] = []
scores['train_std'] = []
scores['test'] = []
scores['test_std'] = []
# create template feature that is product of all input features other than the one to optimize
features_to_mult.remove(feature_to_optimize)
template_feature = None
for feature in features_to_mult:
if template_feature is None:
template_feature = data[feature]
else:
template_feature = template_feature*data[feature]
for p in powers:
# create test feature
test_feature = template_feature*data[feature_to_optimize]**p
param_score_train = []
param_score_test = []
for _ in range(repetitions):
# Perform cross validation and get scores
kf = KFold(folds, shuffle=True)
for train_index, test_index in kf.split(test_feature):
test_feature_train = pd.DataFrame(test_feature.iloc[train_index])
test_feature_test = pd.DataFrame(test_feature.iloc[test_index])
y_train = pd.DataFrame(data.loc[train_index, 'LeafWeight'])
y_test = pd.DataFrame(data.loc[test_index, 'LeafWeight'])
# Fit model
lr = LinearRegression(fit_intercept=fit_intercept)
lr.fit(test_feature_train, y_train)
# Get scores and add to param_score list
train_score = lr.score(test_feature_train, y_train)
test_score = lr.score(test_feature_test, y_test)
param_score_train.append(train_score)
param_score_test.append(test_score)
# Append param_scores to 'master' score
scores['p'].append(p)
scores['train'].append(np.mean(param_score_train))
scores['train_std'].append(np.std(param_score_train))
scores['test'].append(np.mean(param_score_test))
scores['test_std'].append(np.std(param_score_test))
scores = pd.DataFrame(scores)
# Getting bounds of confidence interval
train_upper = scores['train'] + scores['train_std']
train_lower = scores['train'] - scores['train_std']
test_upper = scores['test'] + scores['test_std']
test_lower = scores['test'] - scores['test_std']
# Create plot
figure = plt.figure(figsize=(12,8))
axes = figure.add_subplot(1, 1, 1)
axes.plot(scores["p"], scores["train"], '-', label='train', color ='firebrick')
axes.plot(scores["p"], scores["test"], '-', label='test', color = 'steelblue')
axes.fill_between(scores["p"], train_upper, train_lower, color='firebrick', alpha=0.25)
axes.fill_between(scores["p"], test_upper, test_lower, color='steelblue', alpha=0.25)
axes.legend()
axes.set_xlabel('Power')
axes.set_ylabel(r'$R^2$')
axes.set_title("Validation Curve for Power of " + feature_to_optimize)
axes.set_xticks(powers)
axes.set_ylim([0.0, 1.0])
plt.show()
plt.close()
return scores
Теперь мы можем использовать функцию:
powers = np.linspace(start=-1.0, stop=1.0, num=21) features_to_mult = ['Density', 'Diam1', 'Diam2', 'CanopyHeight'] feature_to_optimize = 'Density' validation_scores = power_validation_plot(mesquite, features_to_mult, feature_to_optimize, powers)
И тестовая кривая, и тренировочная кривая имеют пик около 0, поэтому мы можем сделать вывод, что лучше всего оставить значение `Density` (p) из извлеченного объекта, что приводит нас к следующей модели:
W_L = α₀ + α₁d₁d₂hc
Учитывая перехват
В текущем состоянии наша модель:
WL = α₀ + α₁d₁d₂hc
Помните, что цель здесь - построить размерно обоснованную и объяснимую модель. Итак, что мы знаем об этом коэффициенте α₀? Мы знаем, что он должен иметь единицы измерения M, и мы знаем, что он представляет собой LeafWeight (WL), когда измеренный диаметр или высота навеса равны 0. Другими словами, это LeafWeight (WL) без навеса. Следовательно, с физической точки зрения, α₀ должно быть равно 0 и удалено как коэффициент регрессии.
Давайте посмотрим, как наша модель будет работать при удалении α₀. Модель, которую мы сейчас рассматриваем, выглядит следующим образом:
WL = α₁d₁d₂hc
formula = 'LeafWeight ~ model_2' score_3 = bootstrap_regression(mesquite, formula, LinearRegression(fit_intercept=False), samples = 500) Output: MAE Resample Mean: 144.032 95% Credible Interval for MAE: [128.938,190.931] Sigma Resample Mean: 191.426 95% Credible Interval for sigma: [172.039,250.331] R2 Resample Mean: 0.843 95% Credible Interval for R2: [0.769,0.972]
Эти вероятные интервалы почти полностью перекрывают интервалы с момента включения перехвата. Следовательно, мы, вероятно, сможем убрать перехват из режима. Давайте повторим предыдущую кривую проверки для подтверждения с помощью тестового набора (это также подтвердит, что степень 0 уместна, если мы удалим точку пересечения).
powers = np.linspace(start=-1.0, stop=1.0, num=21) features_to_mult = ['Density', 'Diam1', 'Diam2', 'CanopyHeight'] feature_to_optimize = 'Density' validation_scores = power_validation_plot(mesquite, features_to_mult, feature_to_optimize, powers, fit_intercept=False)
Мы можем увидеть явный обрыв в R², когда мы удалим точку пересечения с модели. Это может быть просто из-за небольшого размера набора данных, но мы пока оставим перехват.
Добавление асимметрии формы к модели
В нашем теоретическом анализе до сих пор мы рассматривали только «идеальные» формы. Однако недостатки, такие как несоответствие между «Diam1» и «Diam2», также могут влиять на «LeafWeight». Мы собираемся включить это, вставив d₁ / d₂ где-нибудь в нашей предыдущей модели. Поскольку это безразмерная величина, у нас есть много вариантов. Для начала рассмотрим:
WL = α₀ + α₁(d₁/d₂)ᵐd₁d₂hc
и оптимизировать м.
mesquite['shape_1'] = mesquite['Diam1']/mesquite['Diam2'] powers = np.linspace(start=-1.0, stop=2.0, num=31) features_to_mult = ['shape_1', 'Diam1', 'Diam2', 'CanopyHeight'] feature_to_optimize = 'shape_1' validation_scores = power_validation_plot(mesquite, features_to_mult, feature_to_optimize, powers)
Похоже, что некоторые выгоды можно получить, добавив эту деталь формы с степенью 1/2. Однако на практике ими можно пренебречь. Нам нужно больше данных, чтобы знать наверняка.
Помните, что это не единственный способ учитывать форму. Существует бесконечное количество возможностей, это лишь одна из них. Попробуем по-другому:
WL = α₀ + α₁(1+d₁/d₂)ᵐd₁d₂hc
Мы будем искать оптимальную м.
mesquite['shape_2'] = 1+mesquite['Diam1']/mesquite['Diam2'] powers = np.linspace(start=-1.0, stop=2.0, num=31) features_to_mult = ['shape_2', 'Diam1', 'Diam2', 'CanopyHeight'] feature_to_optimize = 'shape_2' validation_scores = power_validation_plot(mesquite, features_to_mult, feature_to_optimize, powers)
Опять же, похоже, что здесь можно добиться некоторой выгоды; однако, опять же, прибыль может быть незначительной. Мы могли бы попробовать разные модели для добавления формы, но пока мы продолжим следующие, основываясь на выводах из графиков валидации:
WL = α₀ + α₁(d₁/d₂)^(1/2)d₁d₂hc
mesquite['CanopyHeight_Term'] = (mesquite['Diam1']/mesquite['Diam2'])**(1/2)*mesquite['Diam1']*mesquite['Diam2']*mesquite['CanopyHeight']
Давайте теперь разработаем наши достоверные интервалы для этой модели в качестве нашей новой основы.
formula = 'LeafWeight ~ CanopyHeight_Term' score_4 = bootstrap_regression(mesquite, formula, LinearRegression(), samples = 500) Output: MAE Resample Mean: 137.272 95% Credible Interval for MAE: [123.928,179.249] Sigma Resample Mean: 185.288 95% Credible Interval for sigma: [164.237,242.635] R2 Resample Mean: 0.848 95% Credible Interval for R2: [0.779,0.975]
Добавление общей высоты к модели
До сих пор мы рассматривали в модели только CanopyHeight. По-прежнему возможно, что TotalHeight также влияет на LeafWeight. Опять же, у нас есть много вариантов; нам просто нужно оставаться неизменными по размерам. Попробуем добавить его так же, как CanopyHeight:
WL = α₀ + α₁(d₁/d₂)^(1/2)d₁d₂hc + α₂d₁d₂hT
Опять же, есть много вариантов для изучения, это лишь один из них.
mesquite['TotalHeight_Term'] = mesquite['Diam1']*mesquite['Diam2']*mesquite['TotalHeight'] formula = 'LeafWeight ~ CanopyHeight_Term + TotalHeight_Term' score_5 = bootstrap_regression(mesquite, formula, LinearRegression(), samples = 500) Output: MAE Resample Mean: 121.831 95% Credible Interval for MAE: [110.548,153.799] Sigma Resample Mean: 155.733 95% Credible Interval for sigma: [140.669,196.775] R2 Resample Mean: 0.903 95% Credible Interval for R2: [0.869,0.982]
Это кажется значительным улучшением. Как было сказано ранее, мы могли бы попытаться расширить этот новый термин, чтобы увидеть, можно ли найти другие улучшения, но для краткости мы оставим его как есть.
Подтверждаем, что улучшения были внесены.
difference = {}
difference['mae'] = score_5['resampled_mae'] - score_4['resampled_mae']
difference['sigma'] = score_5['resampled_sigma'] - score_4['resampled_sigma']
difference['R2'] = score_5['resampled_r^2'] - score_4['resampled_r^2']
print("P(new mae < current mae) = {0:.2f}".format(np.mean(difference['mae']<0)))
print("P(new sigma < current sigma) = {0:.2f}".format(np.mean(difference['sigma']<0)))
print("P(new R2 > current R2) = {0:.2f}".format(np.mean(difference['R2']>0)))
Output:
P(new mae < current mae) = 0.72
P(new sigma < current sigma) = 0.78
P(new R2 > current R2) = 0.64
Таким образом, с вероятностью 0,72 мы улучшили нашу оценку средней абсолютной ошибки, с вероятностью 0,78 мы улучшили оценку среднеквадратичной ошибки, а с вероятностью 0,64 мы улучшили нашу оценку R². Так что, вероятно, мы улучшили нашу модель. Наша модель сейчас:
WL = α₀ + α₁(d₁/d₂)^(1/2)d₁d₂hc + α₂d₁d₂hT
Добавление группы
Безусловно, есть еще более объемные термины, которые стоит изучить. Однако теоретически мы могли бы бесконечно исследовать варианты. Поэтому оставим размерные особенности в покое и добавим в `Group`. Обычно мы сначала изучаем EDA, но это не руководство EDA. Мы просто закодируем переменную и добавим ее в качестве взаимодействующей функции. Модель, которую мы сейчас изучаем:
WL = α₀ + (α₁+β₁)(d₁/d₂)^(1/2)d₁d₂hc + (α₂+β₂)d₁d₂hT
Бета-версии отличны от нуля (значение определяется регрессией), если Group = MCD, и 0 в противном случае. Они представляют собой взаимодействие.
# One hot encoding mesquite['MCD'] = mesquite.Group == 'MCD' formula = 'LeafWeight ~ CanopyHeight_Term:MCD + TotalHeight_Term:MCD' score_6 = bootstrap_regression(mesquite, formula, LinearRegression(), samples = 500) Output: MAE Resample Mean: 98.372 95% Credible Interval for MAE: [ 87.01 ,129.802] Sigma Resample Mean: 127.832 95% Credible Interval for sigma: [115.746,160.091] R2 Resample Mean: 0.932 95% Credible Interval for R2: [0.906,0.988]
Это определенно кажется улучшением. Подтверждаем, что это так.
difference = {}
difference['mae'] = score_6['resampled_mae'] - score_5['resampled_mae']
difference['sigma'] = score_6['resampled_sigma'] - score_5['resampled_sigma']
difference['R2'] = score_6['resampled_r^2'] - score_5['resampled_r^2']
print("P(new mae < current mae) = {0:.2f}".format(np.mean(difference['mae']<0)))
print("P(new sigma < current sigma) = {0:.2f}".format(np.mean(difference['sigma']<0)))
print("P(new R2 > current R2) = {0:.2f}".format(np.mean(difference['R2']>0)))
Output:
P(new mae < current mae) = 0.82
P(new sigma < current sigma) = 0.84
P(new R2 > current R2) = 0.60
Опять же, похоже, что мы еще больше улучшили нашу модель. Хотя есть возможность внести дальнейшие улучшения, это будет наша последняя модель:
WL = α₀ + (α₁+β₁)(d₁/d₂)^(1/2)d₁d₂hc + (α₂+β₂)d₁d₂hT
Коэффициенты определяются путем выполнения регрессии.
Отметим, что в этом проекте мы не применяли перекрестную проверку. Обычно мы бы это сделали, но этот набор данных очень мал (n = 46), поэтому наши наборы тестов будут невероятно маленькими. Глядя на построенные нами проверочные кривые, мы понимаем, в чем проблема. Синяя область, представляющая плюс или минус одно стандартное отклонение, огромна. Кроме того, цель проекта состояла в том, чтобы показать, как построить размерно-обоснованные, объяснимые регрессионные модели; повторение перекрестных проверок на каждом этапе отвлечет от цели.
Резюме
При работе с данными, которые включают в себя единицы со смыслом, рассмотрите возможность использования анализа измерений. Путем размерного анализа и постепенного построения регрессионной модели мы можем разработать размерно-надежную, хорошо объяснимую модель. В случае мескитового примера мы смогли создать пространственно обоснованную, объяснимую модель со средним значением R², равным 0,93, улучшив оценки, полученные другими подходами на том же наборе данных.
Цитаты:
Гельман, Эндрю и Дженнифер Хилл. 2007. Анализ данных с использованием регрессии и многоуровневых / иерархических моделей. Аналитические методы социальных исследований. Кембридж: Издательство Кембриджского университета.
