赵工的个人空间

Scipy是一个开源的Python算法库和数学工具包，建立在NumPy之上，扩展了许多科学计算的库函数，如插值、信号处理、优化、线性代数、统计等。参见www.scipy.org网站。
优化函数程序包：提供了解决单变量和多变量的目标函数最小值问题的功能，可以通过算法和方法解决最小化问题。一般情况下，会使用线性回归搜索函数的最大值与最小值。
数值分析程序包：实现了大量的数值分析算法和方法，支持单变量和多变量的插值，以及一维和多维的样条插值，还支持拉格朗日和泰勒多项式插值方法。
统计学模块：支持概率分布函数包括连续变量分布函数、多变量分布函数以及离散变量分布函数。统计函数包括简单的均值到复杂的统计学概念，包括偏度skewness、峰度kurtosis、卡方检验chi-square test等。
聚类程序包与空间算法：聚类是一项流行的数据挖掘技术，生物学、粒子物理学、天文学、生命科学与生物信息学都广泛使用聚类分析解决问题，聚类也用于安全分析与图像处理等，包括K-means聚类、向量化、层次聚类与凝聚式聚类函数等。
图像处理函数：包括基本的图像文件读/写、图像显示，以及简单的处理函数，如裁剪、翻转和旋转，还支持图像过滤函数，包括图像数学变换、平滑、降噪、锐化，还支持图像分割、分类及边缘检测特征提取。
SciPy包含的子模块有：

1. 数学物理常数constants：

sbipy.constants模块包含大量用于科学计算的常数：

2. 特殊函数special：

scipy.special模块包含大量函数，包括基本数学函数和很多特殊函数。

3. 多维图像处理scipy.ndimage：

图像处理和分析可以看作对二维数组的操作。这个包提供了图像处理和分析的各种函数。
为图像加入噪声，然后使用滤波器去除噪声，使用高斯滤波器、中值滤波器和维纳滤波器。
1）图像滤波：
import numpy as np
from scipy import signal
from scipy import misc
from scipy import ndimage
import matplotlib.pyplot as plt
lena=misc.face(gray=True)
noisy_lena=np.copy(lena).astype(np.float)
noisy_lena+=lena.std()*0.5*np.random.standard_normal(lena.shape)
f,axarr=plt.subplots(2,2)
axarr[0,0].imshow(noisy_lena,cmap=plt.cm.gray)
axarr[0,0].axis('off')
axarr[0,0].set_title('Noissy Lena Image')
blurred_lena=ndimage.gaussian_filter(noisy_lena,sigma=3)
axarr[0,1].imshow(blurred_lena,cmap=plt.cm.gray)
axarr[0,1].axis('off')
axarr[0,1].set_title('Blurred Lena')
median_lena=ndimage.median_filter(blurred_lena,size=5)
axarr[1,0].imshow(median_lena,cmap=plt.cm.gray)
axarr[1,0].axis('off')
axarr[1,0].set_title('Median Filter Lena')
wiener_lena=signal.wiener(blurred_lena,(5,5))
axarr[1,1].imshow(wiener_lena,cmap=plt.cm.gray)
axarr[1,1].axis('off')
axarr[1,1].set_title('Wiener Filter Lena')
plt.show()
2）数学形态学：
import numpy as np
from scipy import misc,ndimage
import matplotlib.pyplot as plt
square=np.zeros((32,32)) # 全零数组
square[10:20,10:20]=1 # 其中一部分设为1
x,y=(32*np.random.random((2,15))).astype(np.int) # 随机位置
square[x,y]=1 # 把随机位置设为1
plt.imshow(square)
plt.show()
open_square=ndimage.binary_opening(square) # 开运算
plt.imshow(open_square)
plt.show()
eroded_square=ndimage.binary_erosion(square) # 膨胀运算
plt.imshow(eroded_square)
plt.show()
closed_square=ndimage.binary_closing(square) # 闭运算
plt.imshow(closed_square)
plt.show()
3）图像测量：
import numpy as np
from scipy import misc,ndimage
import matplotlib.pyplot as plt
image=misc.face()
ndimage.measurements.maximum(image)
ndimage.measurements.maximum_position(image)
ndimage.measurements.mean(image)
ndimage.measurements.median(image)
ndimage.measurements.sum(image)
ndimage.measurements.variance(image)
ndimage.measurements.standard_deviation(image)
ndimage.measurements.histogram(image,0,255,256)

