一、Introduction

1. 安装Numpy 与SciPy

（1）使用Python发行包

常见的Python发行包：

Enthought Python Distribution (EPD)

ActivePython (AP)

Python(x,y)

（2）Python Windows

Step1：安装Python

下载https://http://www.sodocs.net/doc/c60f07d3cc7931b764ce153d.html/downloads/release/python-2710/ ，64位操作系统选择Windows x86_64 MSI installer，32位操作系统选择Windows x86 MSI installer

双击安装（最好选择默认路径）

Step 2：安装包管理器pip

a. 下载https://http://www.sodocs.net/doc/c60f07d3cc7931b764ce153d.html/pypi/pip/ ，选择pip-8.1.1.tar.gz，解压。大命令行进入dist\pip-8.1.1

b. 运行python setup.py install

c. 修改环境变量PATH，添加C:\Python27\Scripts

Step 3：使用包管理器安装numpy、scipy

NumPy和SciPy在windows下需要手动安装,否则容易出现意外的错误。过程如下

a. 从http://www.sodocs.net/doc/c60f07d3cc7931b764ce153d.html/~gohlke/pythonlibs/ 下载相应whl包（下载python2.7版本cp27，32位下载win32，64位下载win_amd64）手动安装，手动安装注意依赖关系

b. pip install wheel

c. pip install xxxxx.whl

Step 4：安装其他包

pip install pillow

pip install pandas

pip install -U scikit-learn

pip install matplotlib

pip install ipython

pip install pyreadline

（3）Python OS X

Step 1：安装homebrew

/usr/bin/ruby -e "$(curl -fsSL

https://http://www.sodocs.net/doc/c60f07d3cc7931b764ce153d.html/Homebrew/install/master/install)"

Step 2：重新安装python

brew install python

Step 3：安装pip

sudo easy_install pip

Step 4：安装其他包

pip install numpy

pip install scipy

pip install pillow

pip install pandas

pip install -U scikit-learn

pip install matplotlib

pip install ipython

二、Numpy

1. numpy数组与python列表效率对比

import numpy as np

# 创建大小为10^7 的数组

arr = np.arange(1e7)

larr = arr.tolist()

def list_times(alist, scalar):

for i, val in enumerate(alist):

alist[i] = val * scalar

return alist

# 利用IPython的魔术方法timeit计算运行时间

timeit arr * 1.1

>>> 1 loops, best of 3: 76.9 ms per loop

timeit list_times(larr, 1.1)

>>> 1 loops, best of 3: 2.03 s per loop

2. 创建数组并设置数据类型

（1）从列表转换

alist = [1, 2, 3]

arr = np.array(alist)

（2）np.arange()

arr = np.arange(100)

arr = np.arange(10,100)

（3）np.zeros()

arr = np.zeros(5)

np.zeros((5,5))

cube = np.zeros((5,5,5)).astype(int) + 1

cube = np.ones((5, 5, 5)).astype(np.float16)

arr = np.zeros(2, dtype=int)

arr = np.zeros(2, dtype=np.float32)

（4）reshape()

arr1d = np.arange(1000)

arr3d = arr1d.reshape((10,10,10))

arr3d = np.reshape(arr1s, (10, 10, 10))

（5）revel()

作用与reshape相反

（6）shape

显示数据对象的形状

arr1d.shape

注意：对数据形状结构的改变只是改变了数据的显示形式，即只是改变了数据的引用，对一个数据的改变另一个也会被改变。

3. 记录数组

（1）创建记录数组并赋值

recarr = np.zeros((2,), dtype=('i4,f4,a10'))

#创建大小为2的记录数组，类型为4字节整数、4字节浮点数和10字节字符

recarr[:] = [(1,2.,'Hello'),(2,3.,"World")]

（2）使用zip()

recarr = np.zeros((2,), dtype=('i4,f4,a10'))

col1 = np.arange(2) + 1

col2 = np.arange(2, dtype=np.float32)

col3 = ['Hello', 'World']

recarr[:] = zip(col1, col2, col3)

（3）为每列数据命名

http://www.sodocs.net/doc/c60f07d3cc7931b764ce153d.htmls = ('Integers' , 'Floats', 'Strings')

（4）使用列名访问数据

recarr('Integers')

4. 索引和切片

（1）numpy提供了类似于matlab的索引和切片

alist=[[1,2],[3,4]]

alist[0][1] #python方式

arr = np.array(alist)

