跟我一起学scikit-learn4：Numpy简介

Numpy是Python科学计算的基础库，主要提供了高性能的N维数组实现以及计算能力，还提供了和其他语言如C/C++集成的能力。此外还实现了一些基础的数学算法，如线性代数相关、傅里叶变换以及随机数生成等。

（1）Numpy数组

可以直接用Python列表来创建数组。

In [1]: import numpy as np
In [2]: a = np.array([1,2,3,4])
In [3]: a
Out[3]: array([1, 2, 3, 4])
In [4]: b = np.array([[1,2],[3,4],[5,6]])
In [5]: b
Out[5]:
array([[1, 2],
      [3, 4],
      [5, 6]])

可以查看array的属性，包括数据的维度和类型。

In [6]: b.ndim
Out[6]: 2
In [7]: b.shape
Out[7]: (3, 2)
In [8]: b.dtype
Out[8]: dtype('int64')

也可以使用Numpy提供的函数来创建数组。

In [9]: c = np.arange(10)
In [10]: c
Out[10]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [11]: d = np.linspace(0,2,11)
In [12]: d
Out[12]: array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ])
In [13]: np.ones((3,3))
Out[13]:
array([[1., 1., 1.],
      [1., 1., 1.],
      [1., 1., 1.]])
In [15]: np.zeros((3,6))
Out[15]:
array([[0., 0., 0., 0., 0., 0.],
      [0., 0., 0., 0., 0., 0.],
      [0., 0., 0., 0., 0., 0.]])
In [16]: np.eye(4)
Out[16]:
array([[1., 0., 0., 0.],
      [0., 1., 0., 0.],
      [0., 0., 1., 0.],
      [0., 0., 0., 1.]])
In [17]: np.random.randn(6,4)
Out[17]:
array([[ 1.31627668,  0.09790842, -0.38651855,  0.13092439],
      [ 0.30394956, -1.01969158, -1.34606586,  0.64042307],
      [-0.08192975,  0.95356347,  1.38385924, -1.74049954],
      [ 0.38581386,  0.03418643,  2.13147723, -1.93122617],
      [-1.00341116, -0.0980931 ,  1.59971491, -1.45187833],
      [-1.01829085, -1.33163169,  1.10783757,  0.05001029]])

Numpy提供了灵活的索引机制来访问数组内的元素。

In [18]: a = np.arange(10)
In [19]: a
Out[19]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [20]: a[0],a[3],a[-1]
Out[20]: (0, 3, 9)
In [21]: a[:4]
Out[21]: array([0, 1, 2, 3])
In [22]: a[3:7]
Out[22]: array([3, 4, 5, 6])
In [23]: a[6:]
Out[23]: array([6, 7, 8, 9])
In [24]: a[2:8:2]
Out[24]: array([2, 4, 6])
In [25]: a[2::2]
Out[25]: array([2, 4, 6, 8])
In [26]: a[::3]
Out[26]: array([0, 3, 6, 9])

二维数据的索引分成行和列两个维度，会更灵活一些。

In [27]: a = np.arange(0,51,10).reshape(6,1)+np.arange(6)
In [28]: a
Out[28]:
array([[ 0,  1,  2,  3,  4,  5],
      [10, 11, 12, 13, 14, 15],
      [20, 21, 22, 23, 24, 25],
      [30, 31, 32, 33, 34, 35],
      [40, 41, 42, 43, 44, 45],
      [50, 51, 52, 53, 54, 55]])
In [29]: a[0,0],a[2,-1]
Out[29]: (0, 25)
In [30]: a[0,2:5]
Out[30]: array([2, 3, 4])
In [31]: a[:3,3:]
Out[31]:
array([[ 3,  4,  5],
      [13, 14, 15],
      [23, 24, 25]])