4. 积分scipy.integrate：

scipy.integrate子包中有不同的积分函数。当函数对象给定时，可以用几种积分函数：
·quad：通用积分
·dblquad：通用二重积分
·tplquad：通用三重积分
·nquad：通用N重积分
·fixed_quad：对func(x)做N维高斯积分
·quadrature：在给定容限范围内的高斯积分
·romberg：对函数做Romberg积分
固定样本给定时的积分函数：
·cumtrapz：用梯形积分法计算积分
·simps：用辛氏法则从样本中计算积分
·romb：用Romberg积分法从(2**k+1)个均匀间隔的样本中计算积分
常微分方程中的积分函数：
·odeint：差分方程的通用积分
·ode：用VODE和ZVODE的方式进行ODE积分
·complex_ode：将复数值的ODE转化成实数并积分
1）quad函数执行函数的通用积分，积分范围在负无穷大到正无穷大之间。
import numpy as np
from scipy import special
from scipy import integrate
result=integrate.quad(lambda x:special.jv(4,x), 0, 20)
print(result)
print("Gaussian integrel", np.sqrt(np.pi), integrate.quad(lambda x:np.exp(-x**2), -np.inf, np.inf))
2）如果函数在积分时需要额外的参数，例如变量进行乘或乘方的因子，那么这些参数就可以作为变量传入。
import numpy as np
from scipy.integrate import quad
def integrand(x, a, b, c):
    return a*x*x+b*x+c
a=3
b=4
c=1
result=quad(integrand, 0, np.inf, args=(a,b,c))
print(result)
3）二重积分和三重积分可以用dblquad和tplquad函数来实现。
import numpy as np
from scipy.integrate import quad,dblquad
def integrand1(t,x,n):
    return np.exp(-x*t)/t**n
n=4
result=dblquad(lambda t,x:integrand1(t,x,n),0,np.inf,lambda x:0,lambda x:np.inf)
print(result)
4）定积分可以使用fixed_quad实现。
import numpy as np
from scipy.integrate import quadrature,dblquad,fixed_quad
def integrand(x,a,b):
    return a*x+b
a=2
b=1
fixed_result=fixed_quad(integrand,0,1,args=(a,b))
result=quadrature(integrand,0,1,args=(a,b))
print(fixed_result)
print(result)
5）对于带有任意间隔采样的函数的积分，可以使用simps。辛普森法则可以估计三个相邻点间的函数为抛物线函数：
import numpy as np
from scipy.integrate import simps
def func1(a,x):
    return a*x**2+2
def func2(b,x):
    return b*x**3+4
x=np.array([1,2,4,5,6])
y1=func1(2,x)
intgrl1=simps(y1,x)
print(intgrl1)
y2=func2(3,x)
intgrl2=simps(y2,x)
print(intgrl2)

5. 信号处理scipy.signal：

信号处理工具箱包含卷积、B样条插值、滤波、滤波器设计、Matlab式IIR滤波器设计、连续时间的线性系统、离散时间的线性系统、线性不变系统、信号波形、窗函数、小波分析、谱峰的找寻、光谱分析等函数。
1）设计滤波器，使用signal子模块，带通滤波器使用iirdesign函数。格式：
iirdesign(wp, ws, gpass, gstop, analog=False, ftype='ellip', output='ba')
其中，wp为通带边缘频率；ws为阻带边缘频率；gpass为通带最大增益，单位dB；gstop为阻带最小衰减，单位dB；参数analog为True返回一个模拟滤波器，为False返回数字滤波器；ftype表示滤波器类型的字符串；output表示输出类型的字符串。ftype可选值有：

