python科学计算包numpy使用心得

1 数据结构

numpy使用一种称为ndarray的类似Matlab的矩阵式数据结构管理数据，比python的列表和标准库的array类更为强大，处理数据更为方便。

1.1 数组的生成

在numpy中，生成数组需要指定数据类型，默认是int32，即整数，可以通过dtype参数来指定，一般用到的有int32、bool、float32、uint32、complex，分别代表整数、布尔值、浮点型、无符号整数和复数

一般而言，生成数组的方法有这么几种：

以list列表为参数生成（用tolist方法即可转换回list）：

In[3]: a = array([1, 2, 3])
In[4]: a
Out[4]: array([1, 2, 3])
In[5]: a.tolist()
Out[5]: [1, 2, 3]

指定起点、终点和步长生成等差序列或等比数列：

In[7]: a = arange(1, 10, 2)
In[8]: a
Out[8]: array([1, 3, 5, 7, 9])

In[13]: a = linspace(0, 10, 5)
In[14]: a
Out[14]: array([  0. ,   2.5,   5. ,   7.5,  10. ])

In[148]: a = logspace(0, 3, 10) # 0表示起点为10^0，3表示起点为10^3，基数通过base参数指定
In[149]: a
Out[148]: 
array([    1.        ,     2.15443469,     4.64158883,    10.        ,
          21.5443469 ,    46.41588834,   100.        ,   215.443469  ,
         464.15888336,  1000.        ])

从迭代器中生成：

In[17]: iter = (i for i in range(5))
In[18]: a = fromiter(iter, dtype=int32)
In[19]: a
Out[19]: array([0, 1, 2, 3, 4])

从函数中生成：

In[156]: def f(i, j):
...     return abs(i-j)
...     
In[157]: fromfunction(f, (4, 4))
Out[156]: 
array([[ 0.,  1.,  2.,  3.],
       [ 1.,  0.,  1.,  2.],
       [ 2.,  1.,  0.,  1.],
       [ 3.,  2.,  1.,  0.]])

还可以用zeros、ones、empty等函数快速创建数组。

矩阵视为二维数组：

In[24]: b = array([arange(5), arange(1, 6), arange(2, 7)])
In[25]: b
Out[25]: 
array([[0, 1, 2, 3, 4],
       [1, 2, 3, 4, 5],
       [2, 3, 4, 5, 6]])

根据相同的方法可以拓展到更高维。

另外，我们还可以生成自定义数据格式的数组（称为结构数组），用来记录电子表格或数据库中一行数据的信息：

In[61]: t = dtype([('name', str, 40), ('number', int32), ('score', float32)])
In[62]: t
Out[62]: dtype([('name', '<U40'), ('number', '<i4'), ('score', '<f4')])
In[63]: students = array([('Tom', 10, 80), ('Jenny', 11, 90.5), ('Mike', 9, 98.5)], dtype=t)
In[64]: students
Out[64]: 
array([('Tom', 10, 80.0), ('Jenny', 11, 90.5), ('Mike', 9, 98.5)], 
      dtype=[('name', '<U40'), ('number', '<i4'), ('score', '<f4')])
In[65]: students[1]
Out[65]: ('Jenny', 11, 90.5)

后面我们会看到pandas提供了一种更精致的方法处理记录。

1.2 数组的索引

简单的下标索引：

In[30]: a[2]
Out[30]: 2
In[31]: b[2, 1]
Out[31]: 3

与python一样，索引的起点为0。负数的索引当然也是可以的：

In[32]: a[-1]
Out[32]: 4
In[33]: b[-1, -2]
Out[33]: 5

以整数数组为下标索引，一次性索引多个值：

In[162]: arange(11, 20)[array([2, 4, 8])]
Out[161]: array([13, 15, 19])

还可以通过布尔值来索引：

In[40]: idx = array([True, False, False, True, True])
In[41]: a[idx]
Out[41]: array([0, 3, 4])

这可以应用在高级索引中，比如条件索引：

b[b>3]
Out[42]: array([4, 4, 5, 4, 5, 6])

得到b中所有大于3的元素，以array形式返回，我们能这么写的原因是b>3会返回一个布尔数组，形式与b一致，各位置的值是b中各元素与3比较之后的结果：

In[43]: b>3
Out[43]: 
array([[False, False, False, False,  True],
       [False, False, False,  True,  True],
       [False, False,  True,  True,  True]], dtype=bool)

1.3 数组的切片

ndarray数组支持各种形式的切片，既可以以下标为线索，还可以以值为线索，为了区分二者，重新生成一个数组：

