1.引入
scipy是一个利用NumPy数组和操作来处理科学家和工程师通常面临的标准问题的包：集成、确定一个函数的最大值或最小值、寻找大稀疏矩阵的特征向量、检验两个分布是否相同等等。从优化和数据拟合开始，因为这些是一些最常见的任务，然后通过插值、集成、空间分析、聚类、信号和图像处理、稀疏矩阵和统计。
2.最优化和最小化
线性回归，找到一个函数的最小值和最大值，确定一个函数的根，并找到两个函数相交的位置。
2.1数据建模和拟合
有几种方法可以用线性回归来拟合数据。在本节中，我们将使用curve_fit，这是一种基于χ2的方法（换句话说，是一种最佳拟合方法）。在下面的示例中，我们从一个带噪声的已知函数生成数据，然后用curve_fit拟合有噪声数据。我们将在这个例子中建模的函数是一个简单的线性方程，f(x)=ax+b。
import numpy as np
from scipy.optimize import curve_fit
# Creating a function to model and create data
def func(x, a, b):
return a*x+b
# Generating clean data
x = np.linspace(0, 10, 100)
y = func(x, 1, 2)
# Adding noise to the data
yn = y + 0.9 * np.random.normal(size=len(x))
# Executing curve_fit on noisy data
popt, pcov = curve_fit(func, x, yn)
# popt returns the best fit values for parameters of
# the given model (func).
print(popt)
对a,b拟合的很接近
最小二乘法去拟合高斯分布，一个非线性函数
# Creating a function to model and create data
def func(x, a, b, c):
return a*np.exp(-(x-b)**2/(2*c**2))
# Generating clean data
x = np.linspace(0, 10, 100)
y = func(x, 1, 5, 2)
# Adding noise to the data
yn = y + 0.2 * np.random.normal(size=len(x))
# Executing curve_fit on noisy data
popt, pcov = curve_fit(func, x, yn)
可以拟合一个具有多个高斯剖面的一维数据集
# Two-Gaussian model
def func(x, a0, b0, c0, a1, b1,c1):
return a0*np.exp(-(x - b0) ** 2/(2 * c0 ** 2)) + a1 * np.exp(-(x - b1) ** 2/(2 * c1 ** 2))
# Generating clean data
x = np.linspace(0, 20, 200)
y = func(x, 1, 3, 1, -2, 15, 0.5)
# Adding noise to the data
yn = y + 0.2 * np.random.normal(size=len(x))
# Since we are fitting a more complex function,
# providing guesses for the fitting will lead to
# better results.
guesses = [1, 3, 1, 1, 15, 1]
# Executing curve_fit on noisy data
popt, pcov = curve_fit(func, x, yn, p0=guesses)
2.2对功能的解决方案拟合
使用scipy.optimize.fsolve找到方程的根
from scipy.optimize import fsolve
import numpy as np
line = lambda x:x+3
solution = fsolve(line, -2)
print(solution)
array([-3.])
from scipy.optimize import fsolve
import numpy as np
# Defining function to simplify intersection solution
def findIntersection(func1, func2, x0):
return fsolve(lambda x : func1(x) - func2(x), x0)
# Defining functions that will intersect
funky = lambda x : np.cos(x / 5) * np.sin(x / 2)
line = lambda x : 0.01 * x - 0.5
# Defining range and getting solutions on intersection points
x = np.linspace(0,45,10000)
result = findIntersection(funky, line, [15, 20, 30, 35, 40, 45])
# Printing out results for x and y
print(result, line(result))
[13.40773078 18.11366128 31.78330863 37.0799992 39.84837786 43.8258775 ] [-0.36592269 -0.31886339 -0.18216691 -0.12920001 -0.10151622 -0.06174122]
3. 插值
包含信息的数据通常具有一种功能形式，作为分析人员，希望对其进行建模。给定一组样本数据，获取点之间的中间值有助于理解和预测数据在非采样域中将做什么。科学研究所提供了十几种不同的插值功能，从简单的单变量情况到复杂的多变量情况。当采样数据可能由一个自变量引导时，则使用单变量插值，而多变量插值假设有多个自变量。
插值有两种基本的插值方法：（1)将一个函数匹配到整个数据集，或者(2）适应数据集的不同部分，其中每个函数的连接被平滑地连接。第二种类型被称为样条曲线插值，当数据的函数形式很复杂时，它可以是一个非常强大的工具。
我们将首先展示如何插值一个简单的函数，然后继续到一个更复杂的情况。下面的示例使用具有不同拟合参数的scipy.interpolate.interp1d插值一个正弦函数。第一个参数是“线性”拟合，第二个参数是“二次”拟合。
import numpy as np
from scipy.interpolate import interp1d
# Setting up fake data
x = np.linspace(0, 10 * np.pi, 20)
y = np.cos(x)
# Interpolating data
fl = interp1d(x, y, kind='linear')
fq = interp1d(x, y, kind='quadratic')
# x.min and x.max are used to make sure we do not
# go beyond the boundaries of the data for the
# interpolation.
xint = np.linspace(x.min(), x.max(), 1000)
yintl = fl(xint)
yintq = fq(xint)
使用一个叫做scipy.interpolate.UnivariateSpline的样条拟合函数插值有噪声的数据。
import numpy as np
import matplotlib.pyplot as mpl
from scipy.interpolate import UnivariateSpline
# Setting up fake data with artificial noise
sample = 30
x = np.linspace(1, 10 * np.pi, sample)
y = np.cos(x) + np.log10(x) + np.random.randn(sample) / 10
# Interpolating the data
f = UnivariateSpline(x, y, s=1)
#选项s是平滑系数，在与噪声匹配数据时应使用平滑系数。如果是s=0，则插值将通过所有点，同时忽略噪声
# x.min and x.max are used to make sure we do not
# go beyond the boundaries of the data for the
# interpolation.
xint = np.linspace(x.min(), x.max(), 1000)
yint = f(xint)
scipy.interpolate.griddata函数用于其处理非结构化N维数据的能力。例如，如果您有一个1000×1000像素的图像，然后随机选择1000点，则可以重建图像
import numpy as np
from scipy.interpolate import griddata
# Defining a function
ripple = lambda x, y: np.sqrt(x**2 + y**2)+np.sin(x**2 + y**2)
# Generating gridded data. The complex number defines
# how many steps the grid data should have. Without the
# complex number mgrid would only create a grid data structure
# with 5 steps.
grid_x, grid_y = np.mgrid[0:5:1000j, 0:5:1000j]
# Generating sample that interpolation function will see
xy = np.random.rand(1000, 2)
sample = ripple(xy[:,0]*5, xy[:,1] * 5)
# Interpolating data with a cubic
grid_z0 = griddata(xy * 5, sample, (grid_x, grid_y), method='cubic')
4.积分
积分是数学和科学中的一个关键工具，因为微分和积分是微积分的两个关键组成部分。给定来自函数或数据集的曲线，我们可以计算它下面的面积。在传统的课堂设置中，我们将解析地整合一个函数，但研究环境中的数据很少以这种形式给出，我们需要近似它的明确积分。
scipy有一系列不同的函数整合积分方程式和数据。我们将首先讨论这些函数，然后转向数据解决方案。之后，我们将使用我们之前使用的数据拟合工具来计算定积分解。
4.1分析积分
使用下面所表达的函数。集成起来很简单，而且其解决方案的估计误差也很小。
import numpy as np
from scipy.integrate import quad
# Defining function to integrate
func = lambda x: np.cos(np.exp(x)) ** 2
# Integrating function with upper and lower
# limits of 0 and 3, respectively
solution = quad(func, 0, 3)
print solution
# The first element is the desired value
# and the second is the error.
# (1.296467785724373, 1.397797186265988e-09)
4.2 数值积分
import numpy as np
from scipy.integrate import quad, trapz
# Setting up fake data
x = np.sort(np.random.randn(150)*4+ 4).clip(0,5)
func = lambda x: np.sin(x) * np.cos(x ** 2) + 1
y = func(x)
# Integrating function with upper and lower
# limits of 0 and 5, respectively
fsolution = quad(func, 0, 5)
dsolution = trapz(y, x=x)
print('fsolution='+ str(fsolution[0]))
print('dsolution='+ str(dsolution))
print('The difference is ' + str(np.abs(fsolution[0] - dsolution)))
# fsolution = 5.10034506754
# dsolution = 5.04201628314
# The difference is 0.0583287843989.
四元积分器只能使用一个可调用的函数，而trapz是一个利用数据点的数值积分器。
5.数据分析
Numpy有基本的分析函数，像mean,std,median,argmax,argmin,此外，numpy.arrays有内置的方法，允许我们轻松使用Numpy 分析
import numpy as np
# Constructing a random array with 1000 elements
x = np.random.randn(1000)
# Calculating several of the built-in methods
# that numpy.array has
mean = x.mean()
std = x.std()
var = x.var()
4.1连续分布和离散分布
scipy提供的丰富的函数与numpy.random函数对比，numpy.random局限于均匀分布和高斯分布。
当我们从统计数据中调用一个分布时，我们可以通过几种方式提取其信息：概率密度函数、累积分布函数、随机变量样本、百分比点函数等等。
常见连续分布：
下面代码展示了如何接近如下分布
import numpy as np
import scipy.stats import norm
# Set up the sample range
x = np.linspace(-5,5,1000)
# Here set up the parameters for the normal distribution,
# where loc is the mean and scale is the standard deviation.
dist = norm(loc=0, scale=1)
# Retrieving norm's PDF and CDF
pdf = dist.pdf(x)
cdf = dist.cdf(x)
# Here we draw out 500 random values from the norm.
sample = dist.rvs(500)
连续分布可以被loc和scale不同的取值设置。
在其他情况下，人们将需要一个离散分布，如泊松，二项式，或几何。与连续分布不同，离散分布对于给定数量的事件发生在固定的时间/空间间隔内，事件以已知的平均速率和每个事件无关的事件是有用的问题。
下图是概率质量密度的集合分布
import numpy as np
from scipy.stats import geom
# Here set up the parameters for the geometric distribution.
p = 0.5
dist = geom(p)
# Set up the sample range.
x = np.linspace(0, 5, 1000)
# Retrieving geom's PMF and CDF
pmf = dist.pmf(x)
cdf = dist.cdf(x)
# Here we draw out 500 random values.
sample = dist.rvs(500)
4.2函数
scipy有60多个统计功能，如果你只是好奇有什么是可用的功能，可能会难以消化。思考统计函数的最佳方法是，它们可以分别描述或检验样本，例如，特定值的频率或 Kolmogorov-Smirnov test(KS检验)。在stats包中，有许多函数如kstest和normaltest用与检验样本。这些分布检验非常有助于确定样本是否来自特定分布。
import numpy as np
from scipy import stats
# Generating a normal distribution sample
# with 100 elements
sample = np.random.randn(100)
# normaltest tests the null hypothesis.
out = stats.normaltest(sample)
print('normaltest output')
print('Z-score='+ str(out[0]))
print('P-value='+ str(out[1]))
# kstest is the Kolmogorov-Smirnov test for goodness of fit.
# Here its sample is being tested against the normal distribution.
# D is the KS statistic and the closer it is to 0 the better.
out = stats.kstest(sample, 'norm')
print('\nkstest output for the Normal distribution')
print('D='+ str(out[0]))
print('P-value='+ str(out[1]))
# Similarly, this can be easily tested against other distributions,
# like the Wald distribution.
out = stats.kstest(sample, 'wald')
print('\nkstest output for the Wald distribution')
print('D='+ str(out[0]))
print('P-value='+ str(out[1]))
研究人员通常使用描述性函数来进行统计学。统计包中的一些描述函数包括几何平均值（平均值）、样本偏差（偏斜）和样本中值的频率（项）。
import numpy as np
from scipy import stats
# Generating a normal distribution sample
# with 100 elements
sample = np.random.randn(100)
# The harmonic mean: Sample values have to
# be greater than 0.
out = stats.hmean(sample[sample > 0])
print('Harmonic mean='+ str(out))
# The mean, where values below -1 and above 1 are
# removed for the mean calculation
out = stats.tmean(sample, limits=(-1, 1))
print('\nTrimmed mean='+ str(out))
# Calculating the skewness of the sample偏度
out = stats.skew(sample)
print('\nSkewness='+ str(out))
# Additionally, there is a handy summary function called
# describe, which gives a quick look at the data.
out = stats.describe(sample)
print('\nSize='+ str(out[0]))
print('Min='+ str(out[1][0]))
print('Max='+ str(out[1][1]))
print('Mean='+ str(out[2]))
print('Variance='+ str(out[3]))
print('Skewness='+ str(out[4]))
print('Kurtosis='+ str(out[5]))
5.空间和聚类分析
从生物科学到天体物理科学，空间和聚类分析是识别模式、星群和星系团的关键。例如，在生物学中，不同植物物种的间隔暗示了种子是如何分散、与环境相互作用和生长的。在天体物理学中，这些分析技术被用来寻找和识别星团、星系团和大尺度的细丝（由星系团组成）。在计算机科学领域，识别和映射复杂的节点和信息网络本身是一项至关重要的研究。随着大数据和数据挖掘，识别数据集群正变得越来越重要，以便组织已发现的信息，而不是被它淹没。
scipy提供了空间分析类（scipy.spatial）和集群分析类（scipy.cluster）。空间类包括分析数据点之间距离的函数(例如，k-d树)。集群类提供了两个主要的子类：矢量量化(vq)和层次聚类（层次结构）。
矢量量化集群群，大的数据点（向量）组成的数据集，其中每个组由质心表示。层次结构子类包含构造集群和分析其子结构的函数。
5.1矢量量化
向量量化是一个可以与信号处理、数据压缩和聚类相关联的一般术语。在这里，我们将关注集群组件，从如何将数据到vq包以识别集群。
import numpy as np
from scipy.cluster import vq
# Creating data
c1 = np.random.randn(100, 2) + 5
c2 = np.random.randn(30, 2) - 5
c3 = np.random.randn(50, 2)
# Pooling all the data into one 180 x 2 array
data = np.vstack([c1, c2, c3])
# Calculating the cluster centroids and variance
# from kmeans
centroids, variance = vq.kmeans(data, 3)
# The identified variable contains the information
# we need to separate the points in clusters
# based on the vq function.
identified, distance = vq.vq(data, centroids)
# Retrieving coordinates for points in each vq
# identified core
vqc1 = data[identified == 0]
vqc2 = data[identified == 1]
vqc3 = data[identified == 2]
5.2层次聚类
层次结构聚类是一个强大的工具，用于识别嵌套在较大结构中的结构。但是使用输出可能是很棘手的，因为我们不能像使用kmans技术那样获得干净识别的集群。下面是一个例子9，其中我们生成了一个包含多个集群的系统。若要使用层次结构函数，我们建立了一个距离矩阵，并且输出是一个树状图树。
import numpy as np
import matplotlib.pyplot as mpl
from mpl_toolkits.mplot3d import Axes3D
from scipy.spatial.distance import pdist, squareform
import scipy.cluster.hierarchy as hy
# Creating a cluster of clusters function
def clusters(number = 20, cnumber = 5, csize = 10):
# Note that the way the clusters are positioned is Gaussian randomness.
rnum = np.random.rand(cnumber, 2)
rn = rnum[:,0] * number
rn = rn.astype(int)
rn[np.where(rn < 5 )] = 5
rn[np.where(rn > number/2. )] = round(number / 2., 0)
ra = rnum[:,1] * 2.9
ra[np.where(ra < 1.5)] = 1.5
cls = np.random.randn(number, 3) * csize
# Random multipliers for central point of cluster
rxyz = np.random.randn(cnumber-1, 3)
for i in xrange(cnumber-1):
tmp = np.random.randn(rn[i+1], 3)
x = tmp[:,0] + ( rxyz[i,0] * csize )
y = tmp[:,1] + ( rxyz[i,1] * csize )
z = tmp[:,2] + ( rxyz[i,2] * csize )
tmp = np.column_stack([x,y,z])
cls = np.vstack([cls,tmp])
return cls
# Generate a cluster of clusters and distance matrix.
cls = clusters()
D = pdist(cls[:,0:2])
D = squareform(D)
# Compute and plot first dendrogram.
fig = mpl.figure(figsize=(8,8))
ax1 = fig.add_axes([0.09,0.1,0.2,0.6])
Y1 = hy.linkage(D, method='complete')
cutoff = 0.3 * np.max(Y1[:, 2])
Z1 = hy.dendrogram(Y1, orientation='right', color_threshold=cutoff)
ax1.xaxis.set_visible(False)
ax1.yaxis.set_visible(False)
# Compute and plot second dendrogram.
ax2 = fig.add_axes([0.3,0.71,0.6,0.2])
Y2 = hy.linkage(D, method='average')
cutoff = 0.3 * np.max(Y2[:, 2])
Z2 = hy.dendrogram(Y2, color_threshold=cutoff)
ax2.xaxis.set_visible(False)
ax2.yaxis.set_visible(False)
# Plot distance matrix.
ax3 = fig.add_axes([0.3,0.1,0.6,0.6])
idx1 = Z1['leaves']
idx2 = Z2['leaves']
D = D[idx1,:]
D = D[:,idx2]
ax3.matshow(D, aspect='auto', origin='lower', cmap=mpl.cm.YlGnBu)
ax3.xaxis.set_visible(False)
ax3.yaxis.set_visible(False)
# Plot colorbar.
fig.savefig('cluster_hy_f01.pdf', bbox = 'tight')
使用一个叫做fcluster的函数区分结构，它在某个阈值上提供了每个集群的索引。fcluster的输出将取决于计算链接函数时使用的方法，如完整或单个。分配给集群的截止值作为fcluster函数中的第二个输入提供。在树状图函数中，截止值的默认值为0.7*np.max(Y[：，2])，但在这里我们将使用与上例相同的截止值，扩展器为0.3。
# Same imports and cluster function from the previous example
# follow through here.
# Here we define a function to collect the coordinates of
# each point of the different clusters.
def group(data, index):
number = np.unique(index)
groups = []
for i in number:
groups.append(data[index == i])
return groups
# Creating a cluster of clusters
cls = clusters()
# Calculating the linkage matrix
Y = hy.linkage(cls[:,0:2], method='complete')
# Here we use the fcluster function to pull out a
# collection of flat clusters from the hierarchical
# data structure. Note that we are using the same
# cutoff value as in the previous example for the dendrogram
# using the 'complete' method.
cutoff = 0.3 * np.max(Y[:, 2])
index = hy.fcluster(Y, cutoff, 'distance')
# Using the group function, we group points into their
# respective clusters.
groups = group(cls, index)
# Plotting clusters
fig = mpl.figure(figsize=(6, 6))
ax = fig.add_subplot(111)
colors = ['r', 'c', 'b', 'g', 'orange', 'k', 'y', 'gray']
for i, g in enumerate(groups):
i = np.mod(i, len(colors))
ax.scatter(g[:,0], g[:,1], c=colors[i], edgecolor='none', s=50)
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
fig.savefig('cluster_hy_f02.pdf', bbox = 'tight')
6.信号和图像处理
scipy允许我们轻松读写JPEG,PNG格式的图片，不用考虑吧图像文件结构。
下面展示了从国际空间站制作图片的代码
import numpy as np
from scipy.misc import imread, imsave
from glob import glob
# Getting the list of files in the directory
files = glob('space/*.JPG')
# Opening up the first image for loop
im1 = imread(files[0]).astype(np.float32)
# Starting loop and continue co-adding new images
for i in xrange(1, len(files)):
print i
im1 += imread(files[i]).astype(np.float32)
# Saving img
imsave('stacked_image.jpg', im1)
python环境中的JPG图像是（426,640,3）的ndarray，其中三层分别是红色、绿色和蓝色。
7.稀疏矩阵Sparse Matrices
有了numpy，我们可以在包含106个元素的阵列上以合理的速度操作。一旦我们达到107个元素，操作就会开始变慢，python的内存将取决于计算机内存变得有限，如果这些大规模数组主要包含零，那么就很幸运了，因为这是稀疏矩阵的属性。如果稀疏矩阵处理正确，操作时间和内存使用量可能会大幅下降。下面的简单例子说明了这一点。
import numpy as np
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigh
import scipy.sparse
import time
N = 3000
# Creating a random sparse matrix
m = scipy.sparse.rand(N, N)
# Creating an array clone of it
a = m.toarray()
print('The numpy array data size: ' + str(a.nbytes) + ' bytes')
print('The sparse matrix data size: ' + str(m.data.nbytes) + ' bytes')
# Non-sparse
t0 = time.time()
dt = str(np.round(time.time() - t0, 3)) + ' seconds'
print('Non-sparse operation takes ' + dt)
# Sparse
t0 = time.time()
res2 = eigsh(m)
dt = str(np.round(time.time() - t0, 3)) + ' seconds'
print('Sparse operation takes ' + dt)
稀疏矩阵的内存效率是它的100倍，特征操作大约是非稀疏情况的快150倍。
在二维和三维几何学中，有许多用于工程学、计算流体动力学、计算流体动力学、电磁学、热力学和声学等领域的稀疏数据结构。稀疏矩阵的非几何实例适用于优化、经济建模、数学和统计学以及网络/图论。使用scipy可以读写常见的稀疏矩阵文件格式，如Matrix Market 和 Harwell-Boeing，或者加载MatLab文件。这对于与使用这些数据格式的其他人的协作特别有用。
7.Numpy之外读写文件
numpy提供了一组良好的输入和输出用于ASCII文件。如果只需要共享信息以从一个python环境到另一个环境读取，它的二进制支持很好。但是，那么更普遍使用的二进制文件格式呢？如果您正在使用Matlab或与正在使用它的其他人合作，读写支持Matlab的文件(使用scipy.io.loadmat和scipy.savemat)不是问题。
在天文学、地理和医学等领域，有一种叫做IDL的编程语言。它以二进制格式保存文件，并可以使用一个名为scipy.io.readsav的内置包读取文件。它是一个灵活和快速的模块，但它没有编写能力。
最后，您可以查询、读取和写入Matrix Market文件。这些文件非常常用于共享以智能格式编写的矩阵数据结构。这种格式在其他语言中，如C、格式格式和矩阵语言，都很受支持，所以由于它的普遍性和用户可读性，它是一个很好的使用格式。它也适用于稀疏矩阵。
参考书籍：scipy_and_numpy
扫一扫
718
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