output的类型参数有ba和zpk，ba表示分子/分母型，zpk表示极点-零点型。
设计好滤波器后，调用SciPy的signal子模块中的lfilter函数进行滤波。格式：
lfilter(b, a, axis=-1, zi=None)
其中，b为滤波器的分子系数向量；a为滤波器的分母系数向量；x为需要滤波的信号向量；axis为输入向量的滤波轴线；zi为滤波器延迟初始条件。
示例：
b,a=signal.iirdesign(0.05,0.2,2,10)     # 设计带通滤波器
z=signal.lfilter(b,a,y)                # 利用滤波器对信号y滤波
代码中设计的带通滤波器，通带截止频率为0.05*f0，阻带的起始频率为0.2*f0，通带最大增益2dB，阻带最小衰减为10dB。
detrend函数是一个滤波函数，用于去除矩阵或者向量的均值或者线性趋势。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
t=np.linspace(0,5,100)
x=t+np.random.normal(size=100)
plt.plot(t,x,linewidth=3)
plt.show()
plt.plot(t,signal.detrend(x),linewidth=3)
plt.show()
2）cspline2d函数将一个带有镜像对称边界的二维FIR滤波器应用于样条函数，而convolve2d函数对任意二维滤波器卷积并允许选择镜像对称边界。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal,misc
img=misc.face(gray=True)
splineresult=signal.cspline2d(img,2.0)
arr1=np.array([[-1,0,1],[-2,0,2],[-1,0,1]],dtype=np.float32)
derivative=signal.convolve2d(splineresult,arr1,boundary='symm',mode='same')
plt.figure()
plt.imshow(derivative)
plt.title("image filtered by spline edge filter")
plt.gray()
plt.show()

6. 傅里叶变换scipy.fftpack：

对实数或复数序列的离散傅里叶变换和逆变换可以分别用fft和ifft函数来计算，fft是指快速傅里叶变换。
import numpy as np
from scipy.fftpack import fft,ifft
x=np.random.random_sample(5)
y=fft(x)
print(y)
yinv=ifft(y)
print(yinv)

7. 空间数据结构和算法scipy.spatial：

空间数据指的是地理空间或垂直空间相关的数据对象或元素，这种数据包括点、线、多边形、其他几何和地理特性信息。几何和地理特征信息可以映射为位置信息，用于跟踪或定位各种装置。
KD树是一种空间划分数据结构，它将所有点构建在k维空间中。在数学上，空间划分是将一个空间划分为多个相邻空间的过程，它将空间分成不重叠的区域，空间中的每个点只属于一个区域。scipy可以实现快速近邻查找算法，还可以对初始向量集合进行距离矩阵的计算。
Scipy有支持各种空间计算功能的模块，可以计算Delaunay三角剖分、Voronoi图、N维凸包，还有画图工具，可以将这些计算结果画在二维空间中。三角剖分的示例：
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
arr_pt=np.array([[0,0],[0,1.1],[1,0],[1,1]])
arr1=np.array([0.,0.,1.,1.])
arr2=np.array([0.,1.1,0.,1.])
triangle_result=Delaunay(arr_pt)
plt.triplot(arr1,arr2,triangle_result.simplices.copy())
plt.show()
plt.plot(arr1,arr2,'ro')
plt.show()
凸包包括给定集合中的所有点，可以用ConvexHull函数来计算：
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
randpoints=np.random.rand(25,2)
hull=ConvexHull(randpoints)
plt.plot(randpoints[:,0],randpoints[:,1],'x')
for simplex in hull.simplices:
    plt.plot(randpoints[simplex,0],randpoints[simplex,1],'k')
