Методичні нотатки про відбір інформативних ознак (feature selection)

Всім привіт!

Мене звуть Олексій. Я Data Scientist в компанії Align Technology. У цьому матеріалі я розповім вам про підходах до feature selection, які ми практикуємо в ході експериментів з аналізу даних.

У нашій компанії статистики та інженери machine learning аналізують великі обсяги клінічної інформації, пов'язані з лікуванням пацієнтів. У двох словах зміст цієї статті можна звести до вилучення цінних крупиці знання, що містяться в невеликій частці доступних нам зашумлених і надлишкових гігабайтів даних.

Дана стаття призначена для статистиків, інженерів машинного навчання та фахівців, які цікавляться питаннями виявлення залежностей у наборах даних. Також матеріал, викладений у статті, може бути цікавий широкому колу читачів, небайдужих до data mining. У матеріалі не будуть порушені питання feature engineering і, зокрема, застосування таких методів як аналіз головних компонент.




Встановлено, що набори з великою кількістю вхідних змінних можуть містити надлишкову інформацію і ведуть до переобученности моделей машинного навчання в тому випадку, якщо в моделі не вбудований регуляризатор. Стадія відбору інформативних ознак (ГІП тут і далі) часто є необхідним кроком у предобработке даних в ході експерименту.

  • У першій частині цієї статті ми зробимо огляд деяких методів відбору ознак і розглянемо теоретичні моменти у цьому процесі. Даний розділ є скоріше систематизацією наших поглядів.
  • У другій частині статті на прикладі штучного набору даних ми поекспериментуємо з відбором інформативних ознак і зробимо порівняння результатів.
  • У третій частині ми висвітлимо теорію і практику використання заходів з теорії інформації в застосуванні до обговорюваної задачі. Метод, представлений в цьому розділі, володіє новизною, однак він, як і раніше, вимагає проведення численних перевірок на різних наборах даних.


Експерименти, проведені в статті, не претендують на наукову закінченість в силу того, що аналітичні обґрунтування оптимальності того чи іншого методу наведено не будуть, і читач надсилається до першоджерел за більш детальним і математично точним викладом. Крім цього, disclaimer полягає в тому, що на інших даних цінність того чи іншого методу зміниться, що і робить задачу в цілому інтелектуально привабливою.

В кінці статті будуть зведені результати експериментів і зроблено посилання на повний код R на Git.

Я висловлюю подяку всім людям, хто прочитав матеріал до публікації і зробили його краще, зокрема, Владу Щербінін і Олексію Селезньову.

1) Способи і методи відбору інформативних ознак.

Давайте розглянемо загальний підхід до класифікації методів ГІП, звернувшись до Вікі:

Алгоритми відбору інформативних ознак можуть бути представлені наступними групами: Wrappers (обгорткові), Filters (фільтрувальної) і Embedded (вбудовані в машини). (Я залишу без точного перекладу ці терміни на увазі розмитості їх звучання для російськомовного співтовариства — прим. моє.)

Обгорткові алгоритми створюють поднаборы, використовуючи пошук в просторі можливих вхідних змінних і оцінюють отримані поднаборы входів шляхом навчання повної моделі на наявних даних. Обгорткові алгоритми можуть бути дуже дорогими і ризикують перенавчати модель. (Якщо не використовується валидационная вибірка — прим. моє.)

Фільтрувальної алгоритми схожі на обгорткові в тому, що вони також виконують пошук під наборів вхідних даних, але, замість запуску повної моделі, важливість поднабора для вихідної змінної оцінюється з допомогою більш простого (фильтровочного) алгоритму.

Вбудовані в машини алгоритми оцінюють важливість вхідних ознак за допомогою евристики, заздалегідь закладеної в навчання.


Джерело.

Приклади.

Оберточным алгоритмом ГІП можна назвати комбінацію методів, що включає пошук поднабора вхідних змінних з подальшим навчанням, наприклад, random forest, на відібраних даних та оцінкою його помилки на кроссвалидации. Тобто для кожної ітерації ми будемо навчати цілу машину (вже фактично готовий до подальшого використання).

Фильтровочным алгоритмом ГІП можна назвати перебір вхідних змінних, доповнюється проведенням статистичного тесту на залежність між вибраними змінними та виходом. Якщо наші входи і вихід категориальны, то можливо провести тест хі-квадрат на незалежність між входом (або комбінованим набором входів) і виходом з оцінкою p-value і, як наслідок, бінарним висновком про значущість або незначущості відібраного набору ознак. Іншими прикладами фільтрувальної алгоритмів можна вважати:
  • лінійну кореляцію входу і виходу;
  • статистичний тест на різницю в середніх у разі категоріальних входів і безперервного виходу;
  • F-критерій (дисперсійний аналіз).
Вбудованим алгоритмом ГІП є, наприклад, значення p-values, відповідні коефіцієнтів лінійної регресії. У даному випадку p-value дозволяє також зробити бінарний висновок про значимому отличиии коефіцієнта від нуля. Якщо отмасштабировать всі входи моделі, то модулі ваг можна трактувати як показники важливості. Також можна використовувати R^2 моделі — міру пояснення дисперсії процесу змодельованими значеннями. Ще одним прикладом може служити функція оцінки важливості вхідних змінних, вбудована в random forest. Крім того, можна використовувати модулі ваг, що відповідають входам в штучної нейронної мережі. Цим список не вичерпується.

На цьому етапі важливо зрозуміти, що це розмежування, по суті, вказує на відмінність фітнес-функцій ГІП, тобто міри релевантності знайденого поднабора вхідних ознак по відношенню до розв'язуваної задачі. Згодом ми ще повернемося до питання вибору фітнес-функції.

Тепер, коли ми трохи орієнтуємося в основних групах методів ГІП, пропоную звернути увагу на те, якими методами здійснюється саме перебір під наборів вхідних змінних. Знову звернемося до сторінки Вікі:

Підходи до перебору включають:
  • Повний перебір
  • Перший кращий кандидат

  • Імітація відпалу
  • Генетичний алгоритм
  • Жадібний пошук на включення
  • Жадібний пошук на виключення
  • Particle swarm optimization
  • Targeted projection pursuit
  • Scatter Search
  • Variable Neighborhood Search


Джерело.

Я навмисне не став перекладати назви деяких алгоритмів, так як не знайомий з їх російськомовної інтерпретацією.

Пошук поднабора предикторів — це дискретна задача, так як на виході ми отримуємо вектор цілих чисел, що позначають індекси входів, виду:

входи: 1 2 3 4 5 6… 1000
відбір: 0 0 1 1 1 0… 1

Ми повернемося до цієї особливості пізніше і проілюструємо, до чого вона веде практично.

Від того, як буде налаштований пошук поднабора вхідних ознак, сильно залежить результат всього експерименту. Для того щоб інтуїтивно зрозуміти основна відмінність в цих підходах, я пропоную читачеві поділити їх на дві групи: жадібні (greedy) та нежадібні (non-greedy).

Жадібні алгоритми пошуку.


Вони часто використовуються, так як швидкі і дають хороший результат у багатьох завданнях. Жадібність алгоритму полягає в тому, що якщо був обраний (або виключений) один з кандидатів на входження в фінальний підмножину, в подальшому він у ньому залишається (в разі жадібного включення) чи назавжди відсутня (у разі жадібного винятки). Таким чином, якщо на попередніх ітераціях був обраний кандидат А, на пізніх ітераціях підмножиною завжди буде включати його та інших кандидатів, які разом з А показують поліпшення метрики важливості поднабора для вихідної змінної. Зворотна ситуація для виключення: якщо був прибраний кандидат А, тому що після його виключення метрика важливості найменш постраждала чи покращилася, дослідник не отримає інформації про метриці важливості, де в поднаборе А є й інші, виключені пізніше кандидати.

Якщо проводити паралель з пошуком максимуму (мінімуму) в багатовимірному просторі, жадібний алгоритм буде застрягати в локальному мінімумі, якщо такий є, або швидко знаходити оптимальне рішення, якщо є єдиний мінімум, і він глобальний.

З іншого боку, перебір всіх жадібних варіантів здійснюється відносно швидко і дозволяє врахувати деякі взаємодії між входами.

