Python遗传规划实战用gplearn和DEAP解决符号回归问题符号回归是一种强大的机器学习技术它能够从数据中自动发现隐藏的数学表达式。与传统的回归方法不同符号回归不需要预先指定模型形式而是通过进化计算自动探索可能的数学表达式空间。本文将深入探讨如何使用Python中的gplearn和DEAP库来实现符号回归并提供完整的实战代码。1. 符号回归与遗传规划基础符号回归的核心思想是通过进化算法自动发现描述数据关系的数学表达式。这种方法特别适用于那些我们无法预先确定合适模型形式的问题场景。遗传规划(GP)作为遗传算法的一种扩展专门用于解决这类符号回归问题。在遗传规划中每个候选解都是一个数学表达式树树的节点可以是运算符(如加、减、乘、除)叶子节点则是变量或常数。进化过程通过选择、交叉和变异等操作不断改进这些表达式树。遗传规划的主要优势能够发现人类难以想到的复杂非线性关系结果具有可解释性不像神经网络那样的黑箱不需要特征工程自动处理变量间交互作用提示符号回归特别适合物理、工程领域的问题因为发现的数学表达式往往具有物理意义。2. gplearn库入门与实践gplearn是一个专门为符号回归设计的Python库它提供了与scikit-learn兼容的API使得使用方式非常直观。下面我们通过一个完整示例来演示如何使用gplearn进行符号回归。2.1 安装与环境准备首先安装必要的库pip install gplearn numpy matplotlib scikit-learn2.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_size5000, generations20, stopping_criteria0.01, p_crossover0.7, p_subtree_mutation0.1, p_hoist_mutation0.05, p_point_mutation0.1, max_samples0.9, verbose1, random_state0) 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, namelog, arity1) # 自定义适应度指标 def _mape(y, y_pred, w): diffs np.abs((y - y_pred) / np.maximum(0.001, y)) return 100. * np.average(diffs, weightsw) mape make_fitness(function_mape, greater_is_betterFalse) # 使用自定义函数和适应度 est_gp SymbolicRegressor(function_set[add, sub, mul, div, protected_log], metricmape)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, fitnesscreator.FitnessMin) # 定义原始集合 pset gp.PrimitiveSet(MAIN, arity2) # 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, psetpset, 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, psetpset) # 定义评估函数 def evaluate(individual, points): func toolbox.compile(exprindividual) sqerrors [(func(x[0], x[1]) - y)**2 for x, y in points] return (np.mean(sqerrors),) toolbox.register(evaluate, evaluate, pointslist(zip(X_train, y_train))) toolbox.register(select, tools.selTournament, tournsize3) toolbox.register(mate, gp.cxOnePoint) toolbox.register(expr_mut, gp.genFull, min_0, max_2) toolbox.register(mutate, gp.mutUniform, exprtoolbox.expr_mut, psetpset)3.3 进化主循环实现完整的进化过程def main(): pop toolbox.population(n300) hof tools.HallOfFame(1) stats_fit tools.Statistics(lambda ind: ind.fitness.values) stats_size tools.Statistics(len) mstats tools.MultiStatistics(fitnessstats_fit, sizestats_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, statsmstats, halloffamehof, verboseTrue) return pop, log, hof pop, log, hof main() best_ind hof[0] print(fBest individual: {best_ind}) print(fFitness: {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_size0.2, random_state42) # 配置符号回归器 est_gp SymbolicRegressor(population_size1000, generations20, function_set[add, sub, mul, div, sqrt, log, abs], metricmean absolute error, parsimony_coefficient0.01, random_state0) # 训练和评估 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(exprindividual) func(0, 0) # 测试评估 return True except: return False3. 复杂问题分解 对于高维问题可以尝试以下策略先进行特征选择或降维使用多个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-, labelMinimum Fitness) ax1.set_xlabel(Generation) ax1.set_ylabel(Fitness, colorb) ax2 ax1.twinx() ax2.plot(gen, size_avgs, r-, labelAverage Size) ax2.set_ylabel(Size, colorr) 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, progdot) nx.draw_networkx_nodes(g, pos, node_size900, node_colorw) 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, fitnesscreator.FitnessMulti) # 修改评估函数 def evaluate_mo(individual, points): func toolbox.compile(exprindividual) sqerrors [(func(x[0], x[1]) - y)**2 for x, y in points] return np.mean(sqerrors), len(individual) toolbox.register(evaluate, evaluate_mo, pointslist(zip(X_train, y_train))) toolbox.register(select, tools.selNSGA2) # 运行多目标算法 pop toolbox.population(n100) hof tools.ParetoFront() stats tools.Statistics(lambda ind: ind.fitness.values) stats.register(avg, np.mean, axis0) stats.register(std, np.std, axis0) stats.register(min, np.min, axis0) algorithms.eaMuPlusLambda(pop, toolbox, mu50, lambda_100, cxpb0.7, mutpb0.3, ngen40, statsstats, halloffamehof)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_size0.2, random_state42) # 创建符号分类器 est_clf SymbolicClassifier(population_size500, generations20, function_set[add, sub, mul, div, sqrt, log], metriclog loss, parsimony_coefficient0.01, random_state0) # 训练和评估 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(generations20, population_size200, hall_of_fame100, n_components10, function_set[add, sub, mul, div, sqrt], parsimony_coefficient0.0005, max_samples0.9, random_state0) # 构建管道 model Pipeline([ (symbolic, transformer), (regressor, Ridge()) ]) # 训练和评估 model.fit(X_train, y_train) print(Test score:, model.score(X_test, y_test))在实际项目中我发现将符号回归与其他机器学习方法结合使用往往能取得更好的效果。例如先用符号回归生成一些特征再输入到神经网络或梯度提升树中。这种方法结合了符号回归的可解释性和传统机器学习模型的强大预测能力。
Python遗传规划实战：用gplearn和DEAP解决符号回归问题（附完整代码）
Python遗传规划实战用gplearn和DEAP解决符号回归问题符号回归是一种强大的机器学习技术它能够从数据中自动发现隐藏的数学表达式。与传统的回归方法不同符号回归不需要预先指定模型形式而是通过进化计算自动探索可能的数学表达式空间。本文将深入探讨如何使用Python中的gplearn和DEAP库来实现符号回归并提供完整的实战代码。1. 符号回归与遗传规划基础符号回归的核心思想是通过进化算法自动发现描述数据关系的数学表达式。这种方法特别适用于那些我们无法预先确定合适模型形式的问题场景。遗传规划(GP)作为遗传算法的一种扩展专门用于解决这类符号回归问题。在遗传规划中每个候选解都是一个数学表达式树树的节点可以是运算符(如加、减、乘、除)叶子节点则是变量或常数。进化过程通过选择、交叉和变异等操作不断改进这些表达式树。遗传规划的主要优势能够发现人类难以想到的复杂非线性关系结果具有可解释性不像神经网络那样的黑箱不需要特征工程自动处理变量间交互作用提示符号回归特别适合物理、工程领域的问题因为发现的数学表达式往往具有物理意义。2. gplearn库入门与实践gplearn是一个专门为符号回归设计的Python库它提供了与scikit-learn兼容的API使得使用方式非常直观。下面我们通过一个完整示例来演示如何使用gplearn进行符号回归。2.1 安装与环境准备首先安装必要的库pip install gplearn numpy matplotlib scikit-learn2.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_size5000, generations20, stopping_criteria0.01, p_crossover0.7, p_subtree_mutation0.1, p_hoist_mutation0.05, p_point_mutation0.1, max_samples0.9, verbose1, random_state0) 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, namelog, arity1) # 自定义适应度指标 def _mape(y, y_pred, w): diffs np.abs((y - y_pred) / np.maximum(0.001, y)) return 100. * np.average(diffs, weightsw) mape make_fitness(function_mape, greater_is_betterFalse) # 使用自定义函数和适应度 est_gp SymbolicRegressor(function_set[add, sub, mul, div, protected_log], metricmape)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, fitnesscreator.FitnessMin) # 定义原始集合 pset gp.PrimitiveSet(MAIN, arity2) # 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, psetpset, 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, psetpset) # 定义评估函数 def evaluate(individual, points): func toolbox.compile(exprindividual) sqerrors [(func(x[0], x[1]) - y)**2 for x, y in points] return (np.mean(sqerrors),) toolbox.register(evaluate, evaluate, pointslist(zip(X_train, y_train))) toolbox.register(select, tools.selTournament, tournsize3) toolbox.register(mate, gp.cxOnePoint) toolbox.register(expr_mut, gp.genFull, min_0, max_2) toolbox.register(mutate, gp.mutUniform, exprtoolbox.expr_mut, psetpset)3.3 进化主循环实现完整的进化过程def main(): pop toolbox.population(n300) hof tools.HallOfFame(1) stats_fit tools.Statistics(lambda ind: ind.fitness.values) stats_size tools.Statistics(len) mstats tools.MultiStatistics(fitnessstats_fit, sizestats_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, statsmstats, halloffamehof, verboseTrue) return pop, log, hof pop, log, hof main() best_ind hof[0] print(fBest individual: {best_ind}) print(fFitness: {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_size0.2, random_state42) # 配置符号回归器 est_gp SymbolicRegressor(population_size1000, generations20, function_set[add, sub, mul, div, sqrt, log, abs], metricmean absolute error, parsimony_coefficient0.01, random_state0) # 训练和评估 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(exprindividual) func(0, 0) # 测试评估 return True except: return False3. 复杂问题分解 对于高维问题可以尝试以下策略先进行特征选择或降维使用多个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-, labelMinimum Fitness) ax1.set_xlabel(Generation) ax1.set_ylabel(Fitness, colorb) ax2 ax1.twinx() ax2.plot(gen, size_avgs, r-, labelAverage Size) ax2.set_ylabel(Size, colorr) 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, progdot) nx.draw_networkx_nodes(g, pos, node_size900, node_colorw) 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, fitnesscreator.FitnessMulti) # 修改评估函数 def evaluate_mo(individual, points): func toolbox.compile(exprindividual) sqerrors [(func(x[0], x[1]) - y)**2 for x, y in points] return np.mean(sqerrors), len(individual) toolbox.register(evaluate, evaluate_mo, pointslist(zip(X_train, y_train))) toolbox.register(select, tools.selNSGA2) # 运行多目标算法 pop toolbox.population(n100) hof tools.ParetoFront() stats tools.Statistics(lambda ind: ind.fitness.values) stats.register(avg, np.mean, axis0) stats.register(std, np.std, axis0) stats.register(min, np.min, axis0) algorithms.eaMuPlusLambda(pop, toolbox, mu50, lambda_100, cxpb0.7, mutpb0.3, ngen40, statsstats, halloffamehof)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_size0.2, random_state42) # 创建符号分类器 est_clf SymbolicClassifier(population_size500, generations20, function_set[add, sub, mul, div, sqrt, log], metriclog loss, parsimony_coefficient0.01, random_state0) # 训练和评估 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(generations20, population_size200, hall_of_fame100, n_components10, function_set[add, sub, mul, div, sqrt], parsimony_coefficient0.0005, max_samples0.9, random_state0) # 构建管道 model Pipeline([ (symbolic, transformer), (regressor, Ridge()) ]) # 训练和评估 model.fit(X_train, y_train) print(Test score:, model.score(X_test, y_test))在实际项目中我发现将符号回归与其他机器学习方法结合使用往往能取得更好的效果。例如先用符号回归生成一些特征再输入到神经网络或梯度提升树中。这种方法结合了符号回归的可解释性和传统机器学习模型的强大预测能力。