plt.show()
可以用KDTree来查找点集中离选定点最近的点。
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import KDTree
vertx=np.array([[1,1],[1,2],[1,3],[2,1],[2,2],[2,3],[3,1],[3,2],[3,3]])
tree=KDTree(vertx)
tree.query([1.1,1.1])
x_vals=np.linspace(0.5,3.5,31)
y_vals=np.linspace(0.5,3.5,33)
xgrid,ygrid=np.meshgrid(x_vals,y_vals)
xy=np.c_[xgrid.ravel(),ygrid.ravel()]
plt.pcolor(x_vals,y_vals,tree.query(xy)[1].reshape(33,31))
plt.plot(vertx[:,0],vertx[:,1],'ko')
plt.show()

8. 最优化scipy.optimize:

最优化是查找带一个或多个变量的对象函数，在多个预先定义变量约束条件下的最佳方案的过程。对象函数被当作代价函数以最小化，或者被当作利用函数或利益函数以最大化。
scipy.optimize包提供了标量以及多维函数最小化、曲线拟合和根求解的最有用函数。BFGS算法（Broyden- Fletcher-Goldfarb-Shanno）用目标函数的梯度来快速收敛求解。示例：
import numpy as np
from scipy.optimize import minimize
def rosenbrock(x):
    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0+(1-x[:-1])**2.0)
x0=np.array([1.3,0.7,0.8,1.9,1.2])
def rosen_derivative(x):
    x1=x[1:-1]
    x1_m1=x[:-2]
    x1_p1=x[2:]
    derivative=np.zeros_like(x)
    derivative[1:-1]=200*(x1-x1_m1**2)-400*(x1_p1-x1**2)*x1-2*(1-x1)
    derivative[0]=-400*x[0]*(x[1]-x[0]**2)-2*(1-x[0])
    derivative[-1]=200*(x[-1]-x[-2]**2)
    return derivative
res=minimize(rosenbrock,x0,method='BFGS',jac=rosen_derivative,options={'disp':True})
minizime函数也有一个对多个约束条件最小化算法的接口，一种是SLSQP算法（Sequential Least Square Programming optimization）。示例：
import numpy as np
from scipy.optimize import minimize
def func(x,sign=1.0):
    return sign*(2*x[0]*x[1]+2*x[0]-x[0]**2-2*x[1]**2)
def func_deriv(x,sign=1.0):
    dfdx0=sign*(-2*x[0]+2*x[1]+2)
    dfdx1=sign*(2*x[0]-4*x[1])
    return np.array([dfdx0,dfdx1])
cons=({'type':'eq',
       'fun':lambda x:np.array([x[0]**3-x[1]]),
       'jac':lambda x:np.array([3.0*(x[0]**2.0),-1.0])},
      {'type':'ineq',
       'fun':lambda x:np.array([x[1]-1]),
       'jac':lambda x:np.array([0.0,1.0])})
res=minimize(func,[-1.0,1.0],args=(-1.0,),jac=func_deriv,method='SLSQP',options={'disp':True})
print(res.x)
res=minimize(func,[-1.0,1.0],args=(-1.0,),jac=func_deriv,
             constraints=cons,method='SLSQP',options={'disp':True})
print(res.x)
scipy.optimize还能实现寻找全局最小和局部最小，使用BFGF算法寻找局部最小随网格大小的增加会变慢，对标量函数最好使用Brent算法。
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
def f(x):
    return x**2+10*np.sin(x)
x=np.arange(-10,10,0.1)
plt.plot(x,f(x))
plt.show()
optimize.fmin_bfgs(f,0)
grid=(-10,10,0.1)
optimize.brute(f,(grid,))
optimize.brent(f)
optimize.fminbound(f,0,10)
展示约束最优化的应用：
import numpy as np
from scipy import optimize
def f(x):
    return np.sqrt((x[0]-2)**2+(x[1]-3)**2)
def constraint(x):
    return np.atleast_1d(2.5-np.sum(np.abs(x)))