Приклади жадібних алгоритмів: жадібне включення (forward selection; step forward) і жадібне виняток (backward elimination; step backward). На цьому список не обмежується.

Нежадібні алгоритми пошуку.


Принцип роботи нежадібних алгоритмів передбачає можливість відкинути цілком або частково сформовані поднаборы ознак, комбінації під наборів між собою і внести випадкові зміни в поднаборы з метою уникнення локального мінімуму.

Якщо провести аналогію з пошуком максимуму (мінімуму) значення фітнес-функції в багатовимірному просторі, нежадібні алгоритми розглядають набагато більше сусідніх точок і навіть можуть робити великі стрибки у випадкові області.

Можна представити графічно роботу цих видів алгоритмів.

Спочатку жадібне включення:


Тепер нежадібних — стохастичний — пошук:


В обох випадках потрібно відібрати одну найкращу комбінацію з двох входів, чиї індекси відкладені по осях графіка. Жадібний алгоритм починає з того, що відбирає один кращий вхід, перебираючи індекси по горизонталі. А потім додає до відібраному входу другий вхід так, щоб їх сумарна релевантність була максимальна. Але виходить, що з всіх можливих комбінацій він повністю проходить тільки 1/37 частина. При додаванні ще одного виміру, кількість пройдених осередків стане ще менше: приблизно, 1/37^2.

При цьому можлива практична ситуація, коли кращою є та комбінація, яку знайшов жадібний алгоритм. Це може статися, якщо кожен з двох входів окремо не показує кращу релевантність задачі (і вони не будуть відібрані на першому кроці), але їх взаємодія максимізує метрику релевантності.

Нежадібних алгоритм шукає набагато довше:

(O) = 2^n


і перевіряє більше можливих комбінацій входів. Але він має шанс знайти, можливо, ще більш хороший підмножиною входів завдяки размашистому пошуку відразу у всіх змінах завдання.

Є алгоритми пошуку, що виходять за рамки встановленої дихотомії жадібний / нежадібних.

Прикладом такого відокремленого пошуку буде почерговий перебір входів окремо і оцінка їх окремої важливості для вихідної змінної. З цього, до речі, починається перша хвиля в алгоритмі жадібного включення змінних. Але що ж виходить хорошого з таким перебором, за винятком того, що він дуже швидкий? Кожна вхідна змінна починає існувати «у вакуумі», тобто без урахування впливу інших входів на зв'язок між обраним входом і виходом. При цьому, після завершення роботи алгоритму отриманий список входів із зазначенням їх індивідуальної важливості для виходу всього лише дає інформацію про індивідуальної значущості кожного предиктора. При об'єднанні певної кількості найважливіших предикторів, згідно з таким списком, можна отримати ряд проблем:

  • надмірність (у разі кореляції предикторів між собою);
  • недостатність інформації з-за неврахування взаємодій предикторів на етапі відбору;
  • розмитість межі, вище якої потрібно взяти предиктори.


Як можна бачити, проблеми не самі тривіальні.

Основне питання в задачі ГІП формулюється як оптимальна комбінація методу пошуку поднабора і фітнес-функції.

Розглянемо це твердження детальніше. Наше завдання ГІП може бути описана двома гіпотезами:

а) поверхня помилки проста або складна;
б) у даних є прості залежності або складні.

Залежно від відповідей на ці питання слід зупинити вибір на конкретній зв'язці методу пошуку і методу визначення релевантності відібраних ознак.

Поверхня помилки.

Приклад нескладної поверхні:
image

Источник.

Тут ми вибираємо комбінацію з двох входів, визначаючи їх релевантність виходу і спускаємося по плавній поверхні в напрямку градієнта, майже напевно потрапляючи в точку оптимуму.

Приклад складної поверхні:
image

Джерело.

У цьому випадку, вирішуючи ту ж задачу, ми зустрічаємося з безліччю локальних мінімумів, вибратися з яких жадібний алгоритм не зможе. У той же час, алгоритм зі стохастичним пошуком маємо підвищений шанс знайти більш точне рішення.

Раніше ми згадали, що пошук поднабора предикторів — це дискретна задача. Якщо залежність виходу від входів включає в себе взаємодії, при переході з однієї точки простору в сусідню можна спостерігати різкі скачки значення фітнес-функції. Поверхня помилки в нашому випадку часто не гладка, не диференціюється:
image

Це приклад пошуку поднабора з двох входів і відповідне значення функції релевантності поднабора вихідної змінної. Видно, що поверхня далеко не гладка, з піками, а також включає в себе порізані плато приблизно однакових значень. Нічний кошмар для методів жадібного градієнтного спуску.

Залежності.

З ростом числа вимірювань завдання підвищується теоретичекий шанс того, що залежність вихідної змінної має вельми складну структуру і задіє безліч входів. Крім того, завимость може мати лінійний вигляд, так і бути нелінійною. Якщо залежність увазі взаємодію предикторів та нелінійний вигляд, знайти її можна, тільки враховуючи обидва ці моменту, застосувавши, наприклад, навчання random forest або нейронну мережу. Якщо залежність проста, лінійна, включає в себе лише невелику частину всіх предикторів, підхід до її знаходження — і, як наслідок, до ГІП, — можна звести до почергового включення в модель лінійної регресії по 1 або кілька входів з оцінкою якості моделі.

Приклад простої залежності:
image

В даному випадку залежність значення по осі output від значень input1 і input2 описується площиною у просторі:
output = input1*10 + input2*10
Модель такої залежності дуже проста і може бути апроксимована лінійною регресією.

Приклад складної залежності:
image

Ця нелінійна залежність вже не може бути виявлена побудовою лінійної моделі. Її вигляд такий:
output = input1^2 + input2^2


Також необхідно врахувати розмірність задачі.

Якщо кількість вхідних змінних велике, пошук оптимального поднабора стохастичними методами (нежадными) може бути дуже дорогим в силу того, що загальна кількість всіх можливих під наборів задається відношенням
m = 2 ^ n,
де n — кількість всіх вхідних ознак.


Відповідно, пошук мінімуму в такому різноманітті може бути дуже довгий. З іншого боку, застосування жадібного пошуку дозволить зробити перше наближення за розумний час, навіть якщо це буде локальним мінімумом і дослідник буде знати про це.

В умовах відсутності предметного знання про досліджуваному явищі неможливо заздалегідь сказати, наскільки складною буде залежність вхідних змінних і виходу і яку кількість входів буде оптимально для знаходження наближеного або точного рішення задачі відбору оптимального поднабора входів. Також важко передбачити поверхню помилки при ГІП гладкою і простою або складною і порізаною.

Також ми завжди обмежені в ресурсах і повинні приймати найбільш оптимальні рішення. В якості невеликої допомоги при виробленні підходу до ГІП можна використовувати наступну таблицю:
image

Таким чином, у нас завжди є можливість розглянути декілька комбінацій методів пошуку поднабора входів і фітнес-функції релевантності. Найдорожча та, напевно, найбільш ефективна комбінація — це жадібний пошук і обгортковий (wrapper) фітнес-фукнцію. Дозволяє уникати локального мінімуму, при цьому надаючи найбільш точну міру релевантності відібраних входів, так як на кожній ітерації є навчена модель (і її точність на валідації).

Найменш дорогий, але не завжди наимее ефективний підхід, це жадібний пошук в парі з фільтровочною (filter) функцією, яка може бути статистичним тестом, коефіцієнтом кореляції або кількістю взаємної інформації.

Крім того, вбудовані (embedded) методи дозволяють відразу після навчання моделі відсіяти ряд непотрібних з точки зору алгоритму входів, не втративши істотно в точності моделювання.

Хорошим тоном буде спроба кілька разів розв'язати задачу різними способами і вибрати найкращий.

Загалом кажучи, відбір інформативних ознак — це підбір найкращої комбінації методу пошуку в багатовимірному просторі і найбільш оптимальної фітнес-функції релевантності відібраного поднабора по відношенню до вихідної змінної.

image

Джерело.

2) Експерименти з відбору інформативних ознак на синтетичних даних.



Наш піддослідний набір даних («Стенфордський кролик»):
image

Просто я люблю кроликів.

