结合scipy与matplotlib来绘制曲线拟合图

在做科研论文的时候，常常需要在图中描绘某些实际数据观察的同时，使用一个曲线来拟合这些实际数据。在这里，我基于复杂网络中常用的power-law分布来介绍如何利用python进行这一类图形的绘制。

首先简单介绍一下什么是power-law。 power-law中文称作幂率分布，数学的表达式为P(x) = c*x^(-r)，其中c与r是常数。在自然界与社会生活中存在各种各样性质迥异的幂律分布现象，例如经济学中的20/80法则，表示少数人控制了社会中绝大部分的财富。更多相关的介绍可以参考http://en.wikipedia.org/wiki/Power_law。

为了在图形中直观的表述power-law分布，通常会将其转换成对数形式，即 lnP(x) = -r*lnx + lnc。 使用 x’ = lnx， y’ = lnP(x)，上式既可以轻松转换成一个线性形式： y’ = -r*x’ + lnc。下图描述了根据实际数据所绘制的双对数坐标图，可以清晰看到横坐标与纵坐标在对数转换后是一个线性函数的形式（绘图方法见文matplotlib使用小结（基本篇））。

点击查看原始尺寸

在得到上述的实际观察数据后，为了能够使用一个分布去很好的拟合这些数据，我们通常会使用最小二乘法（least squares）来进行参数估计。最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配（来自维基百科——最小二乘法词条）。利用上面的例子来说明，我们已知x’与y‘，希望能找到一个表达式 y* = -r*x’ + lnc，使得y*与y’的差别最小，即求解r与c，使得sum(y’-y*)^2最小。上述计算可以最终转换成矩阵计算，使用矩阵b表示由y’的列向量，a表示由（x’,1）组成的矩阵，t表示系数（-r,lnc）的列向量，求解方程at = b，或者是求解t使得|b-at|最小，这样我们就可以使用一些数学工具软件来轻松求解了。

接下来，我们使用到了python中的一个用于科学计算的包scipy。在scipy下有一个linalg的包专门用于进行线性代数的计算，在这里我们需要使用到的函数是lstsq，从这个函数名字就可以看出来其是专门用于最小二乘法求解的了。在这个函数中，需要使用到的参数就是前面所提得自变量a和因变量b。lstsq(a, b）返回的结果包括所求系数列向量t,残差，矩阵a的秩与奇异值。

此外另外一个用到的包是很常用的numpy,其函数mat可以将列表转换成矩阵。

整个的代码如下(获取原始数据部分的代码没有贴出来)：

import math

from scipy import linalg

import numpy as np

import matplotlib.pyplot as plt

fig = plt.figure()

ax = fig.add_subplot(111)

ax.loglog(x,y,”o”) # 绘制实际数据的图形 x,y表示原始数据的横坐标与纵坐标，没有贴出生成这些数据的代码

lnx = [math.log(f,math.e) for f in x] # 对x进行对数转换

lny = [math.log(f,math.e) for f in y] # 对y进行对数转换

a = np.mat([lnx,[1]*len(x)]).T # 自变量矩阵a

b = np.mat(lny).T # 因变量矩阵b

(t,res,rank,s) = linalg.lstsq(a,b) # 最小二乘法求系数

print t

结果为[[-2.01203744]

[ 2.27992415]]， 实际上是一个2*1的矩阵

r = t[0][0]

c = t[1][0]

x_ = x

y_ = [math.e**(r*a+c) for a in lnx] # 根据求的的系数求出y*

ax.loglog(x_,y_,”r-“) # 绘制拟合的曲线图

plt.show()

所得图形如下，可以看到拟合的很好：

点击查看原始尺寸

对于其它的分布，可以使用类似的方法进行拟合，只需要首先把分布函数通过适当的方法转换成线性形式即可。

转自：http://flyfeeling.blogbus.com/logs/61319104.html