Python 数据科学 (1)——NumPy(2.Arrays 操作篇)

本文主要讲述 NumPy 中 Arrays 的操作，包括基础属性操作（索引、切片、连接等）；基本计算（计算和广播概念）；聚合计算（最大值、均值等）。这和 Python 中基础数组的操作差不多，但还是有必要进行了解学习，因为很多 Pandas 库里工具都是围绕 NumPy 进行构建的。

Arrays 基础操作

本节主要介绍 Arrays 拥有的特性。为了论述与读者理解的方便性，我们将用随机函数给出具体的两个 Arrays 例子：

In [1]: import numpy as np
In [2]: np.random.seed(0)
In [3]: x1 = np.random.randint(10,size=7)
In [4]: x2 = np.random.randint(10,size=(3,4,5))

Arrays 的基础特性

In [7]: x2.ndim # Arrays的维度数量
Out[7]: 3

In [8]: x2.shape # Arrays每个维度大小
Out[8]: (3, 4, 5)

In [9]: x2.size # Arrays的总大小
Out[9]: 60

In [10]: x2.dtype # Arrays的数据类型
Out[10]: dtype('int64')

In [11]: x2.itemsize # Arrays的单个元素大小(在此例子中我们可以看出'int64'的大小是8 byte)
Out[11]: 8

In [12]: x2.nbytes # Arrays的总字节大小(x2共有60个int64元素，60*8=480 byte)
Out[12]: 480

Arrays 获取元素

与获取 python 原生数组中的元素一致，可以利用下标直接进行获取

In [18]: x1[::2] # x1中所有下标为偶数的元素
Out[18]: array([5, 3, 7, 3])

In [19]: x1[1::2] # x1中所有下标为奇数的元素
Out[19]: array([0, 3, 9]) 

In [20]: x1[::-1] # 反转x1中所有元素
Out[20]: array([3, 9, 7, 3, 3, 0, 5])

In [22]: x2[:2, :3 ,:4] # 获取x2中开始的两个面上的前三行*前四列的数据。
Out[22]:
array([[[5, 2, 4, 7],
       [8, 8, 1, 6],
       [7, 8, 1, 5]], 
      [[3, 5, 0, 2],
       [8, 1, 3, 3],
       [7, 0, 1, 9]]])

以上的获取方法都是引用，如果原 Arrays 中的值改变，那么对应的获取的值将在改变前后的两次中不同。如需使用复制值而不是引用值，则需使用方法 copy():

In [23]: x3 = x1[:3].copy()
In [24]: x3
Out[24]: array([5, 0, 3])

Arrays 的重构(变形)

Arrays 在数学意义上来说属于多维空间向量，但在计算机物理存储结构上，我们任然可以将其当做一维数组。在此基础上理解 Arrays 的重构将给我们带来相当的便利性。使用方法 reshape()，我们可以对数组进行重构，使之变形：

In [29]: x2.reshape(3,20)
Out[29]:
array([[5, 2, 4, 7, 6, 8, 8, 1, 6, 7, 7, 8, 1, 5, 9, 8, 9, 4, 3, 0],
       [3, 5, 0, 2, 3, 8, 1, 3, 3, 3, 7, 0, 1, 9, 9, 0, 4, 7, 3, 2],
       [7, 2, 0, 0, 4, 5, 5, 6, 8, 4, 1, 4, 9, 8, 1, 1, 7, 9, 9, 3]])

事实上，只要是满足所有维度大小的乘积等于数组总元素大小，我们都可以用来作为变形的参数。在本例中，数组 x2 的元素总个数为 345=60，而我们 reshape 后的维度改为了 3*20，即 3 行，每行 20 个元素。

Arrays 的拼接

In [44]: grid = np.array([[1, 2, 3],
...: [4, 5, 6]]) # 创建一个新的数组
In [45]: np.concatenate([grid, grid]) # 维度相同的数组，可以直接进行拼接
Out[45]:
array([[1, 2, 3],
       [4, 5, 6],
       [1, 2, 3],
       [4, 5, 6]])

In [46]: np.concatenate([grid, grid], axis=1) # 沿着第二维度进行拼接
Out[46]:
array([[1, 2, 3, 1, 2, 3],
       [4, 5, 6, 4, 5, 6]])

In [47]: x = np.array([1, 2, 3])
In [48]: np.vstack([x, grid]) # 垂直方向拼接
Out[48]:
array([[1, 2, 3],
       [1, 2, 3],
       [4, 5, 6]])

In [49]: y = np.array([[99],
...: [99]])
...: np.hstack([grid, y]) # 水平方向拼接
Out[49]:
array([[ 1, 2, 3, 99],
       [ 4, 5, 6, 99]])