In [32]: a[2,:]
Out[32]: array([20, 21, 22, 23, 24, 25])
In [33]: a[:,3]
Out[33]: array([ 3, 13, 23, 33, 43, 53])
In [34]: a[:,::2]
Out[34]:
array([[ 0,  2,  4],
      [10, 12, 14],
      [20, 22, 24],
      [30, 32, 34],
      [40, 42, 44],
      [50, 52, 54]])
In [35]: a[::2,::3]
Out[35]:
array([[ 0,  3],
      [20, 23],
      [40, 43]])

另外一个索引的方法是通过布尔数组。

In [36]: a = np.random.randint(10,20,6)
In [37]: a
Out[37]: array([17, 19, 12, 11, 18, 11])
In [38]: a%2==0
Out[38]: array([False, False,  True, False,  True, False])
In [39]: a[a%2==0]
Out[39]: array([12, 18])

Numpy总是试图自动地将结果转换为行向量。

In [40]: a = np.arange(0,51,10).reshape(6,1)+np.arange(6)
In [41]: a
Out[41]:
array([[ 0,  1,  2,  3,  4,  5],
      [10, 11, 12, 13, 14, 15],
      [20, 21, 22, 23, 24, 25],
      [30, 31, 32, 33, 34, 35],
      [40, 41, 42, 43, 44, 45],
      [50, 51, 52, 53, 54, 55]])
In [42]: a[a%2==0]
Out[42]:
array([ 0,  2,  4, 10, 12, 14, 20, 22, 24, 30, 32, 34, 40, 42, 44, 50, 52, 54])

在大部分情况下，Numpy数组是共享内存的，如果需要独立保存，需要显式地备份。可以使用np.may_share_memory()函数来判断两个数组是否共享内存。

In [43]: a = np.arange(6)
In [44]: a
Out[44]: array([0, 1, 2, 3, 4, 5])
In [45]: b = a[2:5]
In [46]: b
Out[46]: array([2, 3, 4])
In [47]: b[1] = 100
In [48]: b
Out[48]: array([  2, 100,  4])
In [49]: a
Out[49]: array([  0,  1,  2, 100,  4,  5])
In [50]: np.may_share_memory(a,b)
Out[50]: True
In [51]: b = a[2:6].copy()
In [52]: b
Out[52]: array([  2, 100,  4,  5])
In [53]: b[1] = 3
In [54]: b
Out[54]: array([2, 3, 4, 5])
In [55]: a
Out[55]: array([  0,  1,  2, 100,  4,  5])
In [56]: np.may_share_memory(a,b)
Out[56]: False

作为一个有趣的例子，我们使用埃拉托斯特尼筛法（Sieve of Eratosthenes）来打印[0,100]之间的所有质数。主要思路：从第一个质数2开始，数据里所有能被2整除的数字都不是质数，即从2开始，以2为步长，每跳经过的数字都能被2整除，把其标识为非质数。接着，从下一个质数3开始，重复上述过程。最终即可算出[0,100]之间的所有质数。

In [57]: import numpy as np
In [58]: a = np.arange(1,101)
In [59]: n_max = int(np.sqrt(len(a)))
In [60]: is_prime = np.ones(len(a),dtype=bool)
In [61]: is_prime[0] = False
In [62]: for i in range(2,n_max):
    ...:    if i in a[is_prime]:
    ...:        is_prime[(i**2-1)::i] = False
    ...:
In [63]: print(a[is_prime])
[ 2  3  5  7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97]

（2）Numpy运算

最简单的数值计算是数组和标量进行计算，计算过程是直接把数组里的元素和标量逐个进行计算。

In [1]: import numpy as np
In [2]: a = np.arange(6)
In [3]: a
Out[3]: array([0, 1, 2, 3, 4, 5])
In [4]: a + 5
Out[4]: array([ 5,  6,  7,  8,  9, 10])
In [5]: b = np.random.randint(1,5,20).reshape(4,5)
In [6]: b
Out[6]:
array([[1, 4, 1, 2, 3],
      [2, 2, 2, 4, 2],
      [3, 4, 2, 4, 4],
      [2, 2, 1, 1, 2]])