optimize.fmin_slsqp(f,np.array([0,2]),ieqcons=[constraint,])
optimize.fmin_cobyla(f,np.array([3,4]),cons=constraint)
有多种方法可以求解多项式的根，二分法求解多项式的根方法示例：
import numpy as np
from scipy import optimize
def polynomial_func(x):
    return np.cos(x)**3+4-2*x
print(optimize.bisect(polynomial_func,1,5))
还可以使用牛顿-拉斐森法求解多项式的根：
import scipy
from scipy import optimize
def polynomial_func(x):
    y=x+2*scipy.cos(x)
    return y
print(optimize.newton(polynomial_func,1))
在数学上，拉格朗日乘子法被用于寻找函数在等量约束的条件下的局部最小和局部最大。
import numpy as np
from scipy.optimize import fsolve
def func_orig(data):
    xval=data[0]
    yval=data[1]
    Multiplier=data[2]
    return xval+yval+Multiplier*(xval++2+yval**2-1)
def deriv_func_orig(data):
    dLambda=np.zeros(len(data))
    step_size=1e-3
    for i in range(len(data)):
        ddata=np.zeros(len(data))
        ddata[i]=step_size
        dLambda[i]=(func_orig(data+ddata)-func_orig(data-ddata))/(2*step_size)
    return dLambda
data1=fsolve(deriv_func_orig,[1,1,0])
print(data1,func_orig(data1))
data2=fsolve(deriv_func_orig,[-1,-1,0])
print(data2,func_orig(data2))

9. 插值scipy.interpolate：

插值，是在给定的已知的离散值集合范围内寻找新数据点的方法。interpolate子包包含插值运算的各种方法，支持用spline函数、univariate函数和multivariate函数进行一维和多维插值，还有拉格朗日和泰勒多项式插值。对于FITPACK和DFITPACK函数，还有封装类。
用线性和立方插值进行一维插值的示例：
import numpy as np
from scipy.interpolate import interp1d
x_val=np.linspace(0,20,10)
y_val=np.cos(-x_val**2/8.0)
f=interp1d(x_val,y_val)
f2=interp1d(x_val,y_val,kind='cubic')
xnew=np.linspace(0,20,25)
import matplotlib.pyplot as plt
plt.plot(x_val,y_val,'o',xnew,f(xnew),'-',xnew,f2(xnew),'--')
plt.legend(['data','linear','cubic'],loc='best')
plt.show()
用griddata函数可对超过150个点的多元数据进行插值，生成四个子图：
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
def func_user(x,y):
    return x*(1-x)*np.cos(4*np.pi*x)*np.sin(4*np.pi*y**2)**2
x,y=np.mgrid[0:1:100j,0:1:200j]
points=np.random.rand(150,2)
values=func_user(points[:,0],points[:,1])
grid_z0=griddata(points,values,(x,y),method='nearest')
grid_z1=griddata(points,values,(x,y),method='linear')
grid_z2=griddata(points,values,(x,y),method='cubic')
f,axarr=plt.subplots(2,2)
axarr[0,0].imshow(func_user(x,y).T,extent=(0,1,0,1),origin='lower')
axarr[0,0].plot(points[:,0],points[:,1],'k',ms=1)
axarr[0,0].set_title('Original')
axarr[0,1].imshow(grid_z0.T,extent=(0,1,0,1),origin='lower')
axarr[0,1].set_title('Nearest')
axarr[1,0].imshow(grid_z1.T,extent=(0,1,0,1),origin='lower')
axarr[1,0].set_title('Linear')
axarr[1,1].imshow(grid_z2.T,extent=(0,1,0,1),origin='lower')
axarr[1,1].set_title('Cubic')
plt.show()

10. 线性代数scipy.linalg：