Ми будемо дивитися на залежність висоти точки (вісь Z) від широти і довготи. При цьому в набір я додаю 10 шумових змінних з розподілом приблизно відповідним суміші двох інформативних входів (X і Y), але ніяк не пов'язаних зі змінною Z.

Подивимося на гістограми щільності розподілу для змінних X, Y, Z і однією з шумових змінних.
image

image

image

image

Видно, що розподілу з довільними параметрами. При цьому всі шумові змінні розподілені так, що у них є невеликий пік в деякій області значень.

Далі набір даних буде випадково поділений на дві частини: навчання і валідація.

Підготовка даних.

Код
library(onion)
data(bunny)

#XYZ point plot
open3d()
points3d(bunny, col=8, size=0.1)

bunny_dat <- as.data.frame(bunny)

inputs <- append(as.numeric(bunny_dat$x),
as.numeric(bunny_dat$y))

for (i in 1:10){

naming <- paste('input_noise_', i, sep = ")
bunny_dat[, eval(naming)] <- inputs[sample(length(inputs), nrow(bunny_dat), replace = T)]


}



Експеримент №1: жадібний пошук поднабора входів з лінійною функцією оцінки важливості (як фітнес-функції буде використовуватися варіант wrapper — оцінка коефіцієнта детермінації навченої моделі на валідаційної вибірці).

Код
### greedy search with filter function

library(FSelector)

sampleA <- bunny_dat[sample(nrow(bunny_dat), nrow(bunny_dat) / 2, replace = F), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

sampleB <- bunny_dat[!row.names(bunny_dat) %in% rownames(sampleA), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

linear_fit <- function(subset){
dat <- sampleA[, c(subset, "z")]

lm_m <- lm(formula = z ~.,
data = dat,
model = T)

lm_predict <- predict(lm_m,
newdata = sampleB)

r_sq_validate <- 1 - sum((sampleB$z - lm_predict) ^ 2) / sum((sampleB$z - mean(sampleB$z)) ^ 2)

print(subset)
print(r_sq_validate)
return(r_sq_validate)
}

#subset <- backward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = linear_fit)
subset <- forward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = linear_fit)



Результат:

> subset
<b>[1] "x" "y" "input_noise_2" "input_noise_5" "input_noise_6" "input_noise_8" "input_noise_9"</b>


Виявилися включені шумові змінні.

Подивимося на навчену модель:

> summary(lm_m)

Call:
lm(formula = z ~ ., data = dat, model = T)

Residuals:
Min 1Q Median 3Q Max 
-0.060613 -0.022650 -0.000173 0.024939 0.048544 

Coefficients:
Estimate Std. Error t value Pr(>|t|) 
(Intercept) 0.0232453 0.0005581 41.651 < 2e-16 ***
<b>x -0.0257686 0.0052998 -4.862 1.17 e-06 ***
y -0.1572786 0.0052585 -29.910 < 2e-16 ***</b>
input_noise_2 -0.0017249 0.0027680 -0.623 0.533 
input_noise_5 -0.0027391 0.0027848 -0.984 0.325 
input_noise_6 0.0032417 0.0027907 1.162 0.245 
input_noise_8 0.0044998 0.0027723 1.623 0.105 
input_noise_9 0.0006839 0.0027808 0.246 0.806 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1

Residual standard error: 0.02742 on 17965 degrees of freedom
Multiple R-squared: 0.04937, Adjusted R-squared: 0.049 
F-statistic: 133.3 on 7 and 17965 DF, p-value: < 2.2 e-16


Бачимо, що насправді статистичну значимість приймають тільки наші оригінальні входи і вільний член рівняння.

А тепер проведемо жадібне виключення змінних.

Код
subset <- backward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = linear_fit)
#subset <- forward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = linear_fit)



Результат:
> subset
<b>[1] "x" "y" "input_noise_2" "input_noise_5" "input_noise_6" "input_noise_8" "input_noise_9"</b>


Модель також включила шуми.

Подивимося на навчену модель:

> summary(lm_m)

Call:
lm(formula = z ~ ., data = dat, model = T)

Residuals:
Min 1Q Median 3Q Max 
-0.060613 -0.022650 -0.000173 0.024939 0.048544 

Coefficients:
Estimate Std. Error t value Pr(>|t|) 
(Intercept) 0.0232453 0.0005581 41.651 < 2e-16 ***
<b>x -0.0257686 0.0052998 -4.862 1.17 e-06 ***
y -0.1572786 0.0052585 -29.910 < 2e-16 ***</b>
input_noise_2 -0.0017249 0.0027680 -0.623 0.533 
input_noise_5 -0.0027391 0.0027848 -0.984 0.325 
input_noise_6 0.0032417 0.0027907 1.162 0.245 
input_noise_8 0.0044998 0.0027723 1.623 0.105 
input_noise_9 0.0006839 0.0027808 0.246 0.806 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1

Residual standard error: 0.02742 on 17965 degrees of freedom
Multiple R-squared: 0.04937, Adjusted R-squared: 0.049 
F-statistic: 133.3 on 7 and 17965 DF, p-value: < 2.2 e-16


Аналогічно, всередині моделі бачимо, що важливі тільки оригінальні входи.

Якщо ж ми навчимо модель тільки на змінних X і Y, ми отримаємо:

> print(subset)
<b>[1] "x" "y"</b>
> print(r_sq_validate)
<b>[1] 0.05185492</b>

> summary(lm_m)

Call:
lm(formula = z ~ ., data = dat, model = T)

Residuals:
Min 1Q Median 3Q Max 
-0.059884 -0.022653 -0.000209 0.024955 0.048238 

Coefficients:
Estimate Std. Error t value Pr(>|t|) 
(Intercept) 0.0233808 0.0005129 45.590 < 2e-16 ***
<b>x -0.0257813 0.0052995 -4.865 1.15 e-06 ***
y -0.1573098 0.0052576 -29.920 < 2e-16 ***</b>
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1

Residual standard error: 0.02742 on 17970 degrees of freedom
Multiple R-squared: 0.04908, Adjusted R-squared: 0.04898 
F-statistic: 463.8 on 2 and 17970 DF, p-value: < 2.2 e-16


Справа в тому, що на валідації R^2 був вищий тоді, коли вимкнені шумові змінні.

Дивний результат! Напевно, через структури даних шуми не мають згубного впливу на модель.

Але ми ще не спробували врахувати взаємодію предикторів.

Код
lm_m <- lm(formula = z ~ x * y,
data = dat,
model = T)

> summary(lm_m)

Call:
lm(formula = z ~ x * y, data = dat, model = T)

Residuals:
Min 1Q Median 3Q Max 
-0.057761 -0.023067 -0.000119 0.024762 0.049747 

Coefficients:
Estimate Std. Error t value Pr(>|t|) 
(Intercept) 0.0196766 0.0006545 30.062 <2e-16 ***
x -0.1513484 0.0148113 -10.218 <2e-16 ***
y -0.1084295 0.0075183 -14.422 <2e-16 ***
x:y 1.3771299 0.1517363 9.076 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1

Residual standard error: 0.02736 on 17969 degrees of freedom
Multiple R-squared: 0.05342, Adjusted R-squared: 0.05327 
F-statistic: 338.1 on 3 and 17969 DF, p-value: < 2.2 e-16




Вийшло непогано. Спільну роботу вели X і Y значимо. А що щодо R^2 на валідації?

> lm_predict <- predict(lm_m,
+ newdata = sampleB)
> 1 - sum((sampleB$z - lm_predict) ^ 2) / sum((sampleB$z - mean(sampleB$z)) ^ 2)
<b>[1] 0.05464066</b>


Це найвище значення, з тих, що ми бачили. На жаль, саме варіант взаємодії не був закладений у фітнес-функцію і ми пропустили цю комбінацію входів.

Експеримент №2: жадібний пошук поднабора входів з лінійною функцією оцінки важливості (як фітнес-функції буде використовуватися варіант embedded — f-статистика навченої моделі на навчальній вибірці).

Код
linear_fit_f <- function(subset){
dat <- sampleA[, c(subset, "z")]

lm_m <- lm(formula = z ~.,
data = dat,
model = T)

print(subset)
print(summary(lm_m)$fstatistic[[1]])
return(summary(lm_m)$fstatistic[[1]])
}

#subset <- backward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = linear_fit_f)
subset <- forward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = linear_fit_f)



