主要利用Scipy进行一些曲线的拟合,基本会用到对线性曲线的拟合和非线性曲线的拟合。
1.利用最小二乘进行拟合
在Scipy的optimize子包中,可以利用leastsq进行最小二乘拟合。
方法1:使用矩阵运算
(这里,假设节点和边已经添加到G中,而且这里的G为有向图)
'''
sort the out degree of the nodes in G
store the unsorted out_degree:nodes in dict_out_degree
store the frequency of out_degree in sorted_out_degree
'''
node_num=G.number_of_nodes()
edge_num=G.number_of_edges()
out_degree=G.out_degree()
dict_out_degree={}
sorted_out_degree=[]
for i in G.nodes():
if((out_degree[i] in dict_out_degree)==True):
dict_out_degree[out_degree[i]]=dict_out_degree[out_degree[i]]+1
else:
dict_out_degree[out_degree[i]]=1
index=0
while index<node_num:
if((index in dict_out_degree)==True):
sorted_out_degree.append(dict_out_degree[index])
else:
sorted_out_degree.append(0)
index=index+1
'''
sort the in_degree of the nodes in G
store the unsorted in_degree:nodes in dict_in_degree
store the frequency of in_degree in sorted_in_degree
'''
in_degree=G.in_degree()
dict_in_degree={}
sorted_in_degree=[]
for i in G.nodes():
if((in_degree[i] in dict_in_degree)==True):
dict_in_degree[in_degree[i]]=dict_in_degree[in_degree[i]]+1
else:
dict_in_degree[in_degree[i]]=1
index=0
while index<node_num:
if((index in dict_in_degree)==True):
sorted_in_degree.append(dict_in_degree[index])
else:
sorted_in_degree.append(0)
index=index+1
xin=range(len(sorted_in_degree))
degree_in=sorted_in_degree
yin=[z/float(sum(degree_in)) for z in degree_in]
'''
store out_degree:frequecy in x:y
'''
x=range(len(sorted_out_degree))
degree=sorted_out_degree
y=[z/float(sum(degree)) for z in degree]
x0=[]
y0=[]
'''
必须把非零元素去掉,因为之后要做log运算
''' for z in x: if z!=0 and y[z]!=0: x0.append(z) y0.append(y[z]) '''这个范围是本人自己定的,可以根据实际的曲线,截取有效的点范围''' xm=x0[11:31] ym=y0[11:31] lnx=[math.log(f,math.e) for f in xm] lny=[math.log(f,math.e) for f in ym] a=np.mat([lnx,[1]*len(xm)]).T b=np.mat(lny).T (t,res,rank,s)=linalg.lstsq(a,b) r=t[0][0]#这里的r即为直线斜率,而c则为直线中的b c=t[1][0] x_=xm y_=[math.e**(r*a+c) for a in lnx] plt.loglog(x,y,'go',x_,y_,'r-',xin,yin,'b+')
plt.xlabel('The degree of the vertices ( n='+str(node_num)+', m='+str(edge_num)+')')plt.legend(['out_degree','$r$','in_degree'])
plt.ylabel('The probability of the degree')
plt.text(8,0.1,'$r$'+'='+str(math.fabs(r)),color="blue",fontsize=12)
plt.savefig('D:/'+str(node_num)+'.png')
方法二:中规中矩的教科书方案。和方法一的区别在哪里,还不太清楚plt.show()
2.利用curve_fit进行非线性最小二乘拟合import numpy as np from scipy.optimize import leastsq #待拟合的函数,x是变量,p是参数 def fun(x, p): a, b = p return a*x + b #计算真实数据和拟合数据之间的误差,p是待拟合的参数,x和y分别是对应的真实数据 def residuals(p, x, y): return fun(x, p) - y #一组真实数据,在a=2, b=1的情况下得出。这里只需要将x1和y1换成上面中的Inx和Iny即可 x1 = np.array([1, 2, 3, 4, 5, 6], dtype=float) y1 = np.array([3, 5, 7, 9, 11, 13], dtype=float) #调用拟合函数,第一个参数是需要拟合的差值函数,第二个是拟合初始值,第三个是传入函数的其他参数 r = leastsq(residuals, [1, 1], args=(x1, y1)) #打印结果,r[0]存储的是拟合的结果,r[1]、r[2]代表其他信息 print r[0]
当拟合的函数中有exp和log出现的时候,发现采用leastsq拟合会报错,原来是exp和log不支持数组运算。不过还好,可以使用curve_fit
from scipy.optimize import curve_fit '''在func中定义了拟合的函数,这个自己设置的''' def func(x,a,b): return (b-a/np.log(x)) def power_law(): x=np.array([1098,1715,3430,8913,9738,14012,15632,22920,25852,31693,41013],dtype=float) y=np.array([4.968,4.5149,4.284,3.612,3.3046,3.3751,3.219,3.02,2.9053,2.76331,2.71217],dtype=float) y0=func(x,16.667,2) popt,pcov=curve_fit(func,x,y) x=np.arange(xnodes[0],xnodes[10],1000) a,b=float(popt[0]),float(popt[1]) tt=plt.plot(xnodes,ylaw,'g+',x,func(x,a,b),color='blue') plt.legend(tt,['origin_data','fitting_curve'],numpoints=1) plt.show() return a,b
相关参考:
http://flyfeeling.blogbus.com/logs/61319104.html
http://marinzou.blogbus.com/logs/65593603.html