SciPy线性代数将变量当作一个对象，这个对象可以转换成一个二维数组并返回一个二维数组。计算矩阵的逆的示例：
import numpy as np
from scipy import linalg
A=np.array([[2,3],[4,5]])
linalg.inv(A)
B=np.array([[3,8]])
print(A*B)
print(A.dot(B.T))
计算矩阵的逆和行列式的值：
import numpy as np
from scipy import linalg
A=np.array([[2,3],[4,5]])
print(linalg.inv(A))
print(linalg.det(A))
可以通过逆矩阵方法速求解线性方程组，也可以使用NumPy的求解线程方程组的方法：
import numpy as np
from scipy import linalg
A=np.array([[2,3],[4,5]])
B=np.array([[5],[6]])
print(linalg.inv(A).dot(B))
print(np.linalg.solve(A,B))
求解一组标量系数并用模型进行数据拟合，其中用到lstsq方法，用来寻找线性矩阵方程的最小平方解。
import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt
coeff_1,coeff_2=5.0,2.0
i=np.r_[1:11]
x=0.1*i
y=coeff_1*np.exp(-x)+coeff_2*x
z=y+0.05*np.max(y)*np.random.randn(len(y))
A=np.c_[np.exp(-x)[:,np.newaxis],x[:,np.newaxis]]
coeff,resid,rank,sigma=linalg.lstsq(A,z)
x2=np.r_[0.1:1.0:100j]
y2=coeff[0]*np.exp(-x2)+coeff[1]*x2
plt.plot(x,z,'x',x2,y2)
plt.axis([0,1,3.0,5.5])
plt.title('Data fitting with linalg.lstsq')
plt.show()
还可以使用linag.svg和linag.diagsvg函数进行奇异值分解：
import numpy as np
from scipy import linalg
A=np.array([[5,4,2],[4,8,7]])
row=2
col=3
U,s,Vh=linalg.svd(A)
Sig=linalg.diagsvd(s,row,col)
U,Vh=U,Vh
print(U)
print(Sig)
print(Vh)
用ARPACK解决稀疏特征值问题：
import numpy as np
from scipy.linalg import eigh
from scipy.sparse.linalg import eigsh
np.set_printoptions(suppress=True)
np.random.seed(0)
random_matrix=np.random.random((75,5))-0.5
random_matrix=np.dot(random_matrix,random_matrix.T)
eigenvalues_all,eigenvectors_all=eigh(random_matrix)
eigenvalues_large,eigenvectors_large=eigsh(random_matrix,3,which='LM')
print(eigenvalues_all[-3:])
print(eigenvalues_large)
print(np.dot(eigenvectors_large.T,eigenvectors_all[:,-3:]))
如果要找最小特征值，可以使用eigenvalues_small, eigenvectors_small =eigsh(random_matrix , 3, which='SM')，如果返回一个没有收敛的错误，在参数中增加传入容限值tol=1E-2，但可能损失精度。也可以增加最大迭代次数，将maxiter=5000传入；还有是用shift-inter模式，加参数sigma=0或2和使用which='LM"。

11. 统计学scipy.stats：

