基因表达式编程(Gene Expression Programming,GEP)
前言
该算法旨在在一组数据点中,用基因表达式编程的方法,根据基因遗传定律,物竞天择、优者生存,劣者淘汰的思想,不断进化种群,找出适宜度最高的染色体来模拟出数据点之间所存在的数学表达式关系。通常该算法用来解决符号回归问题:符号回归(Symbolic Regression)作为一种一种监督学习方法,试图发现某种隐藏的数学公式,以此利用特征变量预测目标变量。符号回归的优点就是可以不用依赖先验的知识或者模型来为非线性系统建立符号模型。符号回归基于进化算法,它的主要目标就是利用进化方法综合出尽可能好的解决用户自定义问题的方法(数学公式,计算机程序,逻辑表达式等)。
1 编码规则
1.1 基因
一段基因由两部分组成:头部和尾部,形如:
其中,红色部分为头部基因,其长度为9,绿色部分为尾部基因,其长度为8。
基因中的每个元素分为两种形式:函数和终点,它们分别从函数集和终点集中进行选择:
- 函数集
- 算术运算符:例如+、-、*、/等
- 初等数学函数,例如:sin、开方、cos等
- 其它一些函数,例如:max、min等
- 布尔运算符:例如与、或、非
- 关系运算符:例如<、>、=等
- 条件运算符:例如 if 等
- 终点集
- 变量数值:例如\(x_1\) 、\(x_2\)
- 常量数值:例如 \(\pi\)、e
- 无参函数:例如rand()等
基因头部长度我们用h表示,尾部长度用t表示,函数集所需最大参数用n表示,比如一个函数集为{+、-、*、/、sin、cos},其中所需参数的最大个数为2,则n为2
那么h、t、n三者之间的关系是:\(t=h*(n-1)+1\)
维持这样的关系,可以保证产生的基因编码都是合法编码。
1.2 基因编码原理
比如下方的数学表达式:
翻译成表达式树为,其中Q
表示开方:
那么可以得到它的基因型为:
下面来看看一个基因型翻译成表达式的过程:
基因型如下:
基因的起点对应表达式树的根:
Q
(开方)需要一个参数,那么往后取一个参数:
乘法需要两个参数,往后取两个参数:
以此类推,可以画出该基因型的表达式树为:
我们可以看出,表达式树的根一定要是函数集的元素,叶子节点一定是终点集的元素。
1.3 染色体
一条染色体可以由多个基因构成,那么我们需要设定基因连接符,比如用+、-、 *、/等符号进行连接,比如:
这三个子树我们可以通过+
进行连接:
其次一条染色体代表一个个体,那么一个种群就可以有多个染色体。
1.4 Dc域
在符号回归问题种,我们使用这样的结构来操纵随机数值常量。例如,如下的染色体:
含有一个附加域 Dc,这个附加域用来对随机数值常量进行编码。该染色体的翻译过程如下图所示:
1.5 进化过程
基因的进化分为以下几种:
- 变异
- 转座
- IS转座
- RIS转座
- 基因转座
- 重组
- 单点重组
- 两点重组
- 基因重组
下面我们对进化操作分别进行讲解:
1.5.1 变异
变异的原则:
- 基因首位元素必须是函数集元素
- 除过首位外的头部基因可以从函数集或者终点集进行选择
- 尾部元素必须是终点集元素
变异过程将随机从种群中选择个体(染色体),再从染色体中随机选择一个基因,从基因中随机选择一个位置。
void Gene::mutation()
{
int index = rand() % gene_len;
char ch;
if (index <= HEAD_LEN)
ch = getRandomElem();
else
ch = getTermElem();
text[index] = ch;
}
1.5.2 转座
1.5.2.1 IS转座
IS转座是从种群中随机选择一个染色体,然后随机从IS转座长度中选择IS长度,然后在染色体中选择IS长度的基因片段,并随机选择基因,插入到除基因首元素之外的头部部分中,如下所示:
void Individual::transposition_is()
{
// 随机选取转座长度
int len_array = sizeof(IS_LEN) / sizeof(int);
int index = rand() % len_array;
int len_is = IS_LEN[index];
// 从整段染色体中随机选取len_is个长度的片段
index = rand() % (len - len_is); // 随机起始索引
string str = content_without_space().substr(index, len_is); // 截取字符串
// 随机选择基因
index = index_rand(); gene[index].transposition(str);
}
1.5.2.2 RIS转座
所有的RIS元素都是从一个函数开始,因此是选自头部分中的序列。因此在头中任选一点,沿基因向后查找,直到发现一个函数为止。该函数成为RIS元素的起始位置。如果找不到任何函数,则变换不执行任何操作。该算子随机选取染色体,需要修饰的基因,RIS元素以及其长度。变化过程如图所示:
void Individual::transposition_ris()
{
// 随机选择基因
int index = index_rand();
// 在头中随机选取一点
int pos = rand() % HEAD_LEN;
// 沿该点向后查找直到发现一个函数 该点即为RIS转座的起始点
int a;
if ((a = gene[index].find_func(pos)) == -1)
return;
pos = a + Gene::length() * index;
// 随机选取RIS转座长度
int len_array = sizeof(RIS_LEN) / sizeof(int);
index = rand() % len_array; int len_ris = RIS_LEN[index];
// 从起始点向后取len_ris个长度的片段 string str = content_without_space().substr(pos, len_ris);
index = index_rand(); gene[index].transposition(str);
}
1.5.2.3 基因转座
基因转座仅仅改变基因在同一染色体中的位置,如图所示:
void Individual::transposition_gene()
{
// 随机选取两个不同位置的基因
if (GENE_NUM > 1)
{
int index1 = 0;
int index2 = 0;
while (index1 == index2)
{
index1 = index_rand();
index2 = index_rand();
}
Gene temp(gene[index1]);
gene[index1] = gene[index2];
gene[index2] = temp;
}
}
1.5.3 重组
1.5.3.1 单点重组
进行单点重组的时候,父代染色体相互配对并在相同的位置切断,两个染色体相互交换重组点之后的部分,如图所示:
void Population::recombination_one()
{
double prob = rand() / double(RAND_MAX);
if (prob < PROB_RECOMBE_ONE)
{
// 随机选取两个不同个体
int index1 = 0;
int index2 = 0;
while (index1 == index2)
{
index1 = rand() % count;
index2 = rand() % count;
}
// 随机选取重组点
int pos = rand() % Individual::length();
string subStr1 = individual[index1].content_without_space().substr(0, pos);
string subStr2 = individual[index2].content_without_space().substr(0, pos);
// 交换两个染色体的片段 individual[index1].recombination(0, pos, subStr2);
individual[index2].recombination(0, pos, subStr1);
}
}
1.5.3.2 两点重组
进行两点重组的时候,父代染色体相互配对,在染色体中随机选择两个点,将染色体切断。两个染色体相互交换重组点之间的部分,形成两个新的子代染色体。如图所示:
void Population::recombination_two()
{
double prob = rand() / double(RAND_MAX);
if (prob < PROB_RECOMBE_TWO)
{
// 随机选取两个不同个体
int index1 = 0;
int index2 = 0;
while (index1 == index2)
{
index1 = rand() % count;
index2 = rand() % count;
}
// 随机选取左、右重组点
int pos_left = rand() % Individual::length();
int pos_right = (rand() % (Individual::length() - pos_left)) + pos_left;
int len = pos_right - pos_left + 1; // 重组点之间的距离
// 重组点之间的字符串,包括左、右重组点
string subStr1 = individual[index1].content_without_space().substr(pos_left, len);
string subStr2 = individual[index2].content_without_space().substr(pos_left, len);
// 交换两个染色体的片段
individual[index1].recombination(pos_left, len, subStr2);
individual[index2].recombination(pos_left, len, subStr1);
}
}
1.5.3.3 基因重组
在 GEP 的第三种重组中,两个染色体中的整个基因相互交换,形成的两个子代染色体含有来自两个父体的基因。
void Population::recombination_gene()
{
double prob = rand() / double(RAND_MAX);
if (prob < PROB_RECOMBE_GENE)
{
// 随机选取两个不同个体 int index1 = 0;
int index2 = 0; while (index1 == index2)
{
index1 = rand() % count; index2 = rand() % count;
}
// 随机选取基因索引
int index = rand() % GENE_NUM;
// 获取两个个体相应索引的基因
string subStr1 = individual[index1].content_without_space().substr(index * Gene::length(), Gene::length());
string subStr2 = individual[index2].content_without_space().substr(index * Gene::length(), Gene::length());
// 交换
individual[index1].recombination(index * Gene::length(), Gene::length(), subStr2);
individual[index2].recombination(index * Gene::length(), Gene::length(), subStr1);
}
}
1.6 适应度函数
个体程序 的适应度 可以用第一个方程来表示,当误差类型选为相对误差的时候,个体程序的适应度 可以用方程(3.1b)来表示:
其中\(M\)是选择范围,\(C_{(i,j)}\)是染色体个体\(i\)对于适应度样本\(j\)(来自\(C_t\)集合中)的返回值。而\(T_j\)是适应度样本\(j\)的目标值。请注意,两个方程中的绝对值项都对应各自的误差(前一个方程中对应于绝对误差,后一个方程中对应于相对误差)。该项称为精度\(p\)。因此,对于一个完美适应的情况, 且 \(C_{(i, j)}=C_t\)且\(f_i=f_{max}=C_tM\)。
1.7 个体的选择
在 GEP 中,个体通过赌盘轮采样策略根据其适应度进行选择。每个个体用圆形赌盘的一块来代表其适应度的比例。赌盘按照群体中个体数的值进行相应次数的旋转,从而始终保持群体的大小不变。
采用这种选择策略,确实有时候会丢失一些最佳个体,而一些普通个体进入到了下一代。但是这样也不一定就不好,因为种群会一代一代地向前推进。而且,因为我们使用每一代中的最佳个体复制策略,所以最佳个体的存在和复制能够得到保证。这样,至少最优的特性从来没有丢失过,并且能够达到不断学习的目的。
还有几种其他最常用的选择算法——赌盘选择、确定性选择、锦标赛选择。
2 符号回归
2.1 算法流程图
2.2 测试函数及数据
2.3 参数设置
参数名称 | 参数值 |
---|---|
函数集 | |
终点集 | a |
头部长度 | 7 |
是否开启Dc域 | True |
Dc域最大值 | 10.0 |
Dc域最小值 | 0.0 |
Dc域元素个数 | 10 |
基因数目 | 4 |
基因连接符 | + |
种群大小 | 20 |
繁衍代数 | 1000 |
变异概率 | 0.043 |
IS转座概率 | 0.1 |
RIS转座概率 | 0.1 |
IS元素长度 | 3 |
RIS元素长度 | 3 |
基因转座概率 | 0.1 |
单点重组概率 | 0.8 |
两点重组概率 | 0.8 |
基因重组概率 | 0.8 |
是否选用绝对误差 | True |
选择范围 | 1000.0 |
2.4 测试结果
经过化简,得到其数学表达式为: 其图像如下所示,其中橙色曲线为函数轨迹,红色点为数据点3 部分代码解读
3.1 基因型转化为表达式树
也称为:K-表达树解码算法
3.1.1 算法流程图
3.1.2 伪代码
- 输入:基因型表达式串K
- 输出:表达式树T
设置队列Q为空;
p->0; // 定义K的下标指针
T->new Node(Kp);
p->p+1;
将T添加到队列Q的尾部;
while(Q != NULL)
{
从Q头部取出一个元素n;
if(n是终结符)
continue;
a->f(n); // 取得函数n的参数个数
for(i->1 to a)
{
构造节点m->new Node(Kp);
将节点m作为节点n的第i个子节点;
p->p+1; 将m添加到Q尾部;
}
}
return T;
3.1.3 代码实现
std::string Gene::decode()
{
queue<char> queue = {};
int index = 0; // 基因片段字符串的下标索引
tree = new BinaryTree(text.at(index));
queue.push(text.at(index));
index++;
// 构造二叉树
while (!queue.empty())
{
char ch = queue.front();
queue.pop();
if (isTerm(ch))
continue;
else
{
int count = paramCount(ch);
BinaryTreeNode* parent = tree->Find(ch);
for (int i = 0; i < count; i++)
{
char value = text.at(index);
if (count == 1)
tree->Insert(parent, value, false);
else
tree->Insert(parent, value);
index++;
queue.push(value);
}
}
}
// 获取模型复杂度
complexity = tree->size();
// 中序遍历二叉树
string result = tree->Output();
delete tree;
return result;
}
3.2 中缀表达式转后缀表达式
规则:从左到右遍历中缀表达式的每个数字和符号,若是数字就输出,即成为后缀表达式的一部分;若是符号,则判断其与栈顶元素的优先级,是右括号或优先级低于栈顶符号则栈顶元素依次出栈并输出,并将当前符号进栈,一直到最终输出后缀表达式为止。
std::queue<char> Gene::infix2postfix(string expression)
{
// 先将字符串表达式依次入队
queue<char> temp;
for (string::iterator it = expression.begin(); it != expression.end(); ++it)
{
temp.push(*it);
}
queue<char> postfix;
stack<char> charStack;
while (!temp.empty())
{
// 弹出字符串队列的队头元素
char ch = temp.front(); temp.pop();
if (isTerm(ch))
{
// 如果是运算数 直接入队
postfix.push(ch);
}
else if (ch == '(')
{
// 如果是左括号 压入堆栈
charStack.push(ch);
}
else if (ch == ')')
{
// 如果是右括号 弹出栈顶运算符并输出 直到遇到左括号(出栈 不输出)
char elem;
while (!charStack.empty() && (elem = charStack.top()) != '(')
{
postfix.push(elem);
charStack.pop();
}
// 弹出左括号
if(!charStack.empty())
charStack.pop();
}
else if (isFunc(ch))
{
// 如果是运算符
// 如果堆栈为空,或栈顶运算符为左括号“(”,则直接将此运算符入栈
// 当优先级大于栈顶运算符时,则把它压栈;
// 当优先级小于或等于栈顶运算符时,将栈顶运算符弹出并输出;
// 再比较新的栈顶运算符,直到该运算符大于栈顶运算符优先级为止,然后将该运算符 压栈
if (charStack.empty() || charStack.top() == '(')
{
charStack.push(ch);
}
else
{
int curPriority = priority(ch);
int topPriority = -1;
while (!charStack.empty() && curPriority <= (topPriority = priority(charStack.top())))
{
char c = charStack.top();
charStack.pop(); postfix.push(c);
}
charStack.push(ch);
}
}
}
// 最后弹出栈内剩余元素
while (!charStack.empty())
{
postfix.push(charStack.top());
charStack.pop();
}
return postfix;
}
3.3 计算后缀表达式
规则:从左到右遍历表达式的每个数字和符号,遇到是数字就进栈,遇到是符号,就将处于栈顶两个数字出栈,进行运算,运算结果进栈,一直到最终获得结果。
double Gene::calculate(queue<char> postfix, map<char, double> value)
{
stack<double> temp; double result = 0.0;
while (!postfix.empty())
{
char ch = postfix.front();
if (isTerm(ch))
{
// 若是数字则进栈
temp.push(value.at(ch));
}
else if (isFunc(ch))
{
// 如果是运算符 则弹出n个其需要的数字
int num = paramCount(ch);
if (num == 1)
{
double value = temp.top();
temp.pop();
result = mathExpression(ch, value);
}
else if (num == 2)
{
double value1 = temp.top();
temp.pop(); double value2 = temp.top();
temp.pop(); result = mathExpression(value1, ch, value2);
}
// 运算结果进栈
temp.push(result);
}
postfix.pop();
}
return result;
}
3.4 轮盘赌算法
3.4.1 基本思想
个体被选中的概率与其适应度函数值成正比。设群体大小为n,个体i的适应度为Fi,则个体i被选中遗传到下一代群体的概率为:
3.4.2 工作过程
设想群体全部个体的适当性分数由一张饼图来代表,如下图所示:
群体中每一染色体指定饼图中一个小块。块的大小与染色体的适应性分数成比例,适应性分数愈高,它在饼图中对应的小块所占面积也愈大。为了选取一个染色体,要做的就是旋转这个轮子,直到轮盘停止时,看指针停止在哪一块上,就选中与它对应的那个染色体。
3.4.3 代码实现
void Population::evolution()
{
// 用轮盘赌算法选择优秀个体
// 计算每个个体在轮盘中占到的概率大小
double allFitness = 0.0; // 总的适应度
vector<double> probFitness;
// 计算适应度的和
for (int i = 0; i < INDIV_NUM; i++)
{
allFitness += individual[i].getFitnessError();
}
// 计算每个适应度占总适应度的概率
for (int i = 0; i < INDIV_NUM; i++)
{
probFitness.push_back(individual[i].getFitnessError() / allFitness);
}
// 左值为左区间 右值为右区间
vector<pair<double,double> > section;
// 将概率展开为长度为1的线段 划分每个个体所占区间
double temp = 0.0;
for (int i = 0; i < INDIV_NUM - 1; i++)
{
section.push_back(pair<double, double>(temp, temp + probFitness[i]));
if(i < INDIV_NUM -2)
temp += probFitness[i];
}
section.push_back(pair<double, double>(temp, 1.0));
vector<string> context = {};
for (int i = 0; i < INDIV_NUM; i++)
{
// 生成0~1间的随机算子
double prob = rand() / double(RAND_MAX);
for (int j = 0; j < INDIV_NUM; j++)
{
// 如果落在该区间内 则加入该个体
if (prob >= section[j].first && prob <= section[j].second)
{
context.push_back(individual[j].content_without_space());
break;
}
}
}
}
4 参考代码
见github地址:Kohirus/GeneExpressionProgramming
Github 上的代码中加入了多目标优化算法(NSGA-II),代码写的有点烂,写完一直没改过,将就点看😂
5 参考文献
- 《基因表达式编程 —— 用一种人工智能方法进行数学建模》