当前位置：首页 > news >正文

Python遗传规划实战：用gplearn和DEAP解决符号回归问题（附完整代码）

news 2026/5/14 15:36:33

Python遗传规划实战：用gplearn和DEAP解决符号回归问题

符号回归是一种强大的机器学习技术，它能够从数据中自动发现隐藏的数学表达式。与传统的回归方法不同，符号回归不需要预先指定模型形式，而是通过进化计算自动探索可能的数学表达式空间。本文将深入探讨如何使用Python中的gplearn和DEAP库来实现符号回归，并提供完整的实战代码。

1. 符号回归与遗传规划基础

符号回归的核心思想是通过进化算法自动发现描述数据关系的数学表达式。这种方法特别适用于那些我们无法预先确定合适模型形式的问题场景。遗传规划(GP)作为遗传算法的一种扩展，专门用于解决这类符号回归问题。

在遗传规划中，每个候选解都是一个数学表达式树，树的节点可以是运算符(如加、减、乘、除)，叶子节点则是变量或常数。进化过程通过选择、交叉和变异等操作不断改进这些表达式树。

遗传规划的主要优势：

能够发现人类难以想到的复杂非线性关系
结果具有可解释性，不像神经网络那样的"黑箱"
不需要特征工程，自动处理变量间交互作用

提示：符号回归特别适合物理、工程领域的问题，因为发现的数学表达式往往具有物理意义。

2. gplearn库入门与实践

gplearn是一个专门为符号回归设计的Python库，它提供了与scikit-learn兼容的API，使得使用方式非常直观。下面我们通过一个完整示例来演示如何使用gplearn进行符号回归。

2.1 安装与环境准备

首先安装必要的库：

pip install gplearn numpy matplotlib scikit-learn

2.2 基础符号回归示例

我们首先生成一个简单的合成数据集来演示基本用法：

import numpy as np from gplearn.genetic import SymbolicRegressor # 生成合成数据 rng = np.random.RandomState(0) X_train = rng.uniform(-1, 1, 100).reshape(50, 2) y_train = X_train[:, 0]**2 - X_train[:, 1]**2 + X_train[:, 1] - 1 # 创建并训练符号回归模型 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, random_state=0) est_gp.fit(X_train, y_train) # 输出最佳程序 print(est_gp._program)

2.3 关键参数解析

gplearn提供了丰富的参数来控制进化过程：

种群参数：

population_size: 种群大小(默认500)
generations: 进化代数(默认20)

遗传操作概率：

p_crossover: 交叉概率(默认0.9)
p_subtree_mutation: 子树变异概率(默认0.01)
p_hoist_mutation: 提升变异概率(默认0.01)
p_point_mutation: 点变异概率(默认0.01)

其他重要参数：

function_set: 可用的函数集合(默认['add','sub','mul','div'])
metric: 适应度指标(默认'mean absolute error')
parsimony_coefficient: 控制程序复杂度的惩罚系数

2.4 自定义函数与适应度

gplearn允许我们自定义函数和适应度指标：

from gplearn.functions import make_function from gplearn.fitness import make_fitness # 自定义函数 def _protected_log(x): return np.log(np.abs(x)) protected_log = make_function(function=_protected_log, name='log', arity=1) # 自定义适应度指标 def _mape(y, y_pred, w): diffs = np.abs((y - y_pred) / np.maximum(0.001, y)) return 100. * np.average(diffs, weights=w) mape = make_fitness(function=_mape, greater_is_better=False) # 使用自定义函数和适应度 est_gp = SymbolicRegressor(function_set=['add', 'sub', 'mul', 'div', protected_log], metric=mape)

3. DEAP库高级遗传规划

DEAP(分布式进化算法框架)是一个更通用的进化计算框架，提供了更大的灵活性。下面介绍如何使用DEAP实现符号回归。

3.1 DEAP基础设置

首先设置基本类型和原始集合：

from deap import base, creator, tools, gp import operator import numpy as np # 创建适应度和个体类型 creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin) # 定义原始集合 pset = gp.PrimitiveSet("MAIN", arity=2) # 2个输入变量 pset.addPrimitive(operator.add, 2) pset.addPrimitive(operator.sub, 2) pset.addPrimitive(operator.mul, 2) pset.addPrimitive(operator.neg, 1) pset.addEphemeralConstant("rand101", lambda: np.random.uniform(-1,1))

3.2 工具箱配置

配置DEAP工具箱来注册各种进化操作：

toolbox = base.Toolbox() toolbox.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=3) toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr) toolbox.register("population", tools.initRepeat, list, toolbox.individual) toolbox.register("compile", gp.compile, pset=pset) # 定义评估函数 def evaluate(individual, points): func = toolbox.compile(expr=individual) sqerrors = [(func(x[0], x[1]) - y)**2 for x, y in points] return (np.mean(sqerrors),) toolbox.register("evaluate", evaluate, points=list(zip(X_train, y_train))) toolbox.register("select", tools.selTournament, tournsize=3) toolbox.register("mate", gp.cxOnePoint) toolbox.register("expr_mut", gp.genFull, min_=0, max_=2) toolbox.register("mutate", gp.mutUniform, expr=toolbox.expr_mut, pset=pset)

3.3 进化主循环

实现完整的进化过程：

def main(): pop = toolbox.population(n=300) hof = tools.HallOfFame(1) stats_fit = tools.Statistics(lambda ind: ind.fitness.values) stats_size = tools.Statistics(len) mstats = tools.MultiStatistics(fitness=stats_fit, size=stats_size) mstats.register("avg", np.mean) mstats.register("std", np.std) mstats.register("min", np.min) mstats.register("max", np.max) pop, log = algorithms.eaSimple(pop, toolbox, 0.5, 0.1, 40, stats=mstats, halloffame=hof, verbose=True) return pop, log, hof pop, log, hof = main() best_ind = hof[0] print(f"Best individual: {best_ind}") print(f"Fitness: {best_ind.fitness.values[0]}")

