networkx使用笔记(三)之好汉篇Scipy(3)

主要利用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.legend(['out_degree','$r$','in_degree'])

plt.xlabel('The degree of the vertices ( n='+str(node_num)+', m='+str(edge_num)+')')

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()

方法二：中规中矩的教科书方案。和方法一的区别在哪里，还不太清楚

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]

2.利用curve_fit进行非线性最小二乘拟合

当拟合的函数中有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      http://z-wang.com/blog/python-least-square/          http://blog.sina.com.cn/s/blog_5f234d4701013ln6.html