符号回归和遗传规划

背景

回归分析是一种常用的统计方法，用来分析自变量和因变量的线性相关关系，在线性回归分析中，变量间的关系形式是确定的，只需要对关系式的系数做出估计。

符号回归是另一种回归分析方法，不同之处在于，符号回归不对变量间的关系形式做出假设，因此，符号回归过程既包括函数形式/运算符的确定，也包括函数中变量系数的确定。符号回归相比于线性回归的优势在于可以挖掘变量间可能存在的非线性关系，但也可能过拟合。

遗传规划（Genetic Program）是一种求解符号回归的方法。本文使用python的gplearn包实现符号回归。

数据准备

本文内容主要参考[1]。假设现在有两个自变量X1,X2，与因变量的关系为

随机生成100个X0,X1,y的随机变量，50个做trainset，50个做testset，用遗传规划的方式来推测变量之间的关系式，作为对比，也用随机森林和决策树的算法做一遍看效果。

上述非线性的方程在X0,X1取值[-1,1]的情况下，图像如下，代码文件获取后台回复“gplearn”。

代码语言：javascript复制

from gplearn.genetic import SymbolicRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.tree import DecisionTreeRegressor
from sklearn.utils.random import check_random_state
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import graphviz
import pydotplus

# plot
x0 = np.arange(-1, 1, 1/10.)
x1 = np.arange(-1, 1, 1/10.)
x0, x1 = np.meshgrid(x0, x1)
y_truth = x0**2 - x1**2   x1 - 1

ax = plt.figure(figsize = (10,7)).gca(projection='3d')
ax.set_xlim(-1, 1)
ax.set_ylim(-1, 1)
surf = ax.plot_surface(x0, x1, y_truth, rstride=1, cstride=1,
                       color='green', alpha=0.5)
plt.show()

在曲面上进行随机数模拟生成训练集和测试集

代码语言：javascript复制

rng = check_random_state(0)

# Training samples
X_train = rng.uniform(-1, 1, 100).reshape(50, 2)
y_train = X_train[:, 0]**2 - X_train[:, 1]**2   X_train[:, 1] - 1

# Testing samples
X_test = rng.uniform(-1, 1, 100).reshape(50, 2)
y_test = X_test[:, 0]**2 - X_test[:, 1]**2   X_test[:, 1] - 1

gplearn 遗传规划

关于gplearn就不过多说明了，个人感觉它是python中快速实现遗传规划最好用的包，类似的包还有deap，但deap如果要做向量的模型，需要改很多代码，放弃了。

用遗传规划的方式实现符号回归总体来说是一个优化问题，首先，最优化的目标的是模型预测值和真实值之间的MSE或MAE或RMSE，跟机器学习算法是一样的。

其次，符号回归需要事先给出一些运算符号，gplearn中内置的运算符包括加减乘除、开方、对数、三角函数等，详细见下图，用户也可以自定义一些运算符

之后需要设定的是遗传规划相关的一些参数，比如适应度函数、变异概率等等，这里适应度函数其实就是上述的优化目标，默认是MAE，在一些特定情境下，也可以改成去最大化spearman相关系数等等，这里就不细加说明了。

最后，实现符号回归的函数是SymbolicRegressor，gplearn中主要有三个训练模型的函数：SymbolicRegressor、SymbolicClassifier、SymbolicTransformer，分别是符号回归、符号分类，还有符号转换，第三个函数主要是用遗传规划构建因子时候用的。

这里我们用的是SymbolicRegressor，用起来也很简单，把训练集丢进去，再把前面说的参数都设置好就可以了，符号函数直接用内置的也够了。代码如下

代码语言：javascript复制

# gp
est_gp = SymbolicRegressor(population_size=5000,
                           generations=20, 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)
est_gp.fit(X_train, y_train)

print(est_gp._program)

这里需要说明的是，gplearn中的print是被重写过的，print之后会输出最终的符号回归形式，上述代码运行后的结果如下

这里Fitness是适应度，默认MAE，随着模型迭代，最终趋近于0。由于模型很简单，所以很快就收敛了。最后一行输出了训练得到的最优表达式，mul是乘法，sub是减法，add是加法，结果来看，仅常数项差了0.001,效果非常好。

此外，gplearn还可以把训练得到的表达式计算过程用树的形式直观的展示出来，结果如下

代码语言：javascript复制

plt.figure(figsize = (10,7))
dot_data = est_gp._program.export_graphviz()
graph = graphviz.Source(dot_data)
graph

最后把符号回归和决策树、随机森林训练的结果做一个对比

代码语言：javascript复制

# 决策树、随机森林
est_tree = DecisionTreeRegressor()
est_tree.fit(X_train, y_train)
est_rf = RandomForestRegressor(n_estimators=10)
est_rf.fit(X_train, y_train)


y_gp = est_gp.predict(np.c_[x0.ravel(), x1.ravel()]).reshape(x0.shape)
score_gp = est_gp.score(X_test, y_test)
y_tree = est_tree.predict(np.c_[x0.ravel(), x1.ravel()]).reshape(x0.shape)
score_tree = est_tree.score(X_test, y_test)
y_rf = est_rf.predict(np.c_[x0.ravel(), x1.ravel()]).reshape(x0.shape)
score_rf = est_rf.score(X_test, y_test)



fig = plt.figure(figsize=(12, 10))

for i, (y, score, title) in enumerate([(y_truth, 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, projection='3d')
    ax.set_xlim(-1, 1)
    ax.set_ylim(-1, 1)
    ax.set_xticks(np.arange(-1, 1.01, .5))
    ax.set_yticks(np.arange(-1, 1.01, .5))
    surf = ax.plot_surface(x0, x1, y, rstride=1, cstride=1, color='green', alpha=0.5)
    points = ax.scatter(X_train[:, 0], X_train[:, 1], y_train)
    if score is not None:
        score = ax.text(-.7, 1, .2, "$R^2 =/ %.6f$" % score, 'x', fontsize=14)
    plt.title(title)

plt.show()

用训练的模型去预测测试集的y，预测的y作图结果如下，符号回归完胜，但也并不是说符号回归一定是最好的，还是要取决于具体的问题场景。