arr[0,1] #单个元素

arr[:,1] #第1列

arr[1,:] #第1行

（2）np.where()

根据条件获取索引号

index = np.where(arr>2)

new_arr = arr[index]

new_arr = np.delete(arr, index)

也可以这样操作：

index = arr > 2 #得到一个逻辑数组

new_arr = arr[index]

注意：第二种方法速度更快，而且可以用“～ index”很容易的得到与index相反的逻辑数组。

5. NumPy数组的布尔操作

NumPy数组元素可以通过逻辑表达式方便的操作

例：

# 创建如Plot A所示的数组

img1 = np.zeros((20, 20)) + 3

img1[4:-4, 4:-4] = 6

img1[7:-7, 7:-7] = 9

# 获取数值大于2且小于6的元素索引

index1 = img1 > 2

index2 = img1 < 6

compound_index = index1 & index2

# 上式与下式结果相同

compound_index = (img1 > 3) & (img1 < 7)

img2 = np.copy(img1)

img2[compound_index] = 0

# 得到Plot B.

# 更复杂的数组逻辑操作

index3 = img1 == 9

index4 = (index1 & index2) | index3

img3 = np.copy(img1)

img3[index4] = 0 # 得到Plot C.

例：

import numpy.random as rand

a = rand.randn(100)

index = a > 0.2

b = a[index]

b = b ** 2 – 2

a[index] = b

6. 读写操作

（1）Python读写文本文件

f = open('somefile.txt', 'r') #以只读方式打开文件，'r'表示读

alist = f.readlines() #将文件内容读入列表，每一行为一个列表元素file f.close() #关闭文件

f = open('newtextfile.txt', 'w') #以可写方式打开文件，'w'表示写

f.writelines(newdata) #写入数据

f.close() #关闭文件

注意：读写完毕之后要将文件关闭

（2）Numpy文件文件读写

Python读写文件文件虽然方便且效率很好，但是不太适合处理极大的文件。当文件内容有结构，且为数字时用NumPy处理，存numpy.ndarray会更合适。

例：

import numpy as np

arr = np.loadtxt('somefile.txt')

np.savetxt('somenewfile.txt')

如果文件各列数据类型不一样，则需要指明数据类型，NumPy用来保存数据的类型为recarray，可以用处理ndarray同样的方法来对元素进行操作。recarray数据类型不能直接保存为文本文件，如果需要的话可以使用matplotlib.mlab实现。

例：

文件example.txt内容如下

XR21 32.789 1

XR22 33.091 2

读入数据

table = np.loadtxt('example.txt',

dtype='names': ('ID', 'Result', 'Type'),

'formats': ('S4', 'f4', 'i2'))

提示：如果文本数据为ASCII格式的，使用Asciitable包读写会更加高效。

（3）二进制文件

文本文件处理简单方便，但是读写速度和文件大小都不能和二进制文件相比，因此大数据处理适合使用二进制文件。

例：

import numpy as np

data = np.empty((1000, 1000)) #创建一个较大的数组

np.save('test.npy', data) #保存数据

np.savez('test.npz', data) #压缩保存数据

newdata = np.load('test.npy') #读入数据

注意：NumPy使用numpy.save和numpy.load来读写二进制文件，但这种二进制文件只能在NumPy下读写，scipy.io可以处理更通用的二进制文件

7. 数学运算

（1）线性代数

NumPy数组间的运算只是相对应元素间的远算，不能用运算符进行矩阵运算，可以使用numpy.dot 和numpy.transpose分别来进行矩阵乘法运算和矩阵转置。其优点在于常规操作时避免了对数据遍历。

NumPy的matrix类型则可以直接用运算符号进行运算。

例：使用matrix解方程组

import numpy as np

A = np.matrix([[3, 6, -5], [1, -3, 2], [5, -1, 4]]) #定义矩阵

B = np.matrix([[12], [-2], [10]])

X = A ** (-1) * B #求方程组

print(X)

例：使用数组解方程组

import numpy as np

a = np.array([[3, 6, -5], [1, -3, 2], [5, -1, 4]])

b = np.array([12, -2, 10])

x = np.linalg.inv(a).dot(b)

print(x)

注意：数组的运算速度更快，而且为了在使用中保持数据类型一致，建议使用数组。

三、SciPy

1. 最优化

（1）数据建模和拟合