4. 实战案例与性能优化

4.1 真实数据集应用

让我们将符号回归应用于一个真实数据集 - 波士顿房价数据集：

from sklearn.datasets import load_boston from sklearn.model_selection import train_test_split # 加载数据 data = load_boston() X, y = data.data, data.target X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42) # 配置符号回归器 est_gp = SymbolicRegressor(population_size=1000, generations=20, function_set=['add', 'sub', 'mul', 'div', 'sqrt', 'log', 'abs'], metric='mean absolute error', parsimony_coefficient=0.01, random_state=0) # 训练和评估 est_gp.fit(X_train, y_train) print("Test score:", est_gp.score(X_test, y_test)) print("Best program:\n", est_gp._program)

4.2 性能优化技巧

1. 并行化计算：

from multiprocessing import Pool # 在评估函数中使用多进程 pool = Pool() toolbox.register("map", pool.map) # 运行进化算法 pop, log, hof = main() # 结束后关闭池 pool.close()

2. 早停机制：

from deap import tools # 在工具箱中添加早停条件 toolbox.decorate("evaluate", tools.DeltaPenalty(check_valid, 1e6)) def check_valid(individual): try: func = toolbox.compile(expr=individual) func(0, 0) # 测试评估 return True except: return False

3. 复杂问题分解：对于高维问题，可以尝试以下策略：

先进行特征选择或降维
使用多个GP模型分别建模，然后集成
分阶段进化：先简单函数集，后逐渐增加复杂度

4.3 结果可视化与分析

可视化进化过程和最佳个体：

import matplotlib.pyplot as plt # 绘制适应度变化 gen = log.select("gen") fit_mins = log.chapters["fitness"].select("min") size_avgs = log.chapters["size"].select("avg") fig, ax1 = plt.subplots() ax1.plot(gen, fit_mins, "b-", label="Minimum Fitness") ax1.set_xlabel("Generation") ax1.set_ylabel("Fitness", color="b") ax2 = ax1.twinx() ax2.plot(gen, size_avgs, "r-", label="Average Size") ax2.set_ylabel("Size", color="r") plt.title("Fitness and Size over Generations") plt.show() # 绘制最佳个体表达式树 nodes, edges, labels = gp.graph(best_ind) import networkx as nx g = nx.Graph() g.add_edges_from(edges) pos = nx.nx_agraph.graphviz_layout(g, prog="dot") nx.draw_networkx_nodes(g, pos, node_size=900, node_color='w') nx.draw_networkx_edges(g, pos) nx.draw_networkx_labels(g, pos, labels) plt.show()

5. 高级主题与扩展应用

5.1 多目标遗传规划

DEAP支持多目标优化，我们可以同时优化模型的准确性和复杂度：

from deap import algorithms, base, creator, tools # 多目标适应度 (最小化误差，最小化复杂度) creator.create("FitnessMulti", base.Fitness, weights=(-1.0, -1.0)) creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMulti) # 修改评估函数 def evaluate_mo(individual, points): func = toolbox.compile(expr=individual) sqerrors = [(func(x[0], x[1]) - y)**2 for x, y in points] return np.mean(sqerrors), len(individual) toolbox.register("evaluate", evaluate_mo, points=list(zip(X_train, y_train))) toolbox.register("select", tools.selNSGA2) # 运行多目标算法 pop = toolbox.population(n=100) hof = tools.ParetoFront() stats = tools.Statistics(lambda ind: ind.fitness.values) stats.register("avg", np.mean, axis=0) stats.register("std", np.std, axis=0) stats.register("min", np.min, axis=0) algorithms.eaMuPlusLambda(pop, toolbox, mu=50, lambda_=100, cxpb=0.7, mutpb=0.3, ngen=40, stats=stats, halloffame=hof)

5.2 符号分类问题

gplearn也可以用于分类问题：

from gplearn.genetic import SymbolicClassifier from sklearn.datasets import load_breast_cancer from sklearn.model_selection import train_test_split # 加载数据 data = load_breast_cancer() X, y = data.data, data.target X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42) # 创建符号分类器 est_clf = SymbolicClassifier(population_size=500, generations=20, function_set=['add', 'sub', 'mul', 'div', 'sqrt', 'log'], metric='log loss', parsimony_coefficient=0.01, random_state=0) # 训练和评估 est_clf.fit(X_train, y_train) print("Test accuracy:", est_clf.score(X_test, y_test))

5.3 符号Transformer特征工程

gplearn的SymbolicTransformer可以用于自动化特征工程：

from gplearn.genetic import SymbolicTransformer from sklearn.pipeline import Pipeline from sklearn.linear_model import Ridge # 创建符号特征转换器 transformer = SymbolicTransformer(generations=20, population_size=200, hall_of_fame=100, n_components=10, function_set=['add', 'sub', 'mul', 'div', 'sqrt'], parsimony_coefficient=0.0005, max_samples=0.9, random_state=0) # 构建管道 model = Pipeline([ ('symbolic', transformer), ('regressor', Ridge()) ]) # 训练和评估 model.fit(X_train, y_train) print("Test score:", model.score(X_test, y_test))

在实际项目中，我发现将符号回归与其他机器学习方法结合使用往往能取得更好的效果。例如，先用符号回归生成一些特征，再输入到神经网络或梯度提升树中。这种方法结合了符号回归的可解释性和传统机器学习模型的强大预测能力。

查看全文

http://www.jsqmd.com/news/516022/