a = arange(11, 20)
In[54]: a
Out[54]: array([11, 12, 13, 14, 15, 16, 17, 18, 19])

根据下标切片：

In[55]: a[1:4]
Out[55]: array([12, 13, 14])
In[56]: a[1:8:2]
Out[56]: array([12, 14, 16, 18])
In[57]: a[1::2]
Out[57]: array([12, 14, 16, 18])
In[58]: a[:8:]
Out[58]: array([11, 12, 13, 14, 15, 16, 17, 18])

方括号中三个参数为别是起点、终点和步长，默认值分别是0、-1、1，注意终点是不被包含的。可以简单地令步长为-1来翻转数组：

In[60]: a[::-1]
Out[60]: array([19, 18, 17, 16, 15, 14, 13, 12, 11])

ndarray也支持多维数组的切片，先生成一个三维数组，可以通过修改一维数组的shape属性或调用其reshape方法来生成：

In[68]: a = arange(0, 24).reshape(2, 3, 4)
In[69]: a
Out[69]: 
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]]])

多维数组的索引其实跟一维区别不大，可以用:代表选取所有：

In[70]: a[:, 0, 0]
Out[70]: array([ 0, 12])
In[71]: a[0, :, 0]
Out[71]: array([0, 4, 8])
In[72]: a[0, 0, :]
Out[72]: array([0, 1, 2, 3])
In[73]: a[0, 0:2, 0:3]
Out[73]: 
array([[0, 1, 2],
       [4, 5, 6]])

多个冒号还可以用...来代替：

In[74]: a[...,3]
Out[74]: 
array([[ 3,  7, 11],
       [15, 19, 23]])

最后，可以使用slice对象来表示切片，它与用1:10:2形式产生切片类似：

In[169]: idx = slice(None, None, 2)
In[171]: a[idx,idx,idx]
Out[170]: 
array([[[ 0,  2],
        [ 8, 10]]])

相当于a[::2, ::2, ::2]

1.4 数组的变换

可以将上述三维数组展平：

In[75]: a.flatten()
Out[75]: 
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])

转置：

In[77]: b.transpose()
Out[77]: 
array([[0, 1, 2],
       [1, 2, 3],
       [2, 3, 4],
       [3, 4, 5],
       [4, 5, 6]])

修改shape属性来改变维度：

In[79]: a.shape = 4, 6
In[80]: a
Out[80]: 
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]])

1.5 数组的组合

首先创建一个与a同大小的数组：

In[83]: b = 2*a

• 1

可以进行多种方式组合，如水平组合：

In[88]: hstack((a, b))
Out[88]: 
array([[ 0,  1,  2,  3,  4,  5,  0,  2,  4,  6,  8, 10],
       [ 6,  7,  8,  9, 10, 11, 12, 14, 16, 18, 20, 22],
       [12, 13, 14, 15, 16, 17, 24, 26, 28, 30, 32, 34],
       [18, 19, 20, 21, 22, 23, 36, 38, 40, 42, 44, 46]])

垂直组合：

In[89]: vstack((a, b))
Out[89]: 
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],
       [ 0,  2,  4,  6,  8, 10],
       [12, 14, 16, 18, 20, 22],
       [24, 26, 28, 30, 32, 34],
       [36, 38, 40, 42, 44, 46]])

用concatenate函数可以同时实现这两种方式，通过指定axis参数，默认为0，使用垂直组合。

还可以进行深度组合：

In[91]: dstack((a, b))
Out[91]: 
array([[[ 0,  0],
        [ 1,  2],
        [ 2,  4],
        [ 3,  6],
        [ 4,  8],
        [ 5, 10]],

       [[ 6, 12],
        [ 7, 14],
        [ 8, 16],
        [ 9, 18],
        [10, 20],
        [11, 22]],

       [[12, 24],
        [13, 26],
        [14, 28],
        [15, 30],
        [16, 32],
        [17, 34]],

       [[18, 36],
        [19, 38],
        [20, 40],
        [21, 42],
        [22, 44],
        [23, 46]]])

就好像将两张二维平面的点数据沿纵轴方向叠在一起一样。

1.6 数组的分割

水平分割：

In[94]: hsplit(a, 3)
Out[94]: 
[array([[ 0,  1],
        [ 6,  7],
        [12, 13],
        [18, 19]]), array([[ 2,  3],
        [ 8,  9],
        [14, 15],
        [20, 21]]), array([[ 4,  5],
        [10, 11],
        [16, 17],
        [22, 23]])]

垂直分割：