SciPy函数curve_fit使用基于卡方的方法进行线性回归分析。下面，首先使用f(x)=ax+b生成带有噪声的数据，然后使用用curv_fit来拟合。

例：线性回归

import numpy as np

from scipy.optimize import curve_fit

# 创建函数f(x)=ax+b

def func(x, a, b):

return a * x + b

# 创建干静数据

x = np.linspace(0, 10, 100)

y = func(x, 1, 2)

# 添加噪声

yn = y + 0.9 * np.random.normal(size=len(x))

# 拟合噪声数据

popt, pcov = curve_fit(func, x, yn)

# 输出最优参数

print(popt)

例：高斯分布拟合

# 创建函数

def func(x, a, b, c):

return a*np.exp(-(x-b)**2/(2*c**2))

# 生成干静数据

x = np.linspace(0, 10, 100)

y = func(x, 1, 5, 2)

# 添加噪声

yn = y + 0.2 * np.random.normal(size=len(x))

# 拟合

popt, pcov = curve_fit(func, x, yn)

（2）函数求解

SciPy的optimize模块中有大量的函数求解工具，fsolve是其中最常用的。例：线性函数求解

from scipy.optimize import fsolve

import numpy as np

line = lambda x: x + 3

solution = fsolve(line, -2)

print solution

例：求函数交叉点

from scipy.optimize import fsolve

import numpy as np

# 用于求解的解函

def findIntersection(func1, func2, x0):

return fsolve(lambda x : func1(x) - func2(x), x0) # 两个函数

funky = lambda x : np.cos(x / 5) * np.sin(x / 2)

line = lambda x : 0.01 * x - 0.5

x = np.linspace(0,45,10000)

result = findIntersection(funky, line, [15, 20, 30, 35, 40, 45]) # 输出结果

print(result, line(result))

2. 插值

（1）interp1d

例：正弦函数插值

import numpy as np

from scipy.interpolate import interp1d

# 创建待插值的数据

x = np.linspace(0, 10 * np.pi, 20) y = np.cos(x)

# 分别用linear和quadratic插值

fl = interp1d(x, y, kind='linear')

fq = interp1d(x, y, kind='quadratic')

xint = np.linspace(x.min(), x.max(), 1000)

yintl = fl(xint)

yintq = fq(xint)

（2）UnivariateSpline

例：噪声数据插值

import numpy as np

from scipy.interpolate import UnivariateSpline

# 创建含噪声的待插值数据

sample = 30

x = np.linspace(1, 10 * np.pi, sample)

y = np.cos(x) + np.log10(x) + np.random.randn(sample) / 10

# 插值，参数s为smoothing factor

f = UnivariateSpline(x, y, s=1)

xint = np.linspace(x.min(), x.max(), 1000)

yint = f(xint)

（3）griddata

例：利用插值重构图片

import numpy as np

from scipy.interpolate import griddata

# 定义一个函数

ripple = lambda x, y: np.sqrt(x**2 + y**2)+np.sin(x**2 + y**2)

# 生成grid数据，复数定义了生成grid数据的step，若无该复数则step为5 grid_x, grid_y = np.mgrid[0:5:1000j, 0:5:1000j]

# 生成待插值的样本数据

xy = np.random.rand(1000, 2)

sample = ripple(xy[:,0] * 5 , xy[:,1] * 5)

# 用cubic方法插值

grid_z0 = griddata(xy * 5, sample, (grid_x, grid_y), method='cubic')

上图中，左侧为原始数据，其中的黑点是待插值的样本，右图为插值后的数据。要想提高质量，生成更大的样本数据即可。

（4）SmoothBivariateSpline

例：利用插值重构图片

import numpy as np

from scipy.interpolate import SmoothBivariateSpline as SBS

# 定义函数

ripple = lambda x, y: np.sqrt(x**2 + y**2)+np.sin(x**2 + y**2)

# 生成待插值样本

xy= np.random.rand(1000, 2)

x, y = xy[:,0], xy[:,1]

sample = ripple(xy[:,0] * 5 , xy[:,1] * 5)

# 插值

fit = SBS(x * 5, y * 5, sample, s=0.01, kx=4, ky=4)

interp = fit(np.linspace(0, 5, 1000), np.linspace(0, 5, 1000))

注意：SmoothBivariateSpline有时候表现比Spline更好一些，但是它对样本数据更敏感一些，相对而言Spline更加健壮。

3. 积分

SicPy中的积分是近似的数值积分，SymPy是一个符号积分的工具包。