Результат послідовного включення змінних — тільки один предиктор Y. Для нього максимизировалась F-Statistic. Тобто це змінна дуже важлива. Але чомусь забута змінна X.

А тепер послідовне виключення змінних.

Результат аналогічний — тільки одна змінна Y.

Можна відзначити, що при максимізації F-Statistic многопеременной моделі всі шуми опинилися за бортом, і модель при цьому вийшла робастної: на валідації коефіцієнт детермінації майже не поступається найкращої моделі з експерименту №1:

> r_sq_validate
<b>[1] 0.05034534</b>


Експеримент №3: почергова оцінка індивідуальної значущості предикторів з допомогою коефіцієнта кореляції Пірсона (цей варіант найбільш простий, він не враховує взаємодії, фітнес-функція також проста — оцінює тільки лінійну зв'язок).

Код
correlation_arr <- data.frame()

for (i in 1:12){

correlation_arr[i, 1] <- colnames(sampleA)[i]
correlation_arr[i, 2] <- cor(sampleA [i], sampleA[, 'z'])

}



Результат:

> correlation_arr
V1 V2
<b>1 x 0.0413782832
2 y -0.2187061876</b>
3 input_noise_1 -0.0097719425
4 input_noise_2 -0.0019297383
5 input_noise_3 0.0002143946
6 input_noise_4 -0.0142325764
7 input_noise_5 -0.0048206943
8 input_noise_6 0.0090877674
9 input_noise_7 -0.0152897433
10 input_noise_8 0.0143477495
11 input_noise_9 0.0027560459
12 input_noise_10 -0.0079526578


Найбільша кореляція Z Y, і на другому місці — X. Однак кореляція X не явно виражена і вимагає проведення статистичного тесту на значущість відмінності коефіцієнта кореляції від нуля для кожної з змінних.

З іншого боку, у всіх проведених 3-х експериментах ми взагалі не враховували взаємодія предикторів (X * Y). Цим можна пояснити те, що оцінка одиничної значущості або лінійне включення предикторів рівняння не дає нам однозначної відповіді.

Такий экспериментик:

> cor(sampleA$x * sampleA$y, sampleA$z)
<b>[1] 0.1211382</b>


Показує, що взаємодія X і Y корелює з Z досить сильно.

Експеримент №4: оцінка важливості предикторів вбудованим в машину алгоритмом (варіант жадібного пошуку і embedded фітнес-функції важливості входів в GBM).

Навчимо Gradient Boosted Trees (gbm) і подивимося на важливості змінних. Хороша стаття, детально описує аспекти застосування GBM: Gradient boosting machines, a tutorial.

Візьмемо параметри навчання зі стелі, встановимо дуже низьку швидкість навчання, щоб уникнути сильного перенавчання. Зауважимо, що будь-які дерева рішень жадібні, а поліпшення моделі шляхом додавання багатьох моделей досягається семплюванням спостережень і входів.

Код
library(gbm)

gbm_dat <- bunny_dat[, c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

gbm_fit <- gbm(formula = z ~.,
distribution = "gaussian",
data = gbm_dat,
n.trees = 500,
interaction.depth = 12,
n.minobsinnode = 100,
shrinkage = 0.0001,
bag.fraction = 0.9,
train.fraction = 0.7,
n.cores = 6)

gbm.перфорація(object = gbm_fit, 
plot.it = TRUE, 
oobag.curve = F, 
overlay = TRUE)

summary(gbm_fit)



Результат:

> summary(gbm_fit)
var rel.inf
<b>y y 69.7919
x x 30.2081</b>
input_noise_1 input_noise_1 0.0000
input_noise_2 input_noise_2 0.0000
input_noise_3 input_noise_3 0.0000
input_noise_4 input_noise_4 0.0000
input_noise_5 input_noise_5 0.0000
input_noise_6 input_noise_6 0.0000
input_noise_7 input_noise_7 0.0000
input_noise_8 input_noise_8 0.0000
input_noise_9 input_noise_9 0.0000
input_noise_10 input_noise_10 0.0000


Цей підхід чудово впорався з поставленим завданням і виділив нешумовые входи, зробивши всі інші входи абсолютно незначущими.

Причому, треба зауважити, що параметри цього експерименту дуже швидка, все працює практично «з коробки». Більш ретельне планування даного експерименту, включаючи кроссвалидацию для отримання оптимальних параметрів навчання — справа більш складна, але при підготовці реальної моделі в production зробити це необхідно.

Експеримент №5: оцінка важливості предикторів, використовуючи стохастичний пошук з лінійною функцією оцінки важливості (нежадібних пошук у просторі входів, а як фітнес-функції буде використовуватися варіант wrapper — оцінка коефіцієнта детермінації навченої моделі на валідаційної вибірці).

В цей раз навчається лінійна модель включає попарні взаємодії між предикторами.

Код
### simulated annealing search with linear model interactions stats

library(scales)
library(GenSA)

sampleA <- bunny_dat[sample(nrow(bunny_dat), nrow(bunny_dat) / 2, replace = F), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

sampleB <- bunny_dat[!row.names(bunny_dat) %in% rownames(sampleA), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

#calculate parameters
predictor_number <- dim(sampleA)[2] - 1
sample_size <- dim(sampleA)[1]
par_v <- runif(predictor_number, min = 0, max = 1)
par_low <- rep(0, times = predictor_number)
par_upp <- rep(1, times = predictor_number)


############### the fitness function
sa_fit_f <- function(par){

indexes <- c(1:predictor_number)

for (i in 1:predictor_number){
if (par[i] >= threshold) {
indexes[i] <- i
} else {
indexes[i] <- 0
}
}

local_predictor_number <- 0
for (i in 1:predictor_number){
if (indexes[i] > 0) {
local_predictor_number <- local_predictor_number + 1
}
}

if (local_predictor_number > 0) {

sampleAf <- as.data.frame(sampleA[, c(indexes[], dim(sampleA)[2])])

lm_m <- lm(formula = z ~ .^2,
data = sampleAf,
model = T)

lm_predict <- predict(lm_m,
newdata = sampleB)

r_sq_validate <- 1 - sum((sampleB$z - lm_predict) ^ 2) / sum((sampleB$z - mean(sampleB$z)) ^ 2)

} else {
r_sq_validate <- 0
}

return(-r_sq_validate)
}


#stimating threshold for variable inclusion

threshold <- 0.5 # do it simply

#run feature selection

start <- Sys.time()

sao <- GenSA(par = par_v, fn = sa_fit_f, lower = par_low, upper = par_upp
, control = list(
#maxit = 10
max.time = 300
, smooth = F
simple.function = F))

trace_ff <- data.frame(sao$trace)$function.value
plot(trace_ff, type = "l")
percent(- sao$value)
final_vector <- c((sao$par >= threshold), T)
names(sampleA)[final_vector]
final_sample <- as.data.frame(sampleA[, final_vector])

Sys.time() - start


# check model
lm_m <- lm(formula = z ~ .^2,
data = sampleA[, final_vector],
model = T)

summary(lm_m)



Що ж вийшло?


<b>[1] "5.53%"</b>
> final_vector <- c((sao$par >= threshold), T)
> names(sampleA)[final_vector]
<b>[1] "x" "y" "input_noise_7" "input_noise_8" "input_noise_9" "z" </b>

> summary(lm_m)

Call:
lm(formula = z ~ .^2, data = sampleA[, final_vector], model = T)

Residuals:
Min 1Q Median 3Q Max 
-0.058691 -0.023202 -0.000276 0.024953 0.050618 

Coefficients:
Estimate Std. Error t value Pr(>|t|) 
(Intercept) 0.0197777 0.0007776 25.434 <2e-16 ***
<b>x -0.1547889 0.0154268 -10.034 <2e-16 ***
y -0.1148349 0.0085787 -13.386 <2e-16 ***</b>
input_noise_7 -0.0102894 0.0071871 -1.432 0.152 
input_noise_8 -0.0013928 0.0071508 -0.195 0.846 
input_noise_9 0.0026736 0.0071910 0.372 0.710 
<b>x:y 1.3098676 0.1515268 8.644 <2e-16 ***</b>
x:input_noise_7 0.0352997 0.0709842 0.497 0.619 
x:input_noise_8 0.0653103 0.0714883 0.914 0.361 
x:input_noise_9 0.0459939 0.0716704 0.642 0.521 
y:input_noise_7 0.0512392 0.0710949 0.721 0.471 
y:input_noise_8 0.0563148 0.0707809 0.796 0.426 
y:input_noise_9 -0.0085022 0.0710267 -0.120 0.905 
input_noise_7:input_noise_8 0.0129156 0.0374855 0.345 0.730 
input_noise_7:input_noise_9 0.0519535 0.0376869 1.379 0.168 
input_noise_8:input_noise_9 0.0128397 0.0379640 0.338 0.735 
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 '' 1

Residual standard error: 0.0274 on 17957 degrees of freedom
Multiple R-squared: 0.05356, Adjusted R-squared: 0.05277 
F-statistic: 67.75 on 15 and 17957 DF, p-value: < 2.2 e-16


Видно, що ми промахиваемся. Включені шуми.

Як бачимо, найкраще значення коефіцієнта детермінації на валідації досягнутий з включенням шумових змінних. При цьому збіжність алгоритму пошуку вичерпна:
image

Спробуємо поміняти фітнес-функцію, збережемо метод пошуку.

Експеримент №6: оцінка важливості предикторів, використовуючи стохастичний пошук з лінійною функцією оцінки важливості (нежадібних пошук у просторі входів, фітнес-функція — embedded p-values, відповідні коефіцієнтам моделі).

Ми будемо відбирати такий набір предикторів, у яких мінімізується середнє значення p-value входять в модель коефіцієнтів.

Код
# fittness based on p-values
sa_fit_f2 <- function(par){

indexes <- c(1:predictor_number)

for (i in 1:predictor_number){
if (par[i] >= threshold) {
indexes[i] <- i
} else {
indexes[i] <- 0
}
}

local_predictor_number <- 0
for (i in 1:predictor_number){
if (indexes[i] > 0) {
local_predictor_number <- local_predictor_number + 1
}
}

if (local_predictor_number > 0) {

sampleAf <- as.data.frame(sampleA[, c(indexes[], dim(sampleA)[2])])

lm_m <- lm(formula = z ~ .^2,
data = sampleAf,
model = T)

mean_val <- mean(summary(lm_m)$coefficients[, 4])

} else {
mean_val <- 1
}

return(mean_val)
}


#stimating threshold for variable inclusion

threshold <- 0.5 # do it simply

#run feature selection

start <- Sys.time()

sao <- GenSA(par = par_v, fn = sa_fit_f2, lower = par_low, upper = par_upp
, control = list(
#maxit = 10
max.time = 60
, smooth = F
simple.function = F))

trace_ff <- data.frame(sao$trace)$function.value
plot(trace_ff, type = "l")
percent(- sao$value)
final_vector <- c((sao$par >= threshold), T)
names(sampleA)[final_vector]
final_sample <- as.data.frame(sampleA[, final_vector])

Sys.time() - start



Результат:

> percent(- sao$value)
<b>[1] "-4.7 e-208%"</b>
> final_vector <- c((sao$par >= threshold), T)
> names(sampleA)[final_vector]
<b>[1] "y" "z"</b>


В цей раз все вийшло. Були відібрані тільки оригінальні предиктори, так як їх p-values дійсно малі.

Збіжність алгоритму хороша за 60 секунд:
image

Експеримент №7: оцінка важливості предикторів, використовуючи жадібний пошук з оцінкою важливості за якістю навченої моделі (це жадібний пошук в просторі входів, фітнес-функція — wrapper, відповідний коефіцієнту детермінації на валідації моделі boosted decision trees).

Код
### greedy search with wrapper GBM validation fitness

library(FSelector)
library(gbm)

sampleA <- bunny_dat[sample(nrow(bunny_dat), nrow(bunny_dat) / 2, replace = F), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

sampleB <- bunny_dat[!row.names(bunny_dat) %in% rownames(sampleA), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

gbm_fit <- function(subset){
dat <- sampleA[, c(subset, "z")]

gbm_fit <- gbm(formula = z ~.,
distribution = "gaussian",
data = dat,
n.trees = 50,
interaction.depth = 10,
n.minobsinnode = 100,
shrinkage = 0.1,
bag.fraction = 0.9,
train.fraction = 1,
n.cores = 7)

gbm_predict <- predict(gbm_fit,
newdata = sampleB,
n.trees = 50)

r_sq_validate <- 1 - sum((sampleB$z - gbm_predict) ^ 2) / sum((sampleB$z - mean(sampleB$z)) ^ 2)

print(subset)
print(r_sq_validate)
return(r_sq_validate)
}

subset <- backward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = gbm_fit)
subset <- forward.search(attributes = names(sampleA)[1:(ncol(sampleA) - 1)], eval.fun = gbm_fit)




Результат при жадібного включення предикторів:

> subset
<b>[1] "x" "y"</b>

> r_sq_validate
<b>[1] 0.2363794</b>


Потрапили в точку!

Результат при жадібного виключення предикторів:

> subset
<b> [1] "x" "y" "input_noise_1" "input_noise_2" "input_noise_3" "input_noise_4" "input_noise_5" "input_noise_6" "input_noise_7" 
[10] "input_noise_9" "input_noise_10"</b>

> r_sq_validate
<b>[1] 0.2266737</b>


Тут все гірше. Включення шумових предикторів в модель не сильно погіршило якість передбачення на валідації. І цьому є пояснення: випадкові лісу рішень мають вбудований регуляризатор і можуть самі ігнорувати неінформативні входи в процесі навчання.

На цьому ми завершимо розділ експериментів з ГІП на стандартних методах. І в наступному розділі я обосную і покажу практичне застосування, на мій погляд, статистично достовірного і робить свою роботу непогано, методу на основі інформаційних показників.

3) Використання теорії інформації для побудови фітнес-функції ГІП.


Ключове питання цього розділу — як описати поняття залежності та сформулювати його в інформаційно-теоретичному сенсі.


image

Источник.

Потрібно почати з поняття інформаційної ентропії. Ентропія (шенноновская) — це синонім невизначеності. Чим більше ми невпевнені щодо значення випадкової величини, тим більше ентропії (ще один синонім — інформації) несе реалізація даної величини. Якщо розглянути приклад з підкиданням монети, то найбільша ентропія буде у симетричної монети з усіх інших варіантів монет, тому що ми маємо найбільшу невизначеність у наступному виході підкидання.

Формула для шенноновской ентропії:

image

Що таке залежність?


Припустимо, що ми підкидаємо монету кілька разів. Чи можемо ми сказати, що наша невпевненість відносного наступного результату кидка зменшилася після того, як ми побачили результат попереднього кидка?

Припустимо, у нас є монета, яка приземляється орлом з імовірністю 2/3 після випадання решки і навпаки — після випадання орла вона приземляється на решку з імовірністю 2/3. При цьому безумовна частота випадінь орла і решки залишається 50/50.

Для такої монети після випадання орла частота випадання решки вже не дорівнює 1/2, також і навпаки. Так, наша невпевненість щодо наступного результату кидка зменшилася (ми вже не очікуємо 50/50).

Для розуміння феномена залежності згадаємо, що теорія імовірностей визначає незалежність так:
p(x,y) == p(x) * p(y)
таким чином, ймовірність спільної реалізації подій дорівнює добутку ймовірностей їх власних реалізацій.

Якщо це спостерігається, події незалежні в математичному сенсі. Якщо ж
p(x,y) != p(x) * p(y)

тоді події математично не незалежні.

Цей же принцип лежить в основі формули для вимірювання залежності між двома (і більше) випадковими величинами в теорії інформації.

Підкреслимо, що тут залежність розуміється в імовірнісному сенсі. Аналіз причинності вимагає більш комплексного огляду, що включає як аналіз помилкових кореляцій і надмірності (через застосування взаємної інформації), так і залучення експертного знання про досліджуваному об'єкті.

Взаємна інформація
image

Також можна вивести взаємну інформацію через ентропії:

image

Якщо сказати простими словами, взаємна інформація це кількість ентропії (невизначеності), яка йде з системи, при наявності предиктивної змінної (або декількох змінних). Наприклад, ентропія випадкової величини дорівнює 3 бітам. Взаємна інформація дорівнює 2 бітам. Значить, на 2/3 невизначеність щодо реалізації випадкової величини компенсується наявністю предиктора.

Взаємна інформація володіє наступними властивостями:

  • симетричність
  • може бути визначена для дискретних і неперервних змінних
  • звертається в нуль, якщо X і Y незалежні
  • нормалізована міра взаємної інформації приймає значення 1, якщо дві змінні повністю детермінують один одного


Взаємна інформація задовольняє деяким вимогам, що пред'являються до ідеальної міру залежності, сформульованим в:

Granger, C, E. Maasoumi e J. Racine, 'A Dependence Metric for Possibly Nonlinear Processes', Journal of Time Series Analysis 25, 2004, 649-669.


Є одна властивість взаємної інформації (ВИ), які корисно знати, застосовуючи дану міру:

  • ВІ може досягати свого максимуму, рівного меншим із значень ентропії для входу та виходу (предиктора і залежної змінної).


Це означає, що якщо у нас вхідна змінна має ентропією 10 біт, а вихідна — 3 біта, то максимум інформації, яка вхідна змінна може передавати вихідний дорівнює 3 бітам. Це максимум ВІ, якою може бути у такій системі.

Інший варіант — вхідна змінна має ентропією 3 біта, а вихід — 10 біт. Максимум інформації, яку вхід може нести на вихід — 3 біта, і це можливий максимум ВІ.

Якщо значення ВІ поділити на ентропію залежної змінної, вийде величина в діапазоні [0, 1], що, за аналогією з коефіцієнтом кореляції, що говорить про те, як, без урахування масштабу значень, вхідна змінна детермінує залежну.

Є ряд аргументів за те, що використання ВІ є кращим, хоча і пов'язане з рядом компромісів:

  • Метод дозволяє знаходити залежності довільного виду (лінійні та нелінійні) у просторі довільної розмірності;
  • Метод виключить всі вхідні змінні, якщо їх сумарна або індивідуальна інформативність нижче статистично значущою;
  • Слабкою стороною методу є той факт, що слабкі залежності — які ледь перевищують рівень статистичного шуму — будуть приховані від очей дослідника після застосування чисельно розрахованого квантиля метрики ВІ;
  • З невеликою кількістю значущих дискретних вхідних змінних дослідник в змозі побудувати людино-читається набір правил, що дозволяють дати інтерпретацію знайденим залежностям;
  • Немає необхідність складати рандомізовані комітети моделей (за аналогією з random forest і busting machines), так як метод у собі пробує безліч можливих рівнів «глибини дерева» для всіх відібраних предикторів та повертає найбільш значиму, у статистичному сенсі, комбінацію інформативних ознак;
  • Метод дуже дорогою, особливо порівняно з іншими методами, сочатающими жадібний градієнтний спуск і більш просту фітнес-функцію, що не вимагає симуляційних поправок;
  • Також варто зазначити, що застосування методу імітації відпалу для пошуку під наборів предикторів, як і інші методи стохастичного спуску, не гарантує знаходження глобального мінімуму за розумне число ітерацій.


Надмірність набору даних.


У наборі вхідних змінних може існувати взаємозалежність, інакше кажучи, залежність одного або декількох вхідних змінних від набору інших вхідних змінних. Якщо не усунути надмірність, модель потраплять зайві входи, які не несуть суттєвої інформації, що може призводити до перенавчання і ефекту коллинеарности.

Для квантування та вимірювання надлишковості є спеціальна метрика з теорії інформації — Мультиинформация (Multiinformation) або Тотальна кореляція (Total Correlation).

Джерело:

Її формула:

image

Дана міра була вперше представлена в:

Watanabe S (1960). Information theoretical analysis of multivariate correlation, IBM Journal of Research and Development 4, 66-82


Мультиинформация (МИ) приймає позитивне значення у випадку, коли ймовірність спільної реалізації всіх вхідних змінних починає відрізнятися від проведення їх маргінальних ймовірностей. Інтуїція цієї метрики в тому, що якщо будь підмножину змінних (від 1 до n) залежить від будь-якого іншого поднабора змінних (від 1 до m), то формула дає нам одним числом уявлення про силу цієї залежності. Максимум мультиинформации відповідає ситуації, коли який-небудь підмножину змінних повністю детермінований іншим поднабором.

Дана метрика корисна, коли потрібно оцінити надмірність вхідних даних у знайденому поднаборе.

Якщо говорити про заходи з теорії інформації в цілому, часто задають питання, чому Шеннон визначив ентропію саме через логарифм ймовірності. Він сам відповідає на питання наступним чином:

  • 1) Вона зручна практично. Параметри, важливі в інженерних додатках — такі, як час, пропускна спроможність, кількість перемикачів і так далі — зазвичай змінюються лінійно при логарифмічній зміні числа можливих варіантів. Приміром, додавання одного перемикача подвоює число можливих станів їх групи, збільшуючи на одиницю його логарифм за основою 2. Збільшення в два рази часу призводить до квадратичним зростання кількості повідомлень, або подвоєння їх логарифма, і так далі.
  • 2) Вона близька до нашого інтуїтивному уявленню про такій мірі. Це тісно пов'язано з попереднім пунктом, так як ми інтуїтивно вимірюємо величини, лінійно порівнюючи їх зі стандартами. Так, нам здається, що на двох перфокартах можна розмістити в два рази більше інформації, а по двом однаковим каналах — передати її в два рази більше.

  • 3) Вона зручна математично. Багато граничні переходи прості в логарифмах, в той час як в термінах числа варіантів вони досить нетривіальні.


Источник.

Для отримання більш повного предствлений про заходи з теорії інформації я запрошую читача звернутися, наприклад, до робіт:
Аківа Мойсейович Яглом, Ісаак Мойсейович Яглом. Вірогідність та інформація. Видання третє, перероблене і доповнене. М., Наука, 1973 — 512 с.
Колмогоров А. Н. — «Теорія інформації та теорія алгоритмів».


Поправка значення МІ і ВІ, отриманих на обмеженій вибірці.


Незважаючи на доступність пакетів і функцій, які без особливих зусиль порахують вам будь-яку метрику з теорії інформації, важливо розуміти, що ці метрики є наївним наближенням істинних значень в силу того, що вони отримані на емпіричному розподілі і вибіркою обмеженого розміру.

Інтуїція такого роду поправок полягає в тому, що знайдені залежності на нашій вибірці розміром n можуть скоректуватися при збільшенні вибірки до N.

Зростання розміру вибірки веде до зменшення залежності. Наведено дані, де апріорі не закладено ніяких залежностей.

Чисельна оцінка наївною взаємної інформації на випадковій вибірці з 100 спостережень.
image

Чисельна оцінка наївною взаємної інформації на випадковій вибірці з 1000 спостережень.
image

Зростання обсягу вибірки це головний фактор, що впливає на шумове кількість ВІ, яка прагне до нуля при гранично великій вибірці.

В ідеалі, ми хотіли б мати міру залежності, яка справедлива для вибірки нескінченного розміру. Це так звана ассимптотическая величина, оцінити яку бралися ряд вчених колективів з самого зародження теорії інформації. Див.:

Джерело 1.
Джерело 2.

Є ряд відомих назв поправок значень ВІ і МІ, а також ентропії, з урахуванням обмеженої вибірки. Їх корисність можна оцінити тільки шляхом проведення експериментів з даними. Ми намагалися оцінити ассимптотическое значення взаємної інформації для невеликої вибірки, не містить жодних залежностей апріорі (дані були синтетичними). При цьому ассимптотические поправки влаштовані так, що навіть на таких шумових даних вони починають показувати наявність залежностей.

Щоб подолати цю особливість інформаційних заходів є шлях набагато довший, але показує стійкі і обґрунтовані результати.

Ми використовуємо випадкове перемішування значень вихідної змінної при зафіксованої структурі вхідних даних. Таким чином ми повністю руйнуємо емпіричні залежності і на безлічі ітерацій цього процесу фіксуємо отримані значення ВІ. В результати є розподіл значень ВІ на даних, в яких відсутні залежності, але збережені розподілу входів та виходу. Ми оцінюємо заданий заздалегідь квантиль отриманого розподілу.

Таку ж операцію можна зробити для корекції значення мультиинформации (заходи надлишковості даних), але на цей раз випадкове перемішування зачіпає всі вхідні змінні і руйнуються зв'язки між будь-якими поднаборами змінних при збереженні власних розподілів змінних.

Проте один момент залишається невирішеним. Це ассимптотическая міра ентропії змінних. Не представляється можливим, без введення додаткових припущень, скорегувати маргінальний розподіл ймовірностей, дане емпірично. Іншими словами, якщо монета при 1 000 кидках випала орлом в 51% випадків, то без знання фізичних властивостей об'єкта складно щось сказати про довірчому інтервалі ймовірності на просторі результатів процесу. Тому в нашому підході ентропія залишається оціненої наївно. Але попрацювати над цим питанням ще можна в майбутньому.

Блок-схема роботи алгоритму на основі скоригованої взаємної інформації і стохастичного пошуку поднабора предикторів.


  • а) Привести набір даних до категориальному увазі одним з відомих способів.
  • б) Оцінити такі параметри набору даних, як кількість рядків, предикторів, середня кількість рівнів предикторів (якщо вони мають різну кількість рівнів) і на основі цих даних порахувати «оптимальна кількість» предикторів у фінальному поднаборе.
  • ) Ініціалізувати вектор типу numeric завдовжки кількість вхідних змінних з допомогою випадкових чисел, рівномірно розподілених в діапазоні [0, 1] і поставити нижню (0) і верхнуюю (1) кордон для значень вектора — це аргументи функції SA.
  • г) Ініціалізувати функції оцінки квантиля мультиинформации, квантиля взаємної інформації, фітнес-функцію, об'єднує всі розрахунки.
  • д) Встановити кількість симуляцій Монте-Карло для оцінки квантилів МІ і ВІ; задати квантили (наприклад, 0.9) для значень шумових МИ і ВИ.
  • е) Задати час або кількість ітерацій роботи алгоритму. Чим більше, тим краще.
  • ж) Запустити алгоритм, дочекатися результатів.