In[97]: vsplit(a, 2)
Out[96]: 
[array([[ 0,  1,  2,  3,  4,  5],
        [ 6,  7,  8,  9, 10, 11]]), array([[12, 13, 14, 15, 16, 17],
        [18, 19, 20, 21, 22, 23]])]

用split函数可以同时实现这两个效果，通过设置其axis参数区别。

类似地，可以通过函数dsplit进行深度分割。

另外可以使用ndarray的一些属性来查看数组的信息：

In[125]: a.ndim # 维数
Out[124]: 2
In[126]: a.size # 元素总个数
Out[125]: 24
In[127]: a.itemsize # 元素在内存中所占的字节
Out[126]: 4
In[128]: a.shape # 维度
Out[127]: (4, 6)
In[130]: a.T # 转置，相当于transponse函数
Out[129]: 
array([[ 0,  6, 12, 18],
       [ 1,  7, 13, 19],
       [ 2,  8, 14, 20],
       [ 3,  9, 15, 21],
       [ 4, 10, 16, 22],
       [ 5, 11, 17, 23]], dtype=int32)

另外多维数组的flat属性可以给出一个”扁平迭代器“——flatiter对象，使我们能像一维数组一样迭代高维数组：

In[134]: for item in array([1, 2, 3, 4]).reshape(2, 2).flat:
...     print(item)
...     
1
2
3
4

flatiter对象可以直接获取多个元素，并直接赋值修改：

In[140]: af = a.flat
In[141]: af[:]
Out[140]: 
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], dtype=int32)
In[143]: af[3] = 15
In[144]: af[:]
Out[143]: 
array([ 0,  1,  2, 15,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23], dtype=int32)

1.7 矩阵的生成

上面提到了可以用二维数组来模拟矩阵，其实，numpy专门提供了一种用于处理矩阵的数据结构——matrix，它通过mat函数构造生成：

In[8]: m = mat('1 2 3; 4 5 6; 7 8 9')
In[9]: m
Out[9]: 
matrix([[1, 2, 3],
        [4, 5, 6],
        [7, 8, 9]])

二维数组与矩阵可以很方便地相互转换：

In[11]: array(m)
Out[11]: 
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
In[12]: mat(_)
Out[12]: 
matrix([[1, 2, 3],
        [4, 5, 6],
        [7, 8, 9]])

用matrix来处理矩阵更方便，有更多方法以供使用，如：

求逆：

In[17]: m.I
Out[17]: 
matrix([[ -4.50359963e+15,   9.00719925e+15,  -4.50359963e+15],
        [  9.00719925e+15,  -1.80143985e+16,   9.00719925e+15],
        [ -4.50359963e+15,   9.00719925e+15,  -4.50359963e+15]])

分块矩阵：

In[25]: I = eye(3)
In[26]: bmat('m I; I m')
Out[26]: 
matrix([[ 1.,  2.,  3.,  1.,  0.,  0.],
        [ 4.,  5.,  6.,  0.,  1.,  0.],
        [ 7.,  8.,  9.,  0.,  0.,  1.],
        [ 1.,  0.,  0.,  1.,  2.,  3.],
        [ 0.,  1.,  0.,  4.,  5.,  6.],
        [ 0.,  0.,  1.,  7.,  8.,  9.]])

2 数据处理

2.1 条件判断和搜索

用where函数可以得到满足条件的索引，便于后期处理：

In[219]: a = arange(24).reshape(4, 6)
In[220]: where(a>8)
Out[219]: 
(array([1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3], dtype=int32),
 array([3, 4, 5, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5], dtype=int32))

用compress函数可以筛选出一维数组中满足条件的值：

In[28]: a[0, :].compress(a[0, :] > 2)
Out[28]: array([3, 4, 5])

2.2 CSV文件读写

CSV（逗号分割值）格式可以简单方便地保存数组或矩阵。相比于python的pickle方法，保存为CSV文件可以用一般文本编辑器随时打开查看。保存和读取CSV文件都很简单。

In[190]: b
Out[189]: 
array([[ 0,  2,  4,  6,  8, 10],
       [12, 14, 16, 18, 20, 22],
       [24, 26, 28, 30, 32, 34],
       [36, 38, 40, 42, 44, 46]])
In[191]: savetxt("b.txt", b, delimiter=",")
In[192]: b1, b2 = loadtxt("b.txt", delimiter=",", usecols=(3, 4), unpack=True)
In[193]: b1, b2
Out[192]: (array([  6.,  18.,  30.,  42.]), array([  8.,  20.,  32.,  44.]))

保存时参数delimiter可选，用来分隔数组各元素，读取时也要相应地指定这个值，读取时也可只读取部分数据，usecols即用来指定选取的列，unpack设置为True时表示将这些列分开存储。

