Python最好用的科学计算库：NumPy快速入门教程（三）
由iquant创建，最终由iquant 被浏览 27 用户
该文章接上一篇Python最好用的科学计算库：NumPy快速入门教程（二）
import numpy as np
%matplotlib inline
深入理解 NumPy
广播
广播操作允许通用函数能够对不是严格相同形状的输入进行处理。
广播的第一个规则是，如果所有输入数组不具有相同数量的维度，则将“1”重复地预先添加到形状较小的数组中，直到所有数组具有相同数量的维度。
广播的第二个规则是，确保有某个维度大小为1的数组，操作起来像该维度与在该维度上具有最大形状的数组一样。假定数组的元素沿着“广播”的那个维度是相同的。
应用了广播规则后，所有数组的维度应该相互匹配。更多请查看Broadcasting
高级索引和技巧
NumPy比Python的内置类型提供了更多的索引工具。除了我们之前提到过了切片和通过整数索引以外，还可以通过整数数组和布尔数组进行索引。
通过数组索引
>>> a = np.arange(12)**2 # 创建1-12个数的平方值
>>> i = np.array( [ 1,1,3,8,5 ] ) # 创建一个用于索引的数组
>>> a[i] # 数组a中位置i的值
array([ 1, 1, 9, 64, 25])
>>> j = np.array( [ [ 3, 4], [ 9, 7 ] ] ) # 创建一个用于索引的二维数组
>>> a[j] # 通过j索引得到与j相同维度的数组
array([[ 9, 16],
[81, 49]])
当被索引的数组a是多维度的，单一的索引数组只对a的第一个维度操作。以下的例子展示的这个特性，该例子将使用一个‘调色板’将一个标签的图片转换为彩色图片。
>>> palette = np.array( [ [0,0,0], # 黑色
... [255,0,0], # 红色
... [0,255,0], # 绿色
... [0,0,255], # 蓝色
... [255,255,255] ] ) # 白色
>>> image = np.array( [ [ 0, 1, 2, 0 ], # 创建一个用于索引palette数组的数组
... [ 0, 3, 4, 0 ] ] )
>>> palette[image] # 生成的图片形状是(2，4，3),因为索引数组的形状是(2，4),被索引数组的
... # 第一个维度中每个元素都是一个形状为(3)的数组
array([[[ 0, 0, 0],
[255, 0, 0],
[ 0, 255, 0],
[ 0, 0, 0]],
[[ 0, 0, 0],
[ 0, 0, 255],
[255, 255, 255],
[ 0, 0, 0]]])
我们也可以使用多维索引数组。多维索引数组的维度必须与被索引数组的维度一致。
>>> a = np.arange(12).reshape(3,4)
>>> a
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> i = np.array( [ [0,1], # 用于索引a的第一维度的数组
... [1,2] ] )
>>> j = np.array( [ [2,1], # 用于索引a的第二维度的数组
... [3,3] ] )
>>> b = a[i,j] # i和j必须具有相同的维度, b[1,0]相当于a[i[1,0],j[1,0]]
>>> b
array([[ 2, 5],
[ 7, 11]])
>>> a[i,2]
array([[ 2, 6],
[ 6, 10]])
>>> a[:,j] # i.e., a[ : , j]
array([[[ 2, 1],
[ 3, 3]],
[[ 6, 5],
[ 7, 7]],
[[10, 9],
[11, 11]]])
我们自然的想到，可以将i和j放在一个python内置类型中（比如list），然后使用list进行索引。
>>> l = [i,j]
>>> a[l] # 相当于a[i,j]
array([[ 2, 5],
[ 7, 11]])
但是，我们不能将i和j组成一个数组，因为用一个数组索引将被理解为对a的第一个维度进行索引。
>>> s = np.array( [i,j] )
>>> a[s] # not what we want
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-30-8891efffcfb5> in <module>()
1 s = np.array( [i,j] )
----> 2 a[s] # not what we want
IndexError: index 3 is out of bounds for axis 0 with size 3
>>> a[tuple(s)] # 相当于a[i,j]
array([[ 2, 5],
[ 7, 11]])
另一个使用索引数组的应用是搜索一个时间序列的最大值。
>>> time = np.linspace(20, 145, 5) # time scale
>>> data = np.sin(np.arange(20)).reshape(5,4) # 4 time-dependent series
>>> time
array([ 20. , 51.25, 82.5 , 113.75, 145. ])
>>> data
array([[ 0. , 0.84147098, 0.90929743, 0.14112001],
[-0.7568025 , -0.95892427, -0.2794155 , 0.6569866 ],
[ 0.98935825, 0.41211849, -0.54402111, -0.99999021],
[-0.53657292, 0.42016704, 0.99060736, 0.65028784],
[-0.28790332, -0.96139749, -0.75098725, 0.14987721]])
>>> ind = data.argmax(axis=0) # 每个序列最大值的索引
>>> ind
array([2, 0, 3, 1])
>>> time_max = time[ind] # 最大值对应的时间
>>> data_max = data[ind, range(data.shape[1])] # 相当于 data[ind[0],0], data[ind[1],1]...
>>> time_max
array([ 82.5 , 20. , 113.75, 51.25])
>>> data_max
array([0.98935825, 0.84147098, 0.99060736, 0.6569866 ])
>>> np.all(data_max == data.max(axis=0)) # 以上操作相当于data.max(axis = 0)
True
你也可以使用索引序列指定赋值的目标：
>>> a = np.arange(5)
>>> a
array([0, 1, 2, 3, 4])
>>> a[[1,3,4]] = 0
>>> a
array([0, 0, 2, 0, 0])
但是，当索引包含重复值，赋值操作被执行多次，并保留最后一次的值。
>>> a = np.arange(5)
>>> a[[0,0,2]]=[1,2,3]
>>> a
array([2, 1, 3, 3, 4])
这很合理，但是，当你想使用Python的+=操作时要小心，因为它可能得到的不是你想要的结果。
>>> a = np.arange(5)
>>> a[[0,0,2]]+=1
>>> a
array([1, 1, 3, 3, 4])
尽管0在列表中出现了两次，但是第0个元素只增加了1次。因为，在python中，“a+=1”等于“a = a + 1”。
用布尔数组做索引
当我们使用整数数组索引时，我们提供了我们想要选择的值的索引。但是使用布尔数组时，方式就变了；我们明确地选择哪个值是我们想要的，哪个不是。
最自然的使用方式就是用一个与原数组形状相同的布尔数组作为索引。
>>> a = np.arange(12).reshape(3,4)
>>> b = a > 4
>>> b #我们得到了一个与a形状一样的布尔数组
array([[False, False, False, False],
[False, True, True, True],
[ True, True, True, True]])
>>> a[b] # 使用b作为索引数组将返回一个一维数组
array([ 5, 6, 7, 8, 9, 10, 11])
这个性质在赋值时非常有用：
>>> a[b] = 0 # 所有大于4的值被赋值为0
>>> a
array([[0, 1, 2, 3],
[4, 0, 0, 0],
[0, 0, 0, 0]])
下面这个例子关于如何使用布尔数组生成Mandelbrot set的图片。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> def mandelbrot( h,w, maxit=20 ):
... """Returns an image of the Mandelbrot fractal of size (h,w)."""
... y,x = np.ogrid[ -1.4:1.4:h*1j, -2:0.8:w*1j ]
... c = x+y*1j
... z = c
... divtime = maxit + np.zeros(z.shape, dtype=int)
...
... for i in range(maxit):
... z = z**2 + c
... diverge = z*np.conj(z) > 2**2 # who is diverging
... div_now = diverge & (divtime==maxit) # who is diverging now
... divtime[div_now] = i # note when
... z[diverge] = 2 # avoid diverging too much
...
... return divtime
>>> plt.imshow(mandelbrot(400,400))
>>> plt.show();
使用布尔数组进行索引的第二种方式与使用整数索引类似；对于数组的每个维度，我们使用一个一维的布尔数组选择我们想要的切片。
>>> a = np.arange(12).reshape(3,4)
>>> b1 = np.array([False,True,True]) # 用于选择第一个维度的值
>>> b2 = np.array([True,False,True,False]) # 用于选择第二个维度的值
>>> a[b1,:] #选择行
array([[ 4, 5, 6, 7],
[ 8, 9, 10, 11]])
>>> a[:,b2] #选择列
array([[ 0, 2],
[ 4, 6],
[ 8, 10]])
>>> a[b1,b2] #很意外的结果
array([ 4, 10])
需要注意的是，你的一维布尔索引的大小必须与你要索引的原数组的对应维度相同。
ix_()函数
ix_函数可用于组合不同的向量，以便获得对每个向量元素进行操作的结果。 例如，如果要计算从每个向量a，b和c中取得的所有a + b * c的值：
>>> a = np.array([2,3,4,5])
>>> b = np.array([8,5,4])
>>> c = np.array([5,4,6,8,3])
>>> ax,bx,cx = np.ix_(a,b,c)
>>> ax
array([[[2]],
[[3]],
[[4]],
[[5]]])
>>> bx
array([[[8],
[5],
[4]]])
>>> cx
array([[[5, 4, 6, 8, 3]]])
>>> ax.shape, bx.shape, cx.shape
((4, 1, 1), (1, 3, 1), (1, 1, 5))
>>> result = ax+bx*cx
>>> result
array([[[42, 34, 50, 66, 26],
[27, 22, 32, 42, 17],
[22, 18, 26, 34, 14]],
[[43, 35, 51, 67, 27],
[28, 23, 33, 43, 18],
[23, 19, 27, 35, 15]],
[[44, 36, 52, 68, 28],
[29, 24, 34, 44, 19],
[24, 20, 28, 36, 16]],
[[45, 37, 53, 69, 29],
[30, 25, 35, 45, 20],
[25, 21, 29, 37, 17]]])
>>> result[3,2,4]
17
>>> a[3]+b[2]*c[4] #对于三个不同维度的数组，我们大大简化了使用三个嵌套的循环来计算所有对应值的a+b*c的操作
17
你也可是使用下面的方式部署：
>>> def ufunc_reduce(ufct, *vectors):
... vs = np.ix_(*vectors)
... r = ufct.identity
... for v in vs:
... r = ufct(r,v)
... return r
然后这样使用
>>> ufunc_reduce(np.add,a,b,c)
array([[[15, 14, 16, 18, 13],
[12, 11, 13, 15, 10],
[11, 10, 12, 14, 9]],
[[16, 15, 17, 19, 14],
[13, 12, 14, 16, 11],
[12, 11, 13, 15, 10]],
[[17, 16, 18, 20, 15],
[14, 13, 15, 17, 12],
[13, 12, 14, 16, 11]],
[[18, 17, 19, 21, 16],
[15, 14, 16, 18, 13],
[14, 13, 15, 17, 12]]])
这个版本的reduce的优点是它利用了广播规则，以避免创建一个参数数组，输出的大小乘以向量的数量。
通过字符串索引
见https://docs.scipy.org/doc/numpy/user/basics.rec.html#structured-arrays
线性代数
继续介绍基本的线性代数运算。
简单数组运算
在numpy文件夹的linalg.py查看更多内容。
>>> import numpy as np
>>> a = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> print(a)
[[1. 2.]
[3. 4.]]
>>> a.transpose() # 数组的转置
array([[1., 3.],
[2., 4.]])
>>> np.linalg.inv(a) #矩阵求逆 ，linalg=linear+algebra
array([[-2. , 1. ],
[ 1.5, -0.5]])
>>> u = np.eye(2) # 创建一个2x2的单位对角矩阵
>>> u
array([[1., 0.],
[0., 1.]])
>>> j = np.array([[0.0, -1.0], [1.0, 0.0]])
>>> np.dot (j, j) # 矩阵内积
array([[-1., 0.],
[ 0., -1.]])
>>> np.trace(u) # 矩阵的迹
2.0
>>> y = np.array([[5.], [7.]])
>>> np.linalg.solve(a, y) #求解由a和y组成的方程组
array([[-3.],
[ 4.]])
>>> np.linalg.eig(j) #计算矩阵的特征向量
(array([0.+1.j, 0.-1.j]),
array([[0.70710678+0.j , 0.70710678-0.j ],
[0. -0.70710678j, 0. +0.70710678j]]))
一些技巧
这里我们列出了一些比较实用的技巧。
自动重塑形状
改变一个数组的形状，你可以忽略其中一个维度，那个维度将自动推断应该是多少：
>>> a = np.arange(30)
>>> a.shape = 2,-1,3 # -1 的意思是“它应该是的值”
>>> a.shape
(2, 5, 3)
>>> a
array([[[ 0, 1, 2],
[ 3, 4, 5],
[ 6, 7, 8],
[ 9, 10, 11],
[12, 13, 14]],
[[15, 16, 17],
[18, 19, 20],
[21, 22, 23],
[24, 25, 26],
[27, 28, 29]]])
向量堆叠
我们如何从一系列同样大小的行向量创建一个二维数组？在matlab中，很简单，如果x和y是同样长度的两个向量，那么我们需要做的就是m=[x;y]。在NumPy中，我们使用column_stack, dstack, hstack 和vstack来实现，取决于你想要在哪个维度上堆叠你的向量。例如：
x = np.arange(0,10,2) # x=([0,2,4,6,8])
y = np.arange(5) # y=([0,1,2,3,4])
m = np.vstack([x,y]) # m=([[0,2,4,6,8],
# [0,1,2,3,4]])
xy = np.hstack([x,y]) # xy =([0,2,4,6,8,0,1,2,3,4])
直方图
NumPy的histogram方法返回两个向量：数组的直方图数据和间隔向量。需要注意的是matplotlib也有用来创建直方图的函数（hist）。两者的不同是，pylab.hist直接将直方图画出来，而numpy.histogram只生成数据。
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # Build a vector of 10000 normal deviates with variance 0.5^2 and mean 2
>>> mu, sigma = 2, 0.5
>>> v = np.random.normal(mu,sigma,10000)
>>> # Plot a normalized histogram with 50 bins
>>> plt.hist(v, bins=50, density=1); # matplotlib version (plot)
>>> plt.show()
>>> #这个是使用NumPy的histogram函数
>>> (n, bins) = np.histogram(v, bins=50, normed=True) # NumPy version (no plot)
>>> plt.plot(.5*(bins[1:]+bins[:-1]), n)
>>> plt.show()