In [7]: b * 3
Out[7]:
array([[ 3, 12,  3,  6,  9],
      [ 6,  6,  6, 12,  6],
      [ 9, 12,  6, 12, 12],
      [ 6,  6,  3,  3,  6]])

使用Numpy的优点是运算速度会比较快，我们可以对比一下使用Python的循环与使用Numpy运算在效率上的差别，从打印的信息可以知道运行效率相差甚远。

In [8]: c = np.arange(10000)
In [9]: %timeit c + 1
4.47 µs ± 21.6 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [10]: %timeit [i+1 for i in c]
1.46 ms ± 2.35 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

另外一种是数组和数组的运算，如果数组的维度相同，那么数组中对应位置元素进行相应的数学运算：

In [15]: a = np.random.randint(1,5,(5,4))
In [16]: a
Out[16]:
array([[4, 1, 2, 2],
      [3, 3, 2, 1],
      [1, 3, 1, 3],
      [3, 2, 3, 4],
      [3, 1, 4, 3]])
In [17]: b = np.ones((5,4),dtype=int)
In [18]: b
Out[18]:
array([[1, 1, 1, 1],
      [1, 1, 1, 1],
      [1, 1, 1, 1],
      [1, 1, 1, 1],
      [1, 1, 1, 1]])
In [19]: a + b
Out[19]:
array([[5, 2, 3, 3],
      [4, 4, 3, 2],
      [2, 4, 2, 4],
      [4, 3, 4, 5],
      [4, 2, 5, 4]])

这里的数组相乘，是对应元素相乘，需要两个数组的维度一致，不是矩阵内积运算

In [20]: c = np.random.randint(1,5,(3,4))
In [21]: c
Out[21]:
array([[4, 2, 1, 2],
      [2, 1, 1, 3],
      [2, 3, 3, 2]])
In [22]: d = np.random.randint(1,5,(3,4))
In [23]: d
Out[23]:
array([[1, 3, 2, 3],
      [2, 3, 2, 4],
      [4, 3, 1, 1]])
In [24]: c * d
Out[24]:
array([[ 4,  6,  2,  6],
      [ 4,  3,  2, 12],
      [ 8,  9,  3,  2]])

矩阵内积使用np.dot()函数，A*B要满足矩阵内积的条件：A的列数=B的行数

In [25]: a = np.random.randint(1,5,(3,2))
In [26]: a
Out[26]:
array([[4, 3],
      [3, 1],
      [4, 2]])
In [27]: b = np.random.randint(1,5,(2,4))
In [28]: b
Out[28]:
array([[4, 2, 2, 4],
      [1, 4, 3, 1]])
In [29]: np.dot(a,b)
Out[29]:
array([[19, 20, 17, 19],
      [13, 10,  9, 13],
      [18, 16, 14, 18]])

如果数组的维度不同，则Numpy会试图使用广播机制来匹配，如果能匹配得上，就进行运算；如果不满足广播条件，则报错。

In [30]: a = np.random.randint(1,5,(5,4))
In [31]: a
Out[31]:
array([[2, 4, 2, 2],
      [3, 3, 3, 1],
      [3, 2, 1, 2],
      [4, 3, 3, 4],
      [3, 4, 4, 4]])
In [32]: b = np.arange(4)
In [33]: b
Out[33]: array([0, 1, 2, 3])
# 5x4数组和1x4数组的加法，满足广播条件
In [34]: a + b
Out[34]:
array([[2, 5, 4, 5],
      [3, 4, 5, 4],
      [3, 3, 3, 5],
      [4, 4, 5, 7],
      [3, 5, 6, 7]])
In [35]: c = np.arange(5)
# 5x4数组和1x5数组的加法，不满足广播条件
In [36]: a + c
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-36-e81e582b6fa9> in <module>()
----> 1 a + c
ValueError: operands could not be broadcast together with shapes (5,4) (5,)
In [37]: c = np.arange(5).reshape(5,1)
In [38]: c
Out[38]:
array([[0],
      [1],
      [2],
      [3],
      [4]])
