python据说功能强大,触角伸到各个领域,网上搜了一下其科学计算和工程计算能力也相当强,具备各种第三方包,除了性能软肋外,其他无可指摘,甚至可以同matlab等专业工具一较高下。

从网上找了一个使用遗传算法实现数据拟合的例子学习了一下,确实Python相当贴合自然语言,终于编程语言也能说人话了,代码整体简洁、优雅。。

代码功能:给出一个隐藏函数 例如 z=x^2+y^2,生成200个数据,利用这200个数据,使用遗传算法猜测这些数据是什么公式生成的。 (说的太直白,一点都不高大上)

代码如下:



Python中的拟合 python怎么拟合数据_matlab

1 # coding=utf-8 
  
  2 from random import random, randint, choice,uniform 
  
  3 from copy import deepcopy 
  
  4 import numpy as np 
  
  5 import matplotlib.pyplot as plt 
  
  6  
  
  7 from random import random, randint, choice 
  
  8 from copy import deepcopy 
  
  9 import numpy as np 
  
 10  
  
 11 # 运算类 
  
 12 class fwrapper: 
  
 13     def __init__(self, function, childcount, name): 
  
 14         self.function = function 
  
 15         self.childcount = childcount 
  
 16         self.name = name 
  
 17  
  
 18 # 节点类 
  
 19 class node: 
  
 20     def __init__(self, fw, children): 
  
 21         self.function = fw.function 
  
 22         self.name = fw.name 
  
 23         self.children = children 
  
 24 #将inp指定的运算符作用到子节点上 
  
 25     def evaluate(self, inp): 
  
 26         # 循环调用子节点的子节点的子节点....的evaluate方法 
  
 27         results = [n.evaluate(inp) for n in self.children] 
  
 28         # 返回运算结果 
  
 29         return self.function(results) 
  
 30 #打印本节点及所属节点的操作运算符 
  
 31     def display(self, indent=0): 
  
 32         print(' ' * indent) + self.name 
  
 33         for c in self.children: 
  
 34             c.display(indent + 1) 
  
 35  
  
 36 #参数节点类,x+y 其中x,y都是参数节点 
  
 37 class paramnode: 
  
 38     def __init__(self, idx): 
  
 39         self.idx = idx 
  
 40 # evaluate方法返回paramnode节点值本身 
  
 41     def evaluate(self, inp): 
  
 42         return inp[self.idx] 
  
 43  
  
 44     def display(self, indent=0): 
  
 45         print '%sp%d' % (' ' * indent, self.idx) 
  
 46  
  
 47 # 常数节点 
  
 48 class constnode: 
  
 49     def __init__(self, v): 
  
 50         self.v = v 
  
 51  
  
 52     def evaluate(self, inp): 
  
 53         return self.v 
  
 54  
  
 55     def display(self, indent=0): 
  
 56         print '%s%d' % (' ' * indent, self.v) 
  
 57  
  
 58 # 操作运算符类 
  
 59 class opera: 
  
 60     # 使用前面定义的fwrapper类生产常用的加减乘除运算,第一个参数是本运算执行方式,第二个参数是本运算接受的参数个数,第三个参数是本运算名称 
  
 61     addw = fwrapper(lambda l: l[0] + l[1], 2, 'add') 
  
 62     subw = fwrapper(lambda l: l[0] - l[1], 2, 'subtract') 
  
 63     mulw = fwrapper(lambda l: l[0] * l[1], 2, 'multiply') 
  
 64  
  
 65     def iffunc(l): 
  
 66         if l[0] > 0: 
  
 67             return l[1] 
  
 68         else: 
  
 69             return l[2] 
  
 70     #定义if运算 
  
 71     ifw = fwrapper(iffunc, 3, 'if') 
  
 72  
  
 73     def isgreater(l): 
  
 74         if l[0] > l[1]: 
  
 75             return 1 
  
 76         else: 
  
 77             return 0 
  
 78     #定义greater运算 
  
 79     gtw = fwrapper(isgreater, 2, 'isgreater') 
  
 80     #构建运算符集合 
  
 81     flist = [addw, mulw, ifw, gtw, subw] 
  
 82  
  
 83     #使用node,paramnode,fwrapper构建一个example 
  
 84     def exampletree(self): 
  
 85         return node(self.ifw, [node(self.gtw, [paramnode(0), constnode(3)]), node(self.addw, [paramnode(1), constnode(5)]), 
  
 86                           node(self.subw, [paramnode(1), constnode(2)]), ]) 
  
 87  
  
 88  
  
 89     # 构建一颗随机运算数,pc为参数(分叉)个数,maxdepth为树的深度,fpr为运算符个数在运算符加节点总数中所占比例,ppr为参数个数在参数加常数个数总数中所占的比例 
  
 90     def makerandomtree(self,pc, maxdepth=4, fpr=0.5, ppr=0.6): 
  
 91         if random() < fpr and maxdepth > 0: 
  
 92             f = choice(self.flist) 
  
 93             # 递归调用makerandomtree实现子节点的创建 
  
 94             children = [self.makerandomtree(pc, maxdepth - 1, fpr, ppr) for i in range(f.childcount)] 
  
 95             return node(f, children) 
  
 96         elif random() < ppr: 
  
 97             return paramnode(randint(0, pc - 1)) 
  
 98         else: 
  
 99             return constnode(randint(0, 10)) 
  