SciPy有很多统计函数用于数据运算，还有用于掩码数组的版本。使用离散二项式随机变量并画出该概率质量函数：
import numpy as np
from scipy.stats import binom
import matplotlib.pyplot as plt
n,p=5,0.4
mean,variance,skewness,kurtosis=binom.stats(n,p,moments='mvsk')
x_var=np.arange(binom.ppf(0.01,n,p),binom.ppf(0.99,n,p))
plt.plot(x_var,binom.pmf(x_var,n,p),'ro',ms=5,label='PMF of binamial')
plt.vlines(x_var,0,binom.pmf(x_var,n,p),colors='r',lw=3,alpha=0.5)
plt.show()
其中，binom.pmf(k)=choose(n,k)*p**k*(1-p)**(n-k)，k取值(0,1,...,n)，n和p是形状参数。
用离散随机变量并画出概率质量函数：
import numpy as np
from scipy.stats import geom
import matplotlib.pyplot as plt
p=0.5
mean,variance,skewness,kurtosis=geom.stats(p,moments='mvsk')
x_var=np.arange(geom.ppf(0.01,p),geom.ppf(0.99,p))
plt.plot(x_var,geom.pmf(x_var,p),'go',ms=5,label='PMF of geomatric')
plt.vlines(x_var,0,geom.pmf(x_var,p),colors='g',lw=3,alpha=0.5)
plt.show()
其中，geom.pmk(k)=(1-p)**(k-1)*p，k≥1，p是形状参数。
广义帕累托连续随机变量的计算并画出其概率密度函数：
import numpy as np
from scipy.stats import genpareto
import matplotlib.pyplot as plt
c=0.1
mean,variance,skewness,kurtosis=genpareto.stats(c,moments='mvsk')
x_var=np.linspace(genpareto.ppf(0.01,c),genpareto.ppf(0.99,c))
plt.plot(x_var,genpareto.pdf(x_var,c),'b-',lw=3,alpha=0.6,label='PDF of Generic Pareto')
plt.show()
其中，genpareto.pdf(x,c)=(1+c*x)**(-1-1/c)，如果c≥0，则x≥0；如果c<0，0≤x≤-1/c。
广义伽马连续型随机变量的计算并画出其概率密度函数：
import numpy as np
from scipy.stats import gengamma
import matplotlib.pyplot as plt
a,c=4.41623854294,3.11930916792
mean,variance,skewness,kurtosis=gengamma.stats(a,c,moments='mvsk')
x_var=np.linspace(gengamma.ppf(0.01,a,c),gengamma.ppf(0.99,a,c),100)
plt.plot(x_var,gengamma.pdf(x_var,a,c),'b-',lw=3,alpha=0.6,label='PDF of Generic Gamma')
plt.show()
其中，gengamma.pdf(x,a,c)=abs(c)*x**(c*a-1)*exp(-x**c)/gamma(a)，这里rx>0，a>0，c!=0，a和c为形状参数。
多元正态随机变量的计算并画出其概率密度函数：
import numpy as np
from scipy.stats import multivariate_normal
import matplotlib.pyplot as plt
x_var=np.linspace(5,25,20,endpoint=False)
y_var=multivariate_normal.pdf(x_var,mean=10,cov=2.5)
plt.plot(x_var,y_var)
plt.show()

12. 聚类scipy.cluster：

聚类，是将一个大的对象集合分成多个组的过程，它用一些参数 来判定一个对象和改组内的对象更相似，而和其他类或组中的对象不相似。
SciPy聚类包有两个模块，向量量化VQ和层次聚类。VQ模块支持K-means聚类和向量量化，层次聚类支持分层聚类和聚合聚类。
·向量量化：VQ是通过原型向量的分布来构建概率密度模型的信号处理技术，它通过将一个大的向量或点的集合分成多个组来完成建模，这些组在其近邻包括几乎相同数量的点，每个组都用中心点来表示。
·K-means：是按照与聚类中心聚类最近的方式，把n个观察点划分为k个聚类
from scipy.cluster.vq import kmeans,vq
from numpy.random import rand
from numpy import vstack,array
import matplotlib.pyplot as plt
data_set=vstack((rand(200,2)+array([.5,.5]),rand(200,2)))
centroids_of_clusters,_=kmeans(data_set,2)
index,_=vq(data_set,centroids_of_clusters)
plt.plot(data_set[index==0,0],data_set[index==0,1],'ob',data_set[index==1,0],data_set[index==1,1],'or')
plt.plot(centroids_of_clusters[:,0],centroids_of_clusters[:,1],'sg',markersize=8)
plt.show()
层次聚类：利用观察值建立一组聚类的层次，结果最终用树状图表示。可以分为两种：
·分离聚类：生成聚类的自上而下的方法，开始于最顶层的类，随着向下移动不断分裂
·聚合聚类：这是一个自下而上的方法，每个观察值是一个类，该方法在向上移动过程中将类不断进行配对
层次聚类包括很多函数，如将层次聚类切分为扁平聚类的函数、聚合聚类函数、聚类可视化函数、数据结构，以及将层次聚类结果表示为树状图的函数、计算聚类层次的函数、检查链接和不规则指标有效性的推断函数，等等。
import numpy
from numpy.random import rand
from scipy.spatial.distance import pdist
from scipy.cluster import hierarchy as sch
import matplotlib.pyplot as plt
x=rand(8,80)
x[0:4,:]*=2
y=pdist(x)
z=sch.linkage(y)
sch.dendrogram(z)
plt.show()