In [62]: np.dstack([x,x])# 沿着第三维度进行拼接
Out[62]:
array([[[1, 1],
        [2, 2],
        [3, 3]]])

Arrays 的分割

In [64]: x = [1, 2, 3, 99, 99, 3, 2, 1]
In [66]: np.split(x,[3,5]) # 遇到3或5则分割数组
Out[66]: [array([1, 2, 3]), array([99, 99]), array([3, 2, 1])] # 数组别分割成三个子数组

与拼接类似，分割有垂直分割 np.vsplit(),水平分割 np.hsplit() 和第三维度分割 np.dsplit()。而作为更通用的方法，或者是更高维度的操作，建议使用 concatenate() 或 splite() 方法时指定 axis 的值。举一拼接例子：

In [71]: grid = np.array([[9, 8, 7],[6, 5, 4]])

In [71]: y = np.array([[99],[99]])
...: np.hstack([grid, y])
Out[72]:
array([[ 9, 8, 7, 99],
       [ 6, 5, 4, 99]])

In [73]: np.concatenate([grid, y], axis=1) # 与上面的hstack方法得到的结果一致
Out[73]:
array([[ 9, 8, 7, 99],
       [ 6, 5, 4, 99]])

基本计算

Python 的动态特性导致原生数组的循环迭代访问特别慢，因为每个元素要单独执行动态解释，不管同一数组里的元素是否是相同的数据类型。

考虑这样的情形，我们有大量的数据需要进行相关计算，并且每个元素都参与其中，我们能使用的方法就是进行循环迭代数据。作为例子，假设我们对一百万个数字进行求导：

In [2]: big_array = np.random.randint(1, 100, size=1000000)
In [3]: def  fun_derivatives(arrays):
...: result = np.empty(len(arrays))
...: for i in  range(len(arrays)):
...: result[i] =  1.0  / arrays[i]
...: return resul
  
In [5]: %timeit fun_derivatives(big_array) # 使用%timeit进行时间测试
1.84 s ± 26.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

当前 cpu MHz : 1606.563，差不多是每秒钟计算 1.6*10^9 次。而在上述例子中，我们把 1.84s 的小数部分忽略，最快也才一秒钟计算 1.0*10^6 次,相对来说还是效率太低，相当于每次计算花费了近千个时钟周期。

In [6]: %timeit 1.0/big_array
3.41 ms ± 9.54 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

我们可以看出，使用 numpy 原生的计算模式，假使把 3.41ms 看做 10ms，最慢也是 1.0*10^8,较之前的循环快了两个数量级。上面两种方法可以保证计算结果是一致的，可以看出，原生的方法效率的提升是很明显的。numpy 为这样的计算提供了静态类型的操作接口，把循环计算直接推入 numpy 的基础编译层，以此来提高效率，并称之为矢量化操作。在 numpy 中矢量化操作是通过 Ufuncs 实现的。

使用 Ufuncs 的 Arrays 运算

1.算数运算

In [4]: x = np.arange(5)
In [5]: x +  8
Out[5]: array([ 8, 9, 10, 11, 12])

In [6]: x /  5
Out[6]: array([ 0. , 0.2, 0.4, 0.6, 0.8])

In [7]: x **  2
Out[7]: array([ 0, 1, 4, 9, 16])

In [8]: -(9*x-4)**3
Out[8]: array([ 64, -125, -2744, -12167, -32768])

以上相当于运算符重载（在 c#中有这样的概念），而在 python 中叫包装。我们亦可以调用 numpy 提供的方法 np.add(x, 2)，和 x + 2 是等价的，而这样的方法有：

2. 绝对值

使用 absolute 或 abs 进行求值。

In [9]: x = np.array([-2, -1, 0, 1, 2])
In [10]: np.abs(x)
Out[10]: array([2, 1, 0, 1, 2])

3. 三角函数

使用 sin、cos、arctan 等三角函数，也可进行运算。

4. 指数和对数

使用 exp、exp2、power、log2 等指数对数函数，也可进行运算。

5. 专业的 Ufuncs 功能

Ufuncs 还有更多的运算功能，如双曲线、位运算，比较运算等。但更多的科学函数功能在 scipy.special 库中可以找到。

In [1]: from scipy import special
In [2]: # 伽马函数 (广义阶乘) 和相关函数
...: x = [1, 5, 10]
...: print("gamma(x) =", special.gamma(x))
...: print("ln|gamma(x)| =", special.gammaln(x))
...: print("beta(x, 2) =", special.beta(x, 2))
gamma(x) = [ 1.00000000e+00  2.40000000e+01  3.62880000e+05]
ln|gamma(x)|  = [ 0. 3.17805383 12.80182748]
beta(x, 2) = [ 0.5  0.03333333  0.00909091]