# 5x4数组和5x1数组的加法，满足广播条件
In [39]: a + c
Out[39]:
array([[2, 4, 2, 2],
      [4, 4, 4, 2],
      [5, 4, 3, 4],
      [7, 6, 6, 7],
      [7, 8, 8, 8]])

从上面的例子可以看出，符合广播的条件是：两个数组必须有一个维度可以扩展，然后在这个维度上进行复制，最终复制出两个相同维度的数组，再进行运算。

作为广播的一个特例，当一个二维数组和一个标量进行运算时，实际上执行的也是广播机制，它有两个维度可以扩展，先在行上进行复制，再在列上进行复制，最红复制出和二维数组相同的维度，再进行运算。

数组还可以直接进行比较，返回一个同维度的布尔数组。针对布尔数组，可以使用all()/any()来返回布尔数组的标量值。

In [40]: a = np.array([1,2,3,4])
In [41]: b = np.array([4,2,2,4])
In [42]: a == b
Out[42]: array([False,  True, False,  True])
In [43]: (a == b).all()
Out[43]: False
In [44]: (a == b).any()
Out[44]: True

Numpy还提供了一些数组运算的内置函数。

In [45]: a = np.arange(6)
In [46]: a
Out[46]: array([0, 1, 2, 3, 4, 5])
In [47]: np.cos(a)
Out[47]:
array([ 1.        ,  0.54030231, -0.41614684, -0.9899925 , -0.65364362, 0.28366219])
In [48]: np.exp(a)
Out[48]:
array([  1.        ,  2.71828183,  7.3890561 ,  20.08553692,
        54.59815003, 148.4131591 ])
In [49]: np.sqrt(a)
Out[49]:
array([0.        , 1.        , 1.41421356, 1.73205081, 2.        , 2.23606798])

Numpy提供了一些基本的统计功能，包括求和、求平均值、求方差等。

In [50]: a = np.random.randint(1,5,6)
In [51]: a
Out[51]: array([3, 2, 2, 1, 4, 1])
In [52]: a.sum()
Out[52]: 13
In [53]: a.mean()
Out[53]: 2.1666666666666665
In [54]: a.std()
Out[54]: 1.0671873729054748
In [55]: a.min()
Out[55]: 1
In [56]: a.max()
Out[56]: 4
In [57]: a.argmin()    # 最小元素的索引
Out[57]: 3
In [58]: a.argmax()    # 最大元素的索引
Out[58]: 4

针对二维数组或者更高维的数组，可以根据行列来计算。

In [59]: b = np.random.randint(1,5,(6,4))
In [60]: b
Out[60]:
array([[4, 2, 4, 4],
      [1, 3, 1, 3],
      [3, 2, 4, 3],
      [2, 2, 3, 3],
      [2, 3, 1, 2],
      [2, 2, 2, 1]])
In [61]: b.sum()
Out[61]: 59
In [63]: b.sum(axis=0)
Out[63]: array([14, 14, 15, 16])
In [64]: b.sum(axis=1)
Out[64]: array([14,  8, 12, 10,  8,  7])
In [65]: b.sum(axis=1).sum()
Out[65]: 59
In [66]: b.min(axis=1)
Out[66]: array([2, 1, 2, 2, 1, 1])
In [67]: b.argmin(axis=1)
Out[67]: array([1, 0, 1, 0, 2, 3], dtype=int64)
In [68]: b.std(axis=1)
Out[68]:
array([0.8660254 , 1.        , 0.70710678, 0.5      , 0.70710678, 0.4330127 ])

其中，axis参数表示坐标轴，0表示按行计算，1表示按列计算。需要特别主要的是，按列计算后，计算结果Numpy会默认转换为行向量显示。