13. 曲线拟合scipy.optimize.curve_fit：

构造一个与数据序列接近的数学函数或曲线的过程为曲线拟合，通常曲线拟合必须满足一些约束条件。曲线拟合也可以用来观察多个变量之间的关系，可以用不同类型的曲线做曲线拟合，如直线、多项式、圆锥曲线、三角函数、圆、椭圆等。
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
xdata=np.random.uniform(0.,50.,80)
ydata=3.*xdata+2.+np.random.normal(0.,5.,80)
plt.plot(xdata,ydata,'.')
def line_func(x,a,b):
    return a*x+b
opt_param,cov_estimated=curve_fit(line_func,xdata,ydata)
errors=np.repeat(5.,80)
plt.errorbar(xdata,ydata,yerr=errors,fmt="none")
opt_param,cov_estimated=curve_fit(line_func,xdata,ydata,sigma=errors)
print('a=',opt_param[0],'+/-',cov_estimated[0,0]**0.5)
print('b=',opt_param[1],'+/-',cov_estimated[1,1]**0.5)
plt.errorbar(xdata,ydata,yerr=errors,fmt="none")
xnew=np.linspace(0.,50.,80)
plt.plot(xnew,line_func(xnew,opt_param[0],opt_param[1]),'r-')
plt.errorbar(xdata,ydata,yerr=errors,fmt="none")
plt.plot(xnew,line_func(xnew,*opt_param),'r-')
plt.show()
代码中先生成一些带噪声的随机数据，然后定义一个函数来表示该模型并做曲线拟合，接下来求解参数a和b，最后画出误差。也可以使用三角函数做曲线拟合，如cos函数：
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
np.random.seed(0)
def f(x,omega,p):
    return np.cos(omega*x+p)
x=np.linspace(0,10,100)
y=f(x,2.5,3)+.1*np.random.normal(size=100)
params,params_cov=curve_fit(f,x,y)
t=np.linspace(0,10,500)
plt.figure(1)
plt.clf()
plt.plot(x,y,'bx')
plt.plot(t,f(t,*params),'r-')
plt.show()

14. 文件输入/输出scipy.io:

SciPy提供了一系列文件格式的读和写的支持，包括MATLAB文件、ALD文件、Matrix Market文件、无格式的FORTRAN文件、WAV声音文件、ARFF文件、NetCDF文件。
网络通用数据格式NetCDF(network Common Data Form)是由美国大学大气研究协会UCAR (University Corporation for Atmospheric Research)的Unidata项目科学家针对科学数据的特点开发的，是一种面向数组型并适于网络共享的数据的描述和编码标准。目前，NetCDF广泛应用于大气科学、水文、海洋学、环境模拟、地球物理等诸多领域。SciPy可以读取这种格式的文件：
import numpy as np
from scipy.io import netcdf
# create a file
f=netcdf.netcdf_file('TextFile.nc','w')  
f.history='Test netCDF File Creation'
f.createDimension('age',12)
age=f.createVariable('age','i',('age',))
age.units='Age of persons in Years'
age[:]=np.arange(12)
f.close()
# read a file
f=netcdf.netcdf_file('TextFile.nc','r')
print(f.history)
age=f.variables['age']
print(age.units)
print(age.shape)
print(age[-1])
f.close()
WEKA机器学习工具将文件存储为ARFF格式，这是一种可能包含数值、字符串和数据值的文本文件。打开读取方法代码：
from scipy.io import arff
file1=open('test.arff')
data,meta=arff.loadarff(file1)
print(data)
print(meta)