В этом уроке я хочу познакомить вас с генетическим программированием на Python с помощью библиотеки gplearn.

Символьная регрессия - это метод машинного обучения, целью которого является определение основного математического выражения, которое лучше всего описывает отношения. Он начинается с создания совокупности наивных случайных формул для представления взаимосвязи между известными независимыми переменными и их целевыми зависимыми переменными с целью прогнозирования новых данных. Каждое последующее поколение программ затем эволюционирует из предыдущего, путем отбора наиболее приспособленных людей из популяции для проведения генетических операций.

1. Введение

Представьте, что вы ученый, работающий в любой области. Вы можете разделить свою работу на три этапа: первый - собрать данные, второй - предложить феноменологическое объяснение этих данных. В-третьих, вы хотели бы объяснить эти данные и феноменологические наблюдения основным принципам.
Представьте, что вы собираетесь собрать данные и открыть для себя гравитацию. В этой области нам следует помнить трех ученых: Тихо Браге (сбор данных), который провел обширные и точные измерения положения планет с течением времени.
Иоганн Кеплер, который на основе измерений Браге вывел аналитические выражения, кратко описывающие движение Солнечной системы (анализ данных).
Наконец, теоретик Исаак Ньютон. Он осознал механизмы, лежащие в основе движения планет вокруг Солнца, которые можно сформулировать в виде универсального закона (происходящего из первых принципов).

Модели машинного обучения (ML) в настоящее время являются предпочтительными инструментами для раскрытия этих физических законов. Несмотря на то, что они показали некоторые многообещающие результаты в прогнозировании свойств материалов, типичные параметризованные модели машинного обучения не позволяют выполнять конечные этапы научного исследования: автоматизацию Ньютона. Примите во внимание, например, алгоритмы машинного обучения / глубокого обучения, которые могут с идеальной оценкой предсказать количество инфицированных людей в глобальной пандемии: если они не могут объяснить количество, они имеют ограниченное применение. Поэтому до сих пор широко используются более простые механистические модели, такие как SIR. Модели машинного обучения могут быть прогнозирующими, но их описания часто слишком подробны (например, модели глубокого обучения с тысячами параметров) или математически ограничены (например, предполагается, что целевая переменная представляет собой линейную комбинацию входных функций).

В этой перспективной статье мы сосредоточимся на альтернативе моделям машинного обучения: символической регрессии. Символьная регрессия одновременно ищет оптимальную форму функции и набор параметров для данной проблемы и является мощным методом регрессии, когда доступно мало, если вообще какое-либо априорное знание структуры / распределения данных.

2. Исходный набор данных и анализ данных

Мы генерируем данные, как в учебнике по регрессии. Единственная разница здесь в том, что мы изменяем количество выборки и коэффициенты, чтобы получить сложную функцию. Мы видим, что результатом является комбинация линейного x, sin (x) и x ** 3.

nsample = 400
sig = 0.2
x = np.linspace(-50, 50, nsample)
X = np.column_stack((x/5, 10*np.sin(x), (x-5)**3, np.ones(nsample)))
beta = [0.01, 1, 0.001, 5.]
y_true = np.dot(X, beta)
y = y_true + sig * np.random.normal(size=nsample)
df = pd.DataFrame()
df['x']=x
df['y']=y

Форма кривой показана ниже:

3. Импорт и реализация GPlearn

Мы импортируем SymbolicRegressor из gplearn, а также дерево решений и регрессор случайного леса из sklearn, из которого мы будем сравнивать результаты.

Мы импортируем из sympy, в частности, мы будем использовать simpify, чтобы сделать выражение результата более читабельным.

Поскольку мы будем сравнивать результаты разных подходов, мы разбиваем наш набор данных на обучение / тест:

X = df[['x']]
y = df['y']
y_true = y
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30, random_state=42)

Наша первая подгонка будет чрезвычайно простой: в наборе функций мы импортируем большую часть, если не все, функции по умолчанию, которые поставляются с библиотекой. Затем мы инициализируем регрессор с некоторыми стандартными параметрами.

С помощью этого кода мы говорим регрессору использовать функцию из списка function_set. Число поколений (эволюции) выражения будет 40, с начальной популяцией 5000. Каждое из этих 5000 будет объединяться между собой, чтобы сформировать выражения, которые будут оцениваться. Функция оценки по умолчанию будет оценивать пригодность, используя по умолчанию метрику MSE. Существует критерий остановки, который по умолчанию будет равен 0,01 в случае, если мы достигнем этой MSE раньше, чем количество поколений. Затем есть три параметра, связанных с эволюцией: p_crossover, p_subtree_mutation и p_hoist_mutation. Подробнее о конкретном использовании этих параметров см. Документацию. Однако вы должны отметить одну вещь: сумма трех параметров должна быть равна или меньше 1. Мы также можем попытаться контролировать общую длину нашего выражения с помощью Parsimony_coefficient. Чем выше будет этот параметр, тем короче регрессор попытается сохранить выражение. С помощью max_samples мы говорим функции фитнеса, чтобы она сравнивала результат с 90% наших данных.