下面通过一个例子来看一下Numpy数值计算的应用，考察的是简单的一维随机漫步算法。比如，两个人用一个均匀的硬币来赌博，硬币抛出正面和反面的概率各占一半。硬币抛出正面时甲方输给乙方一块钱，反面时乙方输给甲方一块钱。我们来考察在这个赌博规则下，随着抛硬币次数的增加，输赢的总金额呈现怎样的分布。

若要解决这个问题，我们可以让足够多的人两两一组来参与这个游戏，然后抛足够多次数的硬币，就可以用统计的方法算出输赢的平均金额。当使用Numpy来实现时，生成多个由-1和1构成的足够长的随机数组，用来代表每次硬币抛出正面和反面的事件。这个二维数组中，每行表示一组参与赌博的人抛出正面和反面的事件序列，然后按行计算这个数组的累加和就是这组输赢的金额。

在实际计算时，先求出每组输赢金额的平方，再求平均值。最后把平方根的值用绿色的点画在二维坐标上，同时画出的红色曲线，来对比两组曲线的重合情况。

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
n_person = 2000
n_times = 500
t = np.arange(n_times)
steps = 2 * np.random.randint(2, size=(n_person, n_times)) - 1
amount = np.cumsum(steps, axis=1)
sd_amount = amount ** 2
mean_sd_amount = sd_amount.mean(axis=0)
plt.figure(figsize=(8, 6))
plt.xlabel(r"$t$", fontsize=16)
plt.tick_params(labelsize=12)
plt.ylabel(r"$\sqrt{\langle (\delta x)^2 \rangle}$", fontsize=24)
plt.plot(t, np.sqrt(mean_sd_amount), 'g.', t, np.sqrt(t), 'r-')

在前面，我们经常使用np.reshape()函数进行数组维度变换，而np.ravel()函数的功能则正好相反，它把多维数组“摊平”成一维向量。

In [1]: import numpy as np
In [2]: a = np.arange(12)
In [3]: a
Out[3]: array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
In [4]: b = a.reshape(4,3)
In [5]: b
Out[5]:
array([[ 0,  1,  2],
      [ 3,  4,  5],
      [ 6,  7,  8],
      [ 9, 10, 11]])
In [6]: b.ravel()
Out[6]: array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])

另外一个常用的方法是使用np.newaxis给数组添加一个维度。

In [7]: a = np.arange(4)
In [8]: a
Out[8]: array([0, 1, 2, 3])
In [10]: a.shape
Out[10]: (4,)
In [11]: b = a[:,np.newaxis]
In [12]: b
Out[12]:
array([[0],
      [1],
      [2],
      [3]])
In [13]: b.shape
Out[13]: (4, 1)
In [14]: c = a[np.newaxis,:]
In [15]: c
Out[15]: array([[0, 1, 2, 3]])
In [16]: c.shape
Out[16]: (1, 4)

Numpy提供了数组排序的功能，可以按行单独排序，也可以按列单独排序。排序时可以返回一个备份，可以直接把排序后的结果保存在当前数组里。

In [17]: a = np.random.randint(1,10,(6,4))
In [18]: a
Out[18]:
array([[4, 4, 3, 3],
      [2, 7, 3, 9],
      [1, 1, 4, 3],
      [3, 2, 4, 7],
      [6, 1, 2, 7],
      [3, 3, 5, 4]])
In [19]: b = np.sort(a,axis=1)
In [20]: b
Out[20]:
array([[3, 3, 4, 4],
      [2, 3, 7, 9],
      [1, 1, 3, 4],
      [2, 3, 4, 7],
      [1, 2, 6, 7],
      [3, 3, 4, 5]])
In [21]: a.sort(axis=0)
In [22]: a
Out[22]:
array([[1, 1, 2, 3],
      [2, 1, 3, 3],
      [3, 2, 3, 4],
      [3, 3, 4, 7],
      [4, 4, 4, 7],
      [6, 7, 5, 9]])

还可以直接计算出排序后的索引，利用排序后的索引可以直接获取到排序后的数组。

