《scikit-learn机器学习常用算法原理及编程实战》—2.3.2 Numpy运算

2.3.2  Numpy运算

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

　　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 [7]: b = np.random.randint(1, 5, 20).reshape(4, 5)

　　In [8]: b

　　Out[8]:

　　array([[3, 3, 2, 3, 4],

　　       [3, 2, 2, 4, 3],

　　       [2, 4, 2, 1, 2],

　　       [2, 4, 4, 4, 2]])

　　In [9]: b * 3        # 数组和标量乘法

　　Out[9]:

　　array([[ 9,  9,  6,  9, 12],

　　       [ 9,  6,  6, 12,  9],

　　       [ 6, 12,  6,  3,  6],

　　       [ 6, 12, 12, 12,  6]])

　　使用Numpy的优点是运行速度会比较快，我们可以对比一下使用Python的循环与使用Numpy运算在效率上的差别，从Log里看到运行效率相差近100倍。

　　In [10]: c = np.arange(10000)

　　In [11]: %timeit c + 1

　　10000 loops, best of 3: 23.7 us per loop

　　In [12]: %timeit [i + 1 for i in c]

　　100 loops, best of 3: 2.61 ms per loop

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

　　In [16]: a = np.random.random_integers(1, 5, (5, 4))

　　In [17]: a

　　Out[17]:

　　array([[3, 1, 3, 2],

　　      [1, 1, 4, 2],

　　      [1, 3, 5, 4],

　　      [3, 3, 3, 4],

　　      [5, 1, 5, 1]])

　　In [23]: b = np.ones((5, 4), dtype=int)

　　In [24]: b

　　Out[24]:

　　array([[1, 1, 1, 1],

　　      [1, 1, 1, 1],

　　      [1, 1, 1, 1],

　　      [1, 1, 1, 1],

　　      [1, 1, 1, 1]])

　　In [25]: a + b            # 数组加法

　　Out[25]:

　　array([[4, 2, 4, 3],

　　      [2, 2, 5, 3],

　　      [2, 4, 6, 5],

　　      [4, 4, 4, 5],

　　      [6, 2, 6, 2]])

　　In [26]: c = np.random.random_integers(1, 5, (3, 4))

　　In [27]: c

　　Out[27]:

　　array([[5, 5, 4, 5],

　　      [2, 2, 1, 3],

　　      [5, 1, 1, 5]])

　　In [28]: d = np.random.random_integers(1, 5, (3, 4))

　　In [29]: d

　　Out[29]:

　　array([[3, 4, 5, 4],

　　      [4, 3, 4, 2],

　　      [1, 4, 5, 4]])

　　In [30]: c * d              # 数组相乘，逐元素相乘，不是矩阵内积运算

　　Out[30]:

　　array([[15, 20, 20, 20],

　　      [ 8,  6,  4,  6],

　　      [ 5,  4,  5, 20]])

　　需要注意的是，乘法是对应元素相乘，不是矩阵内积，矩阵内积使用的是np.dot()函数。

　　In [2]: a = np.random.random_integers(1, 5, (3, 2))

　　In [4]: b = np.random.random_integers(1, 5, (2, 3))

　　In [5]: a

　　Out[5]:

　　array([[3, 1],

　　      [2, 3],

　　      [5, 1]])

　　In [6]: b

　　Out[6]:

　　array([[5, 1, 2],

　　      [3, 1, 4]])

　　In [7]: np.dot(a, b)                    # 矩阵内积

　　Out[7]:

　　array([[18,  4, 10],

　　      [19,  5, 16],

　　      [28,  6, 14]])

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

　　In [30]: a = np.random.random_integers(1, 5, (5, 4))

　　In [31]: a

　　Out[31]:

　　array([[3, 1, 3, 2],

　　      [1, 1, 4, 2],

　　      [1, 3, 5, 4],

　　      [3, 3, 3, 4],

　　      [5, 1, 5, 1]])

　　In [32]: b = np.arange(4)

　　In [33]: b

　　Out[33]: array([0, 1, 2, 3])

　　In [34]: a + b              # 5 ? 4 数组与 1 ? 4 数组的加法，满足广播条件

　　Out[34]:

　　array([[3, 2, 5, 5],

　　      [1, 2, 6, 5],

　　      [1, 4, 7, 7],

　　      [3, 4, 5, 7],

　　      [5, 2, 7, 4]])

　　In [35]: c = np.arange(5)

　　In [36]: c