100  
  
101  
  
102     #变异,变异概率probchange=0.1 
  
103     def mutate(self,t, pc, probchange=0.1): 
  
104         # 变异后返回一颗随机树 
  
105         if random() < probchange: 
  
106             return self.makerandomtree(pc) 
  
107         else: 
  
108             result = deepcopy(t) 
  
109             # 递归调用,给其子节点变异的机会 
  
110             if isinstance(t, node): 
  
111                 result.children = [self.mutate(c, pc, probchange) for c in t.children] 
  
112             return result 
  
113  
  
114     #交叉 
  
115     def crossover(self,t1, t2, probswap=0.7, top=1): 
  
116         # 如果符合交叉概率,就将t2的值返回,实现交叉; 
  
117         if random() < probswap and not top: 
  
118             return deepcopy(t2) 
  
119         else: 
  
120             #如果本层节点未实现交配,递归询问子节点是否符合交配条件 
  
121             #首先使用deepcopy保存本节点 
  
122             result = deepcopy(t1) 
  
123             if isinstance(t1, node) and isinstance(t2, node): 
  
124                 #依次递归询问t1下的各子孙节点交配情况,交配对象为t2的各子孙;t1,t2家族同辈交配 
  
125                 result.children = [self.crossover(c, choice(t2.children), probswap, 0) for c in t1.children] 
  
126             return result 
  
127  
  
128     # random2.display() 
  
129     # muttree=mutate(random2,2) 
  
130     # muttree.display() 
  
131     # cross=crossover(random1,random2) 
  
132     # cross.display() 
  
133  
  
134     #设置一个隐藏函数,也就是原始函数 
  
135     def hiddenfunction(self,x, y): 
  
136         return x ** 2+ y**2 
  
137  
  
138     #依照隐藏函数,生成坐标数据 
  
139     def buildhiddenset(self): 
  
140             rows = [] 
  
141             for i in range(200): 
  
142                 x = randint(0, 10) 
  
143                 x=uniform(-1,1) 
  
144                 y = randint(0, 10) 
  
145                 y=uniform(-1,1) 
  
146                 rows.append([x, y, self.hiddenfunction(x, y)]) 
  
147             print("rows:",rows) 
  
148             return rows 
  
149  
  
150     #拟合成绩函数,判定拟合函数(实际是一颗图灵树)贴近原始函数的程度 
  
151     def scorefunction(self,tree, s): 
  
152         dif = 0 
  
153         # print("tree:",tree) 
  
154         # print("s:",s) 
  
155         for data in s: 
  
156             # print("data[0]:",data[0]) 
  
157             # print("data[1]:",data[1]) 
  
158             v = tree.evaluate([data[0],data[1]]) 
  
159             #累加每个数据的绝对值偏差 
  
160             dif += abs(v - data[2]) 
  
161         return dif 
  
162  
  
163     #返回一个成绩判定函数rankfunction的句柄 
  
164     def getrankfunction(self,dataset): 
  
165         #此函数调用拟合成绩函数,并对成绩排序,返回各个种群的成绩 
  
166         def rankfunction(population): 
  
167             scores = [(self.scorefunction(t, dataset), t) for t in population] 
  
168             scores.sort() 
  
169             return scores 
  
170         return rankfunction 
  
171  
  
172     # hiddenset=buildhiddenset() 
  
173     # scorefunction(random2,hiddenset) 
  
174     # scorefunction(random1,hiddenset) 
  
175  
  
176     def evolve(self,pc, popsize, rankfunction, maxgen=500, mutationrate=0.1, breedingrate=0.4, pexp=0.7, pnew=0.05): 
  
177         #轮盘算法 
  
178         def selectindex(): 
  
179             return int(np.log(random()) / np.log(pexp)) 
  
180         #使用随机树生成第一代各种群 
  