In [23]: a = np.random.randint(1,10,6)
In [24]: a
Out[24]: array([6, 6, 3, 1, 8, 4])
In [25]: idx = a.argsort()
In [26]: idx
Out[26]: array([3, 2, 5, 0, 1, 4], dtype=int64)
In [27]: a[idx]
Out[27]: array([1, 3, 4, 6, 6, 8])

Numpy的高级功能包括多项式求解以及多项式拟合的功能。下面的代码构建了一个二阶多项式：

In [28]: p = np.poly1d([1,-4,3])      # 由一组系数构建一个二阶多项式
In [29]: p(0)                         # x=0时多项式的值
Out[29]: 3
In [30]: p.roots                      # p(x) = 0的根
Out[30]: array([3., 1.])
In [31]: p(p.roots)                   # 根的多项式值=0
Out[31]: array([0., 0.])
In [32]: p.order                      # 多项式的阶数
Out[32]: 2
In [33]: p.coeffs                     # 多项式的系数
Out[33]: array([ 1, -4,  3])

Numpy提供的np.polyfit()函数可以用多项式对数据进行拟合。在下面的例子里，我们生成20个在平方根曲线周围引入随机噪声的点，用一个3阶多项式来拟合这些点。

n_dots = 20
n_order = 3
x = np.linspace(0, 1, n_dots)
y = np.sqrt(x) + 0.2*np.random.rand(n_dots)
p = np.poly1d(np.polyfit(x, y, n_order))
print(p.coeffs)
t = np.linspace(0, 1, 200)
plt.plot(x, y, 'ro', t, p(t), '-')

另外一个有意思的例子是，使用Numpy求圆周率的值。使用的算法是蒙特卡罗方法（Monte Carlo method），其主要思想是，利用正方形的边长画出一个1/4圆的扇形，假设圆的半径为，则正方形的面积是，圆的面积是，它们的面积之比是。

我们在正方形内随机产生足够多的点，计算落在扇形区域内的点的个数与总的点的个数的比值。当产生的随机点足够多时，这个比值和面积比值应该是一致的。这样我们就可以算出的值。判断一个点是否落在扇形区域的方法是计算这个点到圆心的距离，当距离小于半径时，说明落在了扇形区域内。

n_dots = 1000000
x = np.random.random(n_dots)
y = np.random.random(n_dots)
distance = np.sqrt(x ** 2 + y ** 2)
in_circle = distance[distance < 1]
pi = 4 * float(len(in_circle)) / n_dots
print(pi)

最后介绍一下Numpy文件的相关操作。Numpy数组作为文本文件，可以直接保存到文件系统中，也可以直接从文件系统读出数据：

In [1]: import numpy as np
In [2]: a = np.arange(15).reshape(3,5)
In [3]: a
Out[3]:
array([[ 0,  1,  2,  3,  4],
      [ 5,  6,  7,  8,  9],
      [10, 11, 12, 13, 14]])
In [4]: np.savetxt('a.txt',a)   # 写入时把整型转换为浮点型
In [6]: ls 
a.txt
In [8]: b = np.loadtxt('a.txt')
In [9]: b
Out[9]:
array([[ 0.,  1.,  2.,  3.,  4.],
      [ 5.,  6.,  7.,  8.,  9.],
      [10., 11., 12., 13., 14.]])

保存为文本格式的可读性好，但性能较低。也可以直接保存为Numpy特有的二进制格式：

In [10]: a
Out[10]:
array([[ 0,  1,  2,  3,  4],
      [ 5,  6,  7,  8,  9],
      [10, 11, 12, 13, 14]])
In [11]: np.save('a.npy',a)
In [12]: ls -l a.txt a.npy      # 生成的a.npy和a.txt的大小不一样
2019/05/13  23:57  188 a.npy
2019/05/13  23:53  378 a.txt
In [13]: c = np.load('a.npy')   # 读入的仍然是整型
In [14]: c
Out[14]:
array([[ 0,  1,  2,  3,  4],
      [ 5,  6,  7,  8,  9],
      [10, 11, 12, 13, 14]])