# First Test
function_set = ['add', 'sub', 'mul', 'div','cos','sin','neg','inv']
est_gp = SymbolicRegressor(population_size=5000,function_set=function_set,
                           generations=40, stopping_criteria=0.01,
                           p_crossover=0.7, p_subtree_mutation=0.1,
                           p_hoist_mutation=0.05, p_point_mutation=0.1,
                           max_samples=0.9, verbose=1,
                           parsimony_coefficient=0.01, random_state=0,
                          feature_names=X_train.columns)

Мы определяем словарь, называемый преобразователем: его использование будет ясным позже, но, по сути, он нам нужен для преобразования имени функции или строки в эквивалентный математический код Python.

converter = {
    'sub': lambda x, y : x - y,
    'div': lambda x, y : x/y,
    'mul': lambda x, y : x*y,
    'add': lambda x, y : x + y,
    'neg': lambda x    : -x,
    'pow': lambda x, y : x**y,
    'sin': lambda x    : sin(x),
    'cos': lambda x    : cos(x),
    'inv': lambda x: 1/x,
    'sqrt': lambda x: x**0.5,
    'pow3': lambda x: x**3
}

Затем мы подбираем данные, оцениваем нашу оценку R2 и с помощью simpyfy выводим полученное выражение:

est_gp.fit(X_train, y_train)
print('R2:',est_gp.score(X_test,y_test))
next_e = sympify((est_gp._program), locals=converter)
next_e

При verbose = 1 вывод будет примерно таким: в первом столбце указан номер поколения. Затем у нас есть средняя длина и приспособленность всей популяции (5000 единиц). Справа у нас лучшие параметры для человека.

С приведенным выше кодом результат будет примерно таким:

4. Сравнение Gplearn с традиционными подходами к машинному обучению.

Мы достигли хорошего уровня корреляции с исходным выражением. Однако должен сказать, что интерпретировать это уравнение особенно сложно! Давайте посмотрим, как это выглядит в сравнении со стандартными инструментами машинного обучения:

est_tree = DecisionTreeRegressor(max_depth=5)
est_tree.fit(X_train, y_train)
est_rf = RandomForestRegressor(n_estimators=100,max_depth=5)
est_rf.fit(X_train, y_train)
y_gp = est_gp.predict(X_test)
score_gp = est_gp.score(X_test, y_test)
y_tree = est_tree.predict(X_test)
score_tree = est_tree.score(X_test, y_test)
y_rf = est_rf.predict(X_test)
score_rf = est_rf.score(X_test, y_test)

Судя по приведенным выше оценкам, результат явно не впечатляет: оба инструмента машинного обучения показали довольно хорошие результаты по сравнению с нашим аналитическим выражением:

RF score:0.993
DT score: 0.992
GPlearn score: 0.978

Однако, если мы построим график результатов, мы можем выделить, что наша аналитическая функция на самом деле лучше в среднем диапазоне, в то время как она срабатывает в крайних точках.

fig = plt.figure(figsize=(12, 10))
for i, (y, score, title) in enumerate([(y_test, None, "Ground Truth"),
                                       (y_gp, score_gp, "SymbolicRegressor"),
                                       (y_tree, score_tree, "DecisionTreeRegressor"),
                                       (y_rf, score_rf, "RandomForestRegressor")]):
ax = fig.add_subplot(2, 2, i+1)
    points = ax.scatter(X, y_true, color='green', alpha=0.5)
    test = ax.scatter(X_test,y,color='red', alpha=0.5)
    plt.title(title)
plt.show()

5. Улучшение результата

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

from gplearn.functions import make_function
def pow_3(x1):
    f = x1**3
    return f
pow_3 = make_function(function=pow_3,name='pow3',arity=1)
# add the new function to the function_set
function_set = ['add', 'sub', 'mul', 'div','cos','sin','neg','inv',pow_3]
est_gp = SymbolicRegressor(population_size=5000,function_set=function_set,
                           generations=45, stopping_criteria=0.01,
                           p_crossover=0.7, p_subtree_mutation=0.1,
                           p_hoist_mutation=0.05, p_point_mutation=0.1,
                           max_samples=0.9, verbose=1,
                           parsimony_coefficient=0.01, random_state=0,
                          feature_names=X_train.columns)
est_gp.fit(X_train, y_train)
print('R2:',est_gp.score(X_test,y_test))
next_e = sympify((est_gp._program), locals=converter)
next_e

R2 теперь действительно близко к 1, также обратите внимание, насколько проще наше уравнение! Если мы снова сравним с традиционными инструментами машинного обучения, мы действительно ценим улучшения, сделанные на этом этапе.

6. Заключение

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