读写时遇到字符串（如时间）可以通过指定参数converters来转换。

In[252]: 
def datestr2num(s):
    return datetime.datetime.strptime(str(s, encoding="utf-8"), "%Y-%m-%d").date().weekday()

weeks, numbers = loadtxt("b.txt", converters={0:datestr2num}, unpack=True)
In[253]: weeks
Out[252]: array([ 2.,  4.])

2.3 通用函数

用frompyfunc函数可以将一个作用在单一数值的函数映射到作用在数组上的函数：

In[49]: def f(i):
...     return 2*i
...     
In[50]: ff = frompyfunc(f, 1, 1)
In[52]: ff(a)
Out[52]: 
array([[0, 2, 4, 6, 8, 10],
       [12, 14, 16, 18, 20, 22],
       [24, 26, 28, 30, 32, 34],
       [36, 38, 40, 42, 44, 46]], dtype=object)

frompyfunc的两个参数分别定义输入参数和输出参数的个数

另外，numpy提供了一些常用的通用函数，如针对加减乘除的add、subtract、multiply和divide。通用函数都有四个方法：reduce、accumulate、reduceat和outer，以add函数为例：

In[64]: add.reduce(a[0, :])
Out[64]: 15
In[65]: add.accumulate(a[0,:])
Out[65]: array([ 0,  1,  3,  6, 10, 15], dtype=int32)
In[69]: add.reduceat(a[0, :], [0, 5, 2, 4])
Out[69]: array([10,  5,  5,  9], dtype=int32)
In[70]: add.outer(a[0, :], a[1, :])
Out[70]: 
array([[ 6,  7,  8,  9, 10, 11],
       [ 7,  8,  9, 10, 11, 12],
       [ 8,  9, 10, 11, 12, 13],
       [ 9, 10, 11, 12, 13, 14],
       [10, 11, 12, 13, 14, 15],
       [11, 12, 13, 14, 15, 16]])

可见，reduce是将通用函数递归作用在所有元素上，得到最后结果；accumulate也是递归作用在所有元素上，不过它保留中间结果并返回；reduceat则根据指定的起始点进行累积运算，如果终点小于起点，则返回终点处的值；最后outer则是对两个输入数组的所有元素组合进行运算。

3 科学计算

3.1 统计分析

3.1.1 基本统计分析

average函数可以非常方便地计算加权平均值，或者用mean计算算术平均值：

In[204]: a = array([1, 2])
In[205]: average(a, weights=[1,2])
Out[204]: 1.6666666666666667

基本统计分析函数整理如下：

• 中位数：median
• 方差：var
• 标准差：std
• 差分：diff
• 最值：max、min、argmax、argmin（后两个得到最值所在的下标）

3.1.2 随机过程分析

3.2 线性代数

先生成一个各元素是0～1之内的随机数的矩阵：

In[47]: a = mat(fromiter((random.random() for i in range(9)), dtype = float32).reshape(3, 3))
In[48]: a
Out[48]: 
matrix([[ 0.45035544,  0.53587919,  0.57240343],
        [ 0.54386997,  0.16267321,  0.97020519],
        [ 0.6454953 ,  0.38505632,  0.94705021]], dtype=float32)

接下我们可以对它进行各种线性代数的操作， 如：

求逆：

In[49]: a.I
Out[49]: 
matrix([[-10.71426678, -14.01229095,  20.83065987],
        [  5.42686558,   2.7832334 ,  -6.13131571],
        [  5.09620285,   8.41894722, -10.64905548]], dtype=float32)

解线性方程组（用点积验证了结果）：

In[59]: b = fromiter((random.random() for i in range(3)), dtype = float32)
In[60]: b
Out[60]: array([ 0.56506187,  0.99419129,  0.70462942], dtype=float32)
In[61]: linalg.solve(a, b)
Out[61]: array([-5.3072257 ,  1.51327574,  3.74607611], dtype=float32)
In[63]: dot(a, _)
Out[63]: matrix([[ 0.56506193,  0.99419105,  0.70462948]], dtype=float32)

求特征值和特征向量：

In[64]: linalg.eig(a)
Out[64]: 
(array([ 1.78036737, -0.08517434, -0.13511421], dtype=float32),
 matrix([[-0.5075314 , -0.82206506,  0.77804375],
         [-0.56222379,  0.4528676 , -0.57155234],
         [-0.65292901,  0.34513769, -0.26072171]], dtype=float32))

行列式：

In[81]: linalg.det(a)
Out[81]: 0.020488938