Пункт «б» вимагає пояснення. Оптимальна кількість змінних — це умовна величина, яка визначаться за формулою:

optim_var_num < — log(x = sample_size / 100, base = round(mean_levels, 0))


Інтуїція в тому, що, при допущенні незалежності вхідних змінних, середня кількість рівнів змінної потрібно звести в шукану ступінь з тим, щоб отримати загальну кількість унікальних взаємодіючих рівнів таке, щоб на кожному з них в середньому концентрувалося не менше n наблюдей, де n прийнято в експерименті за 100.

Ми не можемо включити занадто багато вхідних змінних з великою кількістю рівнів, так як консервативна оцінка частоти спостережень за вхідними рівнями вийде занадто маленька для того, щоб дати статистичне висновок про залежність вихідної змінної від набору вхідних.

При цьому, встановлюючи поріг, вище якого значення вектора завдовжки кількість вхідних змінних будуть перетворюватися в одиницю (прапор включення індексу змінної), ми робимо імовірнісний розрахунок:

threshold < — 1 — optim_var_num / predictor_number


Суть його зводиться до того, що ми визначаємо максимальне значення ймовірності для відбору обчисленого оптимальної кількості входів. І перевіряється ця логіка застосуванням біноміального розподілу.

Наприклад, візьмемо наші дані: половину всього набору, що використовується для навчання.