181         population = [self.makerandomtree(pc) for i in range(popsize)] 
  
182         #计算每一代各种群的成绩, 
  
183         for i in range(maxgen): 
  
184             scores = rankfunction(population) 
  
185             #打印历代最好成绩 
  
186             print('the best score in genneration ',i,':',scores[0][0]) 
  
187             #成绩完全吻合原函数的话,退出函数 
  
188             if scores[0][0] == 0: 
  
189                 break 
  
190             #创建新一代各种群,成绩前两名的直接进入下一代 
  
191             newpop = [scores[0][1], scores[1][1]] 
  
192             while len(newpop) < popsize: 
  
193                 #pnew为引进随机种群概率,未达此概率的,使用原种群的交配、变异生成新种群 
  
194                 if random() > pnew: 
  
195                     newpop.append( 
  
196                         self.mutate(self.crossover(scores[selectindex()][1], scores[selectindex()][1], probswap=breedingrate), pc, 
  
197                                probchange=mutationrate)) 
  
198                 #引入随机种群 
  
199                 else: 
  
200                     newpop.append(self.makerandomtree(pc)) 
  
201             population = newpop 
  
202             #打印历代最好种群 
  
203             # scores[0][1].display() 
  
204         return scores[0][1] 
  
205  
  
206  
  
207  
  
208 def main(argv): 
  
209     e=opera() 
  
210     def exampletree(): 
  
211         return node(e.ifw,[node(e.gtw,[paramnode(0),constnode(3)]),node(e.addw,[paramnode(1),constnode(5)]),node(e.subw,[paramnode(1),constnode(2)])]) 
  
212  
  
213  
  
214     # random1=e.makerandomtree(2) 
  
215     # random1.evaluate([7,1]) 
  
216     # random1.evaluate([2,4]) 
  
217     # random2=e.makerandomtree(2) 
  
218     # random2.evaluate([5,3]) 
  
219     # random2.evaluate([5,20]) 
  
220     # random1.display() 
  
221     # random2.display() 
  
222  
  
223     # exampletree = e.exampletree() 
  
224     # exampletree.display() 
  
225     # print(exampletree.evaluate([6, 1])) 
  
226     # print('after evaluate:') 
  
227     # exampletree.display() 
  
228     # exampletree.evaluate([2, 3]) 
  
229     # exampletree.evaluate([5, 3]) 
  
230     # exampletree.display() 
  
231  
  
232     a=opera() 
  
233     row2=a.buildhiddenset() 
  
234     # fig=plt.figure() 
  
235     # ax=fig.add_subplot(1,1,1) 
  
236     # plt.plot(np.random.randn(1000).cumsum()) 
  
237     # plt.show() 
  
238  
  
239  
  
240  
  
241     from mpl_toolkits.mplot3d import Axes3D 
  
242     fig = plt.figure() 
  
243     ax = fig.add_subplot(111, projection='3d') 
  
244     X = [1, 1, 2, 2] 
  
245     Y = [3, 4, 4, 3] 
  
246     Z = [1, 2, 1, 1] 
  
247     rx=[] 
  
248     ry=[] 
  
249     rz=[] 
  
250     for i in row2: 
  
251         rx.append(i[0]) 
  
252         ry.append(i[1]) 
  
253         rz.append(i[2]) 
  
254  
  
255     ax.plot_trisurf(rx, ry, rz) 
  
256     rz2=[] 
  
257     rf = a.getrankfunction(row2) 
  
258     final = a.evolve(2, 100, rf, mutationrate=0.2, breedingrate=0.1, pexp=0.7, pnew=0.1,maxgen=500) 
  
259     print('__________________is it?_________________________') 
  
260     final.display() 
  
261     for j in range(0,len(rx)): 
  
262         rz2.append(final.evaluate([rx[j],ry[j]])) 
  
263     fig2 = plt.figure() 
  
264     ax2 = fig2.add_subplot(111, projection='3d') 
  
265     ax2.plot_trisurf(rx, ry, rz2) 
  
266  
  
267     plt.show() 
  
268  
  
269  
  
270     # print(rf) 
  
271     # final = a.evolve(2, 500, rf, mutationrate=0.2, breedingrate=0.1, pexp=0.7, pnew=0.1) 
  
272     # print("final:",final) 
  
273     # print(final.evaluate([1,8])) 
  
274     # print(final.evaluate([2,9])) 
  
275  
  
276  
  
277  
  
278  
  
279  
  
280 if __name__=="__main__": 
  
 
 
281     main(0)

 

看懂不一定写的出来,这是这次写这个程序最大的体会, 得定时拿出来复习复习。