　　Out[36]: array([0, 1, 2, 3, 4])

　　In [37]: a + c              # 5 ? 4 数组与 1 ? 5 数组加法，不满足广播条件，报错

　　-----------------------------------------------------------------------

　　ValueError                                Traceback (most recent call last)

　　<ipython-input-37-ca57d551b7f3> in <module>()

　　----> 1 a + c

　　ValueError: operands could not be broadcast together with shapes (5,4) (5,)

　　In [38]: c = np.arange(5).reshape(5, 1)      # 转换为 5 ? 1 列向量

　　In [39]: c

　　Out[39]:

　　array([[0],

　　      [1],

　　      [2],

　　      [3],

　　      [4]])

　　In [40]: a + c            # 5 ? 4 数组与 5 ? 1 数组的加法，满足广播条件

　　Out[40]:

　　array([[3, 1, 3, 2],

　　      [2, 2, 5, 3],

　　      [3, 5, 7, 6],

　　      [6, 6, 6, 7],

　　      [9, 5, 9, 5]])

　　从上面的例子可以看到，符合广播的条件是两个数组必须有一个维度可以扩展，然后在这个维度上进行复制，最终复制出两个相同维度的数组，再进行运算。具体可参阅上面代码中的注释。作为广播的一个特例，当一个二维数组和一个标量进行运算时，实际上执行的也是广播机制，它有两个维度可扩展，先在行上进行复制，再在列上进行复制，最终复制出和待运算的二维数组维度相同的数组后，再进行运算。

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

　　In [42]: a = np.array([1, 2, 3, 4])

　　In [43]: b = np.array([4, 2, 2, 4])

　　In [44]: a == b

　　Out[44]: array([False,  True, False,  True], dtype=bool)

　　In [45]: a > b

　　Out[45]: array([False, False,  True, False], dtype=bool)

　　In [46]: (a == b).all()

　　Out[46]: False

　　In [47]: (a == b).any()

　　Out[47]: True

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

　　In [48]: a = np.arange(6)

　　In [49]: a

　　Out[49]: array([0, 1, 2, 3, 4, 5])

　　In [50]: np.cos(a)

　　Out[50]:

　　array([ 1.        ,  0.54030231, -0.41614684, -0.9899925 , -0.65364362,

　　        0.28366219])

　　In [52]: np.exp(a)

　　Out[52]:

　　array([   1.        ,    2.71828183,    7.3890561 ,   20.08553692,

　　         54.59815003,  148.4131591 ])

　　In [53]: np.sqrt(a)

　　Out[53]:

　　array([ 0.        ,  1.        ,  1.41421356,  1.73205081,  2.        ,

　　        2.23606798])

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

　　In [68]: a = np.random.random_integers(1, 5, 6)

　　In [69]: a

　　Out[69]: array([2, 1, 4, 5, 4, 1])

　　In [70]: a.sum()

　　Out[70]: 17

　　In [71]: a.mean()

　　Out[71]: 2.8333333333333335

　　In [72]: a.std()

　　Out[72]: 1.5723301886761005

　　In [73]: a.min()

　　Out[73]: 1

　　In [74]: a.max()

　　Out[74]: 5

　　In [75]: a.argmin()            # 最小值元素所在的索引

　　Out[75]: 1

　　In [76]: a.argmax()            # 最大值元素所在的索引

　　Out[76]: 3

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

　　In [77]: b = np.random.random_integers(1, 5, (6, 4))

　　In [78]: b

　　Out[78]:

　　array([[3, 2, 4, 2],

　　      [4, 5, 1, 1],

　　      [4, 4, 1, 4],

　　      [3, 3, 4, 4],

　　      [3, 2, 4, 5],

　　      [3, 5, 2, 5]])

　　In [85]: b.sum()

　　Out[85]: 78

　　In [86]: b.sum(axis=0)

　　Out[86]: array([20, 21, 16, 21])

　　In [87]: b.sum(axis=1)

　　Out[87]: array([11, 11, 13, 14, 14, 15])

　　In [88]: b.sum(axis=1).sum()

　　Out[88]: 78

　　In [94]: b.min(axis=1)

　　Out[94]: array([2, 1, 1, 3, 2, 2])

　　In [95]: b.argmin(axis=1)

　　Out[95]: array([1, 2, 2, 0, 1, 2])

　　In [96]: b.std(axis=1)

　　Out[96]:

　　array([ 0.8291562 ,  1.78535711,  1.29903811,  0.5       ,  1.11803399,

　　        1.29903811])

　　其中，axis参数表示坐标轴，0表示按行计算，1表示按列计算。需要特别注意的是，按列计算后，计算结果Numpy会默认转换为行向量。这个行为和MATLAB/Octave等数值计算软件有较大的差异。

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

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

图2-2  硬币***示意图

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

　　%matplotlib inline

　　import matplotlib.pyplot as plt

　　import numpy as np

　　n_person = 2000          # 总共创建 2000 组参加***

　　n_times = 500            # 抛硬币 500 次

　　t = np.arange(n_times)

　　# 创建只包含 -1 和 1 两种类型元素的数组来表示输赢序列

　　steps = 2 * np.random.random_integers(0, 1, (n_person, n_times)) - 1

　　amount = np.cumsum(steps, axis=1)            # 计算每组的输赢总额

　　sd_amount = amount ** 2                      # 计算平方

　　mean_sd_amount = sd_amount.mean(axis=0)      # 所有参加***的组求平均值

　　# 画出数据（绿色），同时画出平方根的曲线（红色）

　　plt.xlabel(r"$t$")

　　plt.ylabel(r"$\sqrt{\langle (\delta x)^2 \rangle}$")

　　plt.plot(t, np.sqrt(mean_sd_amount), 'g.', t, np.sqrt(t), 'r-');

　　上面代码可以直接在IPython notebook环境下运行，可以看到两根线基本重叠。即一维随机漫步算法的***法则，输赢的金额和***的次数之间基本呈现平方根曲线的分布，如图2-3所示。

图2-3  一维随机漫步算法

　　感兴趣的读者，可以参阅Wikipedia上二维随机漫步的情况en.wikipedia.org/wiki/ Random_walk，也可以思考一下如何使用Numpy来实现。

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

　　In [11]: a = np.arange(12)

　　In [12]: a

　　Out[12]: array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])

　　In [14]: b = a.reshape(4, 3)          # 转换为 4?3 的二维数组

　　In [15]: b

　　Out[15]:

　　array([[ 0,  1,  2],

　　      [ 3,  4,  5],

　　      [ 6,  7,  8],

　　      [ 9, 10, 11]])

　　In [16]: b.ravel()                        # 变为一维向量

　　Out[16]: array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])

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

　　In [20]: a = np.arange(4)

　　In [21]: a

　　Out[21]: array([0, 1, 2, 3])

　　In [22]: a.shape

　　Out[22]: (4,)

　　In [23]: b = a[:, np.newaxis]            # 在列上添加一个维度，变成 4?1 数组

　　In [24]: b

　　Out[24]:

　　array([[0],

　　      [1],

　　      [2],

　　      [3]])

　　In [25]: b.shape

　　Out[25]: (4, 1)

　　In [26]: c = a[np.newaxis, :]          # 在行上添加一个维度，变成 1?4 数组

　　In [27]: c

　　Out[27]: array([[0, 1, 2, 3]])

　　In [28]: c.shape

　　Out[28]: (1, 4)

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

　　In [36]: a = np.random.random_integers(1, 10, (6, 4))

　　In [37]: a

　　Out[37]:

　　array([[ 1,  4,  8, 10],

　　      [10,  9,  6,  2],

　　      [ 1,  4, 10,  5],

　　      [ 5,  7,  1,  1],

　　      [ 5,  2,  2,  8],

　　      [ 6, 10, 10,  7]])

　　In [40]: b = np.sort(a, axis=1)        # 按行独立排序，返回一个备份

　　In [41]: b

　　Out[41]:

　　array([[ 1,  4,  8, 10],

　　      [ 2,  6,  9, 10],

　　      [ 1,  4,  5, 10],

　　      [ 1,  1,  5,  7],

　　      [ 2,  2,  5,  8],

　　      [ 6,  7, 10, 10]])

　　In [42]: a.sort(axis=0)                # 按列排序，直接把结果保存到当前数组

　　In [43]: a

　　Out[43]:

　　array([[ 1,  2,  1,  1],

　　      [ 1,  4,  2,  2],

　　      [ 5,  4,  6,  5],

　　      [ 5,  7,  8,  7],

　　      [ 6,  9, 10,  8],

　　      [10, 10, 10, 10]])

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

　　In [52]: a = np.random.random_integers(1, 10, 6)

　　In [53]: a

　　Out[53]: array([10,  5,  3,  8,  4,  8])