（1）解析积分

例：

import numpy as np

from scipy.integrate import quad

# 定义被积函数

func = lambda x: np.cos(np.exp(x)) ** 2

solution = quad(func, 0, 3)

print solution

# 输出结果中第一个数字为积分值，第二个为误差

# (1.296467785724373, 1.397797186265988e-09)

（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.

4. 统计

SciPy中有包括mean, std, median, argmax, 及argmin等在内的基本统计函数，而且numpy.arrays类型中内置了大部分统计函数，以便于使用。

import numpy as np

# 创建大小为1000的随机数组

elements x = np.random.randn(1000)

mean = x.mean() #均值

std = x.std() #标准差

var = x.var() #方差

SciPy中还包括了各种分布、函数等工具。

（1）连续和离散分布

SciPy的scipy.stats包中包含了大概80种连续分布和10种离散分布。下图是其中的20种连续分布的概率密度函数。这些分布函数其实都依赖于numpy.random函数。

有几种方法来使用scipy.stats中的分布时：概率密度函数（PDFs）、累积分布函数（CDFs）、随机变量样本（RVSs）、百分比点函数（PPFs）等。下面基于标准正态分布函数，来演示如何使用这些分布。

例：

import numpy as np

import scipy.stats import norm

# 创建样本区间

x = np.linspace(-5,5,1000)

# 设置正态分布参数，loc为均值，scale为标准差

dist = norm(loc=0, scale=1)

# 得到正态分布的PDF 和CDF

pdf = dist.pdf(x)

cdf = dist.cdf(x)

# 根据分布生成500个随机数

sample = dist.rvs(500)

可以基于SciPy.stats中的任何连续分布生成随机数，有需要请查看文档。除此外，如泊松分布、二项分布、几何分布等离散分布的使用也很简单。下式为几何分布的概率分布函数（PMF）：

例：

import numpy as np

from scipy.stats import geom

# 设置几何分布的参数

p = 0.5

dist = geom(p)

# 设置样本区间

x = np.linspace(0, 5, 1000)

# 得到几何分布的PMF 和CDF

pmf = dist.pmf(x)

cdf = dist.cdf(x)

# 生成500个随机数

sample = dist.rvs(500)

（2）函数

SciPy中有超过60种统计函数。stats包中包括了诸如kstest 和normaltest等的样本测试函数，用来检测样本是否服从某种分布。提示：在使用这些工具前，要对数据有较好的理解，否则可能会误读它们的结果。

例：样本分布检验

import numpy as np

from scipy import stats

# 生成包括100个服从正态分布的随机数样本

sample = np.random.randn(100)

# 用normaltest检验原假设

out = stats.normaltest(sample)

print('normaltest output')

print('Z-score = ' + str(out[0]))

print('P-value = ' + str(out[1]))

# kstest 是检验拟合度的Kolmogorov-Smirnov检验，这里针对正态分布进行检验，# D是KS统计量的值，越接近0越好

out = stats.kstest(sample, 'norm')

print('\nkstest output for the Normal distribution')

print('D = ' + str(out[0]))

print('P-value = ' + str(out[1]))

# 类似地可以针对其他分布进行检验，例如Wald分布

out = stats.kstest(sample, 'wald')

print('\nkstest output for the Wald distribution')

print('D = ' + str(out[0]))

print('P-value = ' + str(out[1]))

SciPy的stats模块中还提供了一些描述函数，如几何平均（gmean）、偏度（skew）、样本频数（itemfreq）等。

例：

import numpy as np

from scipy import stats

#生成包括100个服从正态分布的随机数样本

sample = np.random.randn(100)

# 调和平均数，样本值须大于0

out = stats.hmean(sample[sample > 0])

print('Harmonic mean = ' + str(out))

# 计算-1到1之间样本的均值

out = stats.tmean(sample, limits=(-1, 1))

print('\nTrimmed mean = ' + str(out))

# 计算样本偏度

out = stats.skew(sample)

print('\nSkewness = ' + str(out))

# 函数describe可以一次给出样本的多种描述统计结果

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

SciPy的stats模块中还有很多统计工具，可以满足绝大多数需要。还可以用RPy，通过它能够在Python中调用R语言进行统计分析。此外，Pandas是python的一个强大的工具包，它可以在大数据上进行快速的统计分析。

5. 空间和聚类分析