У нас 17 973 рядка, 12 предикторів, в кожному з 5 рівнів. В результаті застосування вищеописаних формул ми отримуємо, що оптимальна кількість предикторів одно 3,226.

Застосовуючи формулу порогу для включення предиктора в набір, отримуємо 0,731.

Яке найбільш ймовірне кількість відібраних змінних буде отримано на биномиальном розподілі?
image

Найбільш вірогідним є 3 спостереження. Якщо бути точним, то 5 ^ 3,226 дасть нам 178 рівнів, на яких в середньому розміститься по 100 спостережень.

Проведемо експеримент з використанням наведених міркувань.

Експеримент №8 оцінка важливості предикторів, використовуючи стохастичний (нежадібних) пошук з оцінкою важливості за допомогою взаємної інформації (це нежадібних пошук в просторі входів, фітнес-функція — скоригована на шум і надмірність взаємна інформація).

Експеримент буде запущений при симуляції 0,9 квантиля ВІ (на шумі) при 10 итерциях Монте-Карло, і квантиля 1 при 1 ітерації для коригування МІ (шум). Час поки поставимо 10 хвилин.

Код
############## simlated annelaing search with mutual information fitness function

library(infotheo) # measured in nats, converted to bits
library(scales)
library(GenSA)

sampleA <- bunny_dat[sample(nrow(bunny_dat), nrow(bunny_dat) / 2, replace = F), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

sampleB <- bunny_dat[!row.names(bunny_dat) %in% rownames(sampleA), c("x",
"y",
"input_noise_1", 
"input_noise_2",
"input_noise_3",
"input_noise_4",
"input_noise_5",
"input_noise_6", 
"input_noise_7",
"input_noise_8",
"input_noise_9",
"input_noise_10",
"z")]

# discretize all variables

dat <- sampleA

disc_levels <- 5

for (i in 1:13){

naming <- paste(names(dat[i]), 'discrete', sep = "_")
dat[, eval(naming)] <- discretize(dat[, eval(names(dat[i]))], disc = "equalfreq", nbins = disc_levels)[,1]

}

sampleA <- dat[, 14:26]

#calculate parameters
predictor_number <- dim(sampleA)[2] - 1
sample_size <- dim(sampleA)[1]
par_v <- runif(predictor_number, min = 0, max = 1)
par_low <- rep(0, times = predictor_number)
par_upp <- rep(1, times = predictor_number)


#load functions to memory
shuffle_f_inp <- function(x = data.frame(), iterations_inp, quantile_val_inp){

mutins <- c(1:iterations_inp)

for (count in 1:iterations_inp){

xx <- data.frame(1:dim(x)[1])

for (count1 in 1:(dim(x)[2] - 1)){

y <- as.data.frame(x[, count1])
y$count <- sample(1 : dim(x)[1], dim(x)[1], replace = F)
y <- y[order(y$count), ]

xx <- cbind(xx, y[, 1])

}

mutins[count] <- multiinformation(xx[, 2:dim(xx)[2]])

}

quantile(mutins, probs = quantile_val_inp)

}


shuffle_f <- function(x = data.frame(), iterations, quantile_val){

height <- dim(x)[1]

mutins <- c(1:iterations)

for (count in 1:iterations){

x$count <- sample(1 : height, height, replace = F)

y <- as.data.frame(c(x[dim(x)[2] - 1], x[dim(x)[2]]))
y <- y[order(y$count), ]

x[dim(x)[2]] <- NULL
x[dim(x)[2]] <- NULL
x$dep <- y[, 1]

rm(y)

receiver_entropy <- entropy(x[, dim(x)[2]])
received_inf <- mutinformation(x[, 1 : dim(x)[2] - 1], x[, dim(x)[2]])
corr_ff <- received_inf / receiver_entropy

mutins[count] <- corr_ff

}

quantile(mutins, probs = quantile_val)

}

############### the fitness function
fitness_f <- function(par){

indexes <- c(1:predictor_number)

for (i in 1:predictor_number){
if (par[i] >= threshold) {
indexes[i] <- i
} else {
indexes[i] <- 0
}
}

local_predictor_number <- 0
for (i in 1:predictor_number){
if (indexes[i] > 0) {
local_predictor_number <- local_predictor_number + 1
}
}

if (local_predictor_number > 1) {

sampleAf <- as.data.frame(sampleA[, c(indexes[], dim(sampleA)[2])])

pred_entrs <- c(1:local_predictor_number)

for (count in 1:local_predictor_number){

pred_entrs[count] <- entropy(sampleAf[count])

}

max_pred_ent <- sum(pred_entrs) - max(pred_entrs)

pred_multiinf <- multiinformation(sampleAf[, 1:dim(sampleAf)[2] - 1])

pred_multiinf <- pred_multiinf - shuffle_f_inp(sampleAf, iterations_inp, quantile_val_inp)

if (pred_multiinf < 0){

pred_multiinf <- 0

}

pred_mult_perc <- pred_multiinf / max_pred_ent

inf_corr_val <- shuffle_f(sampleAf, iterations, quantile_val)

receiver_entropy <- entropy(sampleAf[, dim(sampleAf)[2]])
received_inf <- mutinformation(sampleAf[, 1:local_predictor_number], sampleAf[, dim(sampleAf)[2]])

if (inf_corr_val - (received_inf / receiver_entropy) < 0){

fact_ff <- (inf_corr_val - (received_inf / receiver_entropy)) * (1 - pred_mult_perc)
} else {

fact_ff <- inf_corr_val - (received_inf / receiver_entropy)

}

} else if (local_predictor_number == 1) {

sampleAf<- as.data.frame(sampleA[, c(indexes[], dim(sampleA)[2])])

inf_corr_val <- shuffle_f(sampleAf, iterations, quantile_val)

receiver_entropy <- entropy(sampleAf[, dim(sampleAf)[2]])
received_inf <- mutinformation(sampleAf[, 1:local_predictor_number], sampleAf[, dim(sampleAf)[2]])

fact_ff <- inf_corr_val - (received_inf / receiver_entropy)

} else {
fact_ff <- 0
}
return(fact_ff)
}


########## estimating threshold for variable inclusion

iterations = 10
quantile_val = 0.9

iterations_inp = 1
quantile_val_inp = 1

levels_arr <- numeric()
for (i in 1:predictor_number){
levels_arr[i] <- length(unique(sampleA [i]))
}

mean_levels <- mean(levels_arr)

optim_var_num <- log(x = sample_size / 100, base = round(mean_levels, 0))

if (optim_var_num / predictor_number < 1){

threshold <- 1 - optim_var_num / predictor_number

} else {

threshold <- 0.5

}


#run feature selection

start <- Sys.time()

sao <- GenSA(par = par_v, fn = fitness_f, lower = par_low, upper = par_upp
, control = list(
#maxit = 10
max.time = 300
, smooth = F
simple.function = F))

trace_ff <- data.frame(sao$trace)$function.value
plot(trace_ff, type = "l")
percent(- sao$value)
final_vector <- c((sao$par >= threshold), T)
names(sampleA)[final_vector]
final_sample <- as.data.frame(sampleA[, final_vector])

Sys.time() - start



Результат:

> percent(- sao$value)
<b>[1] "18.1%"</b>
> final_vector <- c((sao$par >= threshold), T)
> names(sampleA)[final_vector]
<b>[1] "x_discrete" "y_discrete" "input_noise_2_discrete" "z_discrete" </b> 
> final_sample <- as.data.frame(sampleA[, final_vector])
> 
> Sys.time() - start
Time difference of 10.00453 mins


Майже вийшло. Однак збіжності за цей час алгоритм не показує. Багато даних — потрібно більше часу.

Спробуємо запустити на 20 хвилин, і при цьому зменшимо кількість ітерацій Монте-Карло до 1 для прискорення обрахунку.

Результат:

> percent(- sao$value)
<b><b>[1] "18.2%"</b></b>
> final_vector <- c((sao$par >= threshold), T)
> names(sampleA)[final_vector]
<b><b>[1] "x_discrete" "y_discrete" "input_noise_1_discrete" </b>"z_discrete" </b> 
> final_sample <- as.data.frame(sampleA[, final_vector])
> 
> Sys.time() - start
Time difference of 20.00585 mins


Також немає ідеального рішення — до хорошим входів доданий один шумовий предиктор.

При цьому отримана взаємна інформація на 18% компенсує ентропію залежної змінної. Непогане значення.

Збіжність непогана (графік майже повністю пішов на рівень хороших енергетичних значень):
image

Minimum-redundancy-maximum-relevance (mRMR)

Ідея максимизиации релевантності набору фичей та мінімізації їх надмірності не нова. Метод з відповідною назвою був представлений в:

«Feature selection based on mutual information: criteria of max-dependency, max-relevance, and min-redundancy,»
Hanchuan Peng, Fuhui Long and Chris Ding
IEEE Transactions on Pattern Analysis and Intelligence Machine,
Vol. 27, No. 8, pp.1226-1238, 2005


Незважаючи на раціональність методу і хороший привід для його впровадження в практику, у нас (а також статистичного спільноти) є ряд методологічних претензій. Источник.

Реалізація методу автором передбачає завмер релевантності як (а) усереднення взаємної інформації між кожним з відібраних предикторів і виходом, а завмер надлишковості як (б) усереднення взаємної інформації між усіма вхідними в підмножиною предикторами.

За пунктом (а) можна зробити законне зауваження про те, що вся інформація, передана від предикторів на вихід, не може бути коректно оцінена через просте середнє значення. Залишаються непоміченими ті залежності, які формуються на взаємодії предикторів та входу.

За пунктом (б) також має силу зауваження про те, що міра взаємозв'язку набору предикторів не може бути коректно виведена взяттям середнього значення. Як мінімум, залишаться без уваги під наборів залежності від під наборів (в поднаборе).

Як правильно зазначається у статті на Вікіпедії, ефективність даного алгоритму в ряді випадків під великим питанням, до того ж легко зрозуміти, що оптимальною ситуацією, де алгоритм працює добре, є відсутність взаємодій і повна (попарная і погрупная) незалежність у наборі предикторів, що майже ніколи не виконується. Так само як у випадку вимог до даних лінійної регресії, mRMR пред'являє солідний набір припущень, без яких результати його роботи будуть спотворені.

Пункт (в) критики даного підходу полягає в тому, що не зрозумілим залишається спосіб корекції наївної оцінки заходів залежності в умовах вибірки обмеженого розміру — так звана limited sampling bias.

На противагу наведеним слабких місць в mRMR ми перевели ідею максимізації релевантності та мінімізації надлишковості у площину коректної практичної реалізації.

Підіб'ємо підсумок наших експериментів і зведемо результати в одну таблицю.
image

Якщо підвести якісь узагальнені підсумки, то у нас вийшли гарні результати, застосовуючи важливість предикторів у моделі GBM, а також при мінімізації p-values в лінійної моделі.

Застосування нежадного пошуку спільно з метрикою взаємної інформації дало майже ідеальний результат, який включає один штучно введений шумовий предиктор. Проте розслідування показало, що включення цього предиктора дійсно максимізує метрику, незважаючи на те, що результат без цього предиктора майже такий же. Виходить, що деяку додаткову інформацію про вихід комбінація X-Y-Noise1 дійсно несе. Нам подобається цей алгоритм саме тим, що він дозволяє репортовать про важливість відібраних вхідних ознак на заздалегідь заданому рівні довірі (через оцінку квантиля), що відповідає духу роботи фахівців за статистикою в Align Technology. Цей аспект вносить набагато більше прозорості в результати, на відміну від, наприклад, коефіцієнтів відносної важливості в Gradient Boosted Trees (Джерело.).

Звичайно, ми можемо розрахувати статистичну значимість і для інших методів, наприклад, для коефіцієнтів лінійної моделі або для коефіцієнтів кореляції. Але ми звертаємо увагу читача на те, що, на відміну від згаданих заходів лінійної залежності, взаємна інформація охоплює довільні види залежності, чим і є більш привабливою в своїй універсальності.

Застосування коефіцієнта кореляції дало змішаний результат, де складно було вирішити, чи включати предиктор чи ні. Потрібні додаткові статистичні тести на значимість кореляції. При цьому на інших наборах, де може взагалі бути відсутнім лінійна зв'язок, даний метод не дасть доброго результату.

На цьому я закінчую статтю. Не всі методи були випробувані, не всі тонкощі розкриті, але загальне уявлення про експерименти з відбором інформативних ознак, я сподіваюся, у вас з'явилося. Якщо у вас будуть запитання та обговорення, я завжди постараюся відповісти в коментарях.

Код експериментів в одному файлі на github: Git

Використані джерела інформації:
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.