　　In [55]: idx = a.argsort()

　　In [56]: idx

　　Out[56]: array([2, 4, 1, 3, 5, 0])

　　In [57]: a[idx]

　　Out[57]: array([ 3,  4,  5,  8,  8, 10])

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

　　In [80]: p = np.poly1d([1, -4, 3])      # 二阶多项式的系数

　　In [81]: p(0)                            # 自变量为 0 时的多项式的值

　　Out[81]: 3

　　In [82]: p.roots                            # 多项式的根

　　Out[82]: array([ 3.,  1.])

　　In [83]: p(p.roots)                        # 多项式根处的值肯定是 0

　　Out[83]: array([ 0.,  0.])

　　In [84]: p.order                            # 多项式的阶数

　　Out[84]: 2

　　In [85]: p.coeffs                        # 多项式的系数

　　Out[85]: array([ 1, -4,  3])

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

　　%matplotlib inline

　　import matplotlib.pyplot as plt

　　import numpy as np

　　n_dots = 20

　　n_order = 3

　　x = np.linspace(0, 1, n_dots)                    # [0, 1] 之间创建 20 个点

　　y = np.sqrt(x) + 0.2*np.random.rand(n_dots)

　　p = np.poly1d(np.polyfit(x, y, n_order))      # 用 3 阶多项式拟合

　　print p.coeffs

　　# 画出拟合出来的多项式所表达的曲线以及原始的点

　　t = np.linspace(0, 1, 200)

　　plt.plot(x, y, 'ro', t, p(t), '-');

　　运行效果如图2-4所示。

图2-4  多项式拟合

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

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

图2-5  蒙特卡罗方法计算圆周率

　　import numpy as np

　　# 假设圆的半径为1，圆心在原点

　　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        # 计算出 PI 的值

　　print pi

　　如果我们取一百万个点，计算出来的圆周率大概是3.142376，读者可以修改n_dots的值来看不同的点的个数对计算结果的精度影响。蒙特卡罗方法算法还有非常广泛的应用，感兴趣的读者可以在网络上搜索一下这个算法，详细了解其应用场景。

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

　　In [92]: a = np.arange(15).reshape(3,5)

　　In [93]: a

　　Out[93]:

　　array([[ 0,  1,  2,  3,  4],

　　      [ 5,  6,  7,  8,  9],

　　      [10, 11, 12, 13, 14]])

　　In [94]: np.savetxt('a.txt', a)

　　In [95]: ls

　　a.txt

　　In [96]: cat a.txt                  # 写入时把整型转换为浮点型了

　　0.000000000000000000e+00 1.000000000000000000e+00 2.000000000000000000e+00 

　　3.000000000000000000e+00 4.000000000000000000e+00

　　5.000000000000000000e+00 6.000000000000000000e+00 7.000000000000000000e+00 

　　8.000000000000000000e+00 9.000000000000000000e+00

　　1.000000000000000000e+01 1.100000000000000000e+01 1.200000000000000000e+01 

    1.300000000000000000e+01 1.400000000000000000e+01

　　In [97]: b = np.loadtxt('a.txt')

　　In [98]: b                          # 读出时也是浮点型

　　Out[98]:

　　array([[  0.,   1.,   2.,   3.,   4.],

　　      [  5.,   6.,   7.,   8.,   9.],

　　      [ 10.,  11.,  12.,  13.,  14.]])

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

　　In [116]: a

　　Out[116]:

　　array([[ 0,  1,  2,  3,  4],

　　      [ 5,  6,  7,  8,  9],

　　      [10, 11, 12, 13, 14]])

　　In [117]: np.save('a.npy', a)

　　In [118]: ls -l          # 新生成的 a.npy 与 a.txt 的大小不同

　　total 8

　　-rw-r--r-- 1 kamidox 200 Mar 27 21:27 a.npy

　　-rw-r--r-- 1 kamidox 375 Mar 27 21:00 a.txt

　　In [119]: c = np.load('a.npy')

　　In [120]: c                  # 读入的是整型，不是浮点型

　　Out[120]:

　　array([[ 0,  1,  2,  3,  4],

　　      [ 5,  6,  7,  8,  9],

　　      [10, 11, 12, 13, 14]])

　　本节只是Numpy的快速入门介绍，详细的信息可以参阅Numpy的官方网站www.numpy.  org。本节介绍的几个例子可参阅随书代码ch02.02.ipynb。