SciPy包括scipy.spatial类和scipy.cluster类分别用于空间和聚类分析。前者用于分析数据点之间的距离，后者包括两个子类矢量量化（vq）和层次聚类（hierarchy）。

（1）矢量量化（Vector Quantization）

矢量量化是信号处理、数据压缩和聚类等领域通用的术语。这里仅关注其在聚类中的应用。

例：k均值聚类

import numpy as np

from scipy.cluster import vq

# 生成数据

c1 = np.random.randn(100, 2) + 5

c2 = np.random.randn(30, 2) - 5

c3 = np.random.randn(50, 2)

# 将所有数据放入一个180 x 2 的数组

data = np.vstack([c1, c2, c3])

# 利用k均值方法计算聚类的质心和方差

centroids, variance = vq.kmeans(data, 3)

# 变量identified中存放关于数据聚类的信息

identified, distance = vq.vq(data, centroids)

# 获得各类别的数据

vqc1 = data[identified == 0]

vqc2 = data[identified == 1]

vqc3 = data[identified == 2]

（2）层次聚类

层次聚类是一种重要的聚类方法，但其输出结果比较复杂，不能像k均值那样给出清晰的聚类结果。下面是一个层次聚类的例子，输入一个距离矩阵，输出为一个树状图。例：

# coding:utf-8

import numpy as np

import matplotlib.pyplot as mpl

from scipy.spatial.distance import pdist, squareform

import scipy.cluster.hierarchy as hy

# 用于生成聚类数据的函数

def clusters(number=20, cnumber=5, csize=10):

# 聚类服从高斯分布

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

# 创建待聚类数据及矩离矩阵

cls = clusters()

D = pdist(cls[:, 0:2])

D = squareform(D)

# 绘制左侧树状图

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='left', color_threshold=cutoff) ax1.xaxis.set_visible(False)

ax1.yaxis.set_visible(False)

#绘制顶部树状图

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)

# 显示距离矩阵

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)

# 保存图片，显示图片

fig.savefig('cluster_hy_f01.pdf', bbox = 'tight')

mpl.show()

在上图虽然计算了数据点之间的距离，但是还是难以将各类区分开。函数fcluster 可以根据阈值来区分各类，其输出结果依赖于linkage函数所采用的方法，如complete 或single等，它的第二个参数即是阈值。dendrogram函数中默认的阈值是0.7 * np.max(Y[:, 2])，这里还使用0.3。

例：

# 导入的包同上例一致, 函数cluster同上例

# 获得不同类别数据点的坐标

def group(data, index):

number = np.unique(index)

groups = []

for i in number:

groups.append(data[index == i])

return groups

# 创建数据

cls = clusters()

# 计算kinkage矩阵

Y = hy.linkage(cls[:,0:2], method='complete')

# 从层次数据结构中, 用fcluster函数将层次结构的数据转为flat clusters

cutoff = 0.3 * np.max(Y[:, 2])

index = hy.fcluster(Y, cutoff, 'distance')

# 使用grooup函数将数据划分类别

groups = group(cls, index)

# 绘制数据点

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

mpl.show()

6. 信号和图像处理

SimPy可以高效的处理JPEG和PNG等格式的图片，下面利用SimPy将多副图片叠加成一副图片。

本部分实例更多内容参见：

http://www.sodocs.net/doc/c60f07d3cc7931b764ce153d.html/questions/9251580/stacking-astronomy-images-with-python 例：

import numpy as np

from scipy.misc import imread, imsave

from glob import glob

# 获取文件夹中全部JPG图片

files = glob('space/*.JPG')

# 打开第一个图片

im1 = imread(files[0]).astype(np.float32)

# 叠加图片

for i in xrange(1, len(files)):

print i

im1 += imread(files[i]).astype(np.float32)

# 保存图片

imsave('stacked_image.jpg', im1)

例：

import numpy as np

from scipy.misc import imread, imsave

from glob import glob

# 该函数对比两副图片,选择两副图片中相同位置上较亮的点生成新的图片def chop_lighter(image1, image2):

s1 = np.sum(image1, axis=2)

s2 = np.sum(image2, axis=2)

index = s1 < s2

image1[index, 0] = image2[index, 0]

image1[index, 1] = image2[index, 1]

image1[index, 2] = image2[index, 2]

return image1

files = glob('space/*.JPG')

im1 = imread(files[0]).astype(np.float32)

im2 = np.copy(im1)

for i in xrange(1, len(files)):

print i