#误差函数（高斯积分）
#它的补码和反码
In [5]: x = np.array([0, 0.3, 0.7, 1.0])
...: print("erf(x) =", special.erf(x))
...: print("erfc(x) =", special.erfc(x))
...: print("erfinv(x) =", special.erfinv(x))
erf(x) = [ 0. 0.32862676 0.67780119 0.84270079]
erfc(x) = [ 1. 0.67137324 0.32219881 0.15729921]
erfinv(x) = [ 0. 0.27246271 0.73286908 inf]

6. 高级 Ufuncs 特性

指定输出到

In [6]: x = np.arange(5)
...: y = np.empty(5)
...: np.multiply(x, 10, out=y) # 指定输出到y
Out[6]: array([ 0., 10., 20., 30., 40.])

In [7]: y = np.zeros(10)
...: np.power(2, x, out=y[::2]) # 将x的平方输出到y的偶数项上，可用`y[::2] = 2 ** x`替换
Out[7]: array([ 1., 2., 4., 8., 16.])
In [8]: y
Out[8]: array([ 1., 0., 2., 0., 4., 0., 8., 0., 16., 0.])

聚合

In [9]: x = np.arange(1, 6)
...: np.add.reduce(x) #和聚合
Out[9]: 15

In [10]: np.multiply.reduce(x) #积聚合
Out[10]: 120

累积（相当于存储了聚合的中间结果）

In [11]: np.add.accumulate(x)
Out[11]: array([ 1, 3, 6, 10, 15])

In [12]: np.multiply.accumulate(x)
Out[12]: array([ 1, 2, 6, 24, 120])

外积

In [13]: x = np.arange(1, 6)
...: np.multiply.outer(x, x)
Out[13]:
array([[ 1, 2, 3, 4, 5],
       [ 2, 4, 6, 8, 10],
       [ 3, 6, 9, 12, 15],
       [ 4, 8, 12, 16, 20],
       [ 5, 10, 15, 20, 25]])

广播

不知读者是否发现，在上述的例子中，存在这样的通用模型，即一个标量与一个矩阵的四则运算的结果还是一个形状相同的矩阵。例如，5+[1,2,3,4]=[6,7,8,9]。这一式子的等价式子为 [5,5,5,5]+[1,2,3,4]=[6,7,8,9]。这相当于变量 5，在沿着一维方向上进行了扩展延伸成向量[5,5,5,5]，然后再与矩阵[1,2,3,4]进行计算。在 numpy，这种自动延伸扩展的行为我们称之为广播。其真实本意是为了让不同形状的矩阵能够参与相互计算。

NumPy 中的广播遵循一套严格的规则来确定两个数组之间的交互：

1. 如果两个阵列的尺寸数不同，则尺寸较小的阵列的形状将在其前（左）侧填充 。
2. 如果两个数组的形状在任何维度上都不匹配，则该维度中形状等于 1 的数组将被拉伸以匹配其他形状。
3. 如果在任何维度中，大小不一致且都不等于 1，则会引发错误。

In [9]: a = np.ones((2,3))  # [[ 1.,  1.,  1.], [ 1.,  1.,  1.]]                                               
In [10]: b = np.arange(3)   # [0, 1, 2]                                                
In [11]: a+b                                                                    
Out[11]: 
array([[ 1.,  2.,  3.],
       [ 1.,  2.,  3.]])

根据规则一，小矩阵 b（3）将在左侧被填充为（1,3），根据规则二，（1,3）将被拉伸为（2,3），此时，a 与 b 便能发生计算。假使 b = np.arange(2) 或者 b = np.arange(4) 无论如何，都不可能与 a 形状一致，不满足规则三，则无法计算。

聚合计算

聚合计算旨在获取统计性的摘要内容，是对大量数据的一种概览性观测方法。比如海量数据中，相对于每一个具体值，我们可能更关心最大值、均值、标准差等。在 numpy 中提供了这样的功能。

求和

In [2]: L = np.random.random(100)
In [3]: np.sum(L) # 等价于L.sum()
Out[3]: 50.340203568224226

最小/大值

In [4]: np.min(L) # 等价于L.min()
Out[4]: 0.027618831239017205

In [5]: np.max(L) # 等价于L.max()
Out[5]: 0.9887438369150463

更多的功能函数参照下表（方法对多维数组适用）：