作者：方跃文
Email: fyuewen@gmail.com
时间：始于2017年9月12日， 结束写作于
第四章笔记始于2017年10月17日，结束于2018年1月6日
时间： 2017年10月17日早晨
Numpy，即 numerical python的简称，是高性能科学计算和数据分析的基础包。它是本书所介绍的几乎所有高级工具的构建基础。其部分功能如下：
ndarray，一个具有矢量算数运算和复杂广播能力的快速且节省空间的多维数组
在不需要循环的情况下，用于对数组快速运算的标准数学函数
用于读写磁盘数据的工具以及用于操作内存映射文件的工具
线性代数、随机数生成以及傅里叶变化
用于集成由 C、C++、Fortran 等语言编写的代码的工具
Numpy 本身功能不复杂，但是理解 Numpy 有助于更高效地使用诸如 Pandas 之类的工具。
原书作者主要从事数据分析，所以他关注的功能主要集中于：
用于数据整理和清理、子集构造和过滤、转换等快速的矢量化数组运算
常用的数组算法，如排序、唯一化、集合运算等。
高效地描述统计和数据聚合/摘要运算
用于异构数据集的合并/连接运算的数据和关系型数据运算
将条件逻辑表述为数组表达式（而不是带有if-elif-else分支的循环）
数据的分组运算（聚合、转换、函数应用等）第五章将对此进行详细解释。
注：建议总是使用 import numpy as np； 而不是用 from numpy import *
时间： 2017年10月18日晚
Numpy 一个重要特点就是其 N 维数组对象，即 ndarray，该对象是一个快速而灵活的数据集容器。我们可以利用这种数组对整块数据进行一些运算，它的语法跟标量元素之间的运算相同：
In [2]:
import numpy.random as nrandom
data = nrandom.randn(3,2)
In [3]:
data
Out[3]:
In [4]:
data*10
Out[4]:
In [5]:
data + data
Out[5]:
ndarray 是 同构数据多维容器，that is to say, 所有元素必须是同类型的。
每个数组都有一个 shape （一个表示各维度大小的元祖）和一个 dtype （一个用于说明数组数据类型的对象）：
In [6]:
data.shape # 数组的维数，即行数和列数
Out[6]:
In [7]:
data.dtype #数组中元素的类型
Out[7]:
In [9]:
data.size #数组的大小
Out[9]:
In [14]:
dataconversion = data.astype('int8')
print('data is: ', data)
print('\n dataconversion is ', dataconversion)
虽然大多数数据分析工作不需要深入理解Numpy，但是精通面向数组的编程和思维方式是成为 Python 科学计算达人的一大步骤。
注意：第一版翻译版本中有个批注，说“本书中的数组、Numpy数组、ndarray 基本指的都是同一样东西，即 ndarray 对象”
In [15]:
import numpy as np
data1 = [2,3,3,5,6,9]
array1 = np.array(data1)
In [16]:
print('data1 is ', type(data1))
print('array1 is ', type(array1))
In [18]:
data1[:]
Out[18]:
In [21]:
array1
Out[21]:
In [22]:
print(array1)
print(array1.dtype)
print(array1.shape)
print(array1.size)
嵌套序列（比如由一组等长列表组成的列表），将会被转换为一个多维数组：
In [23]:
import numpy as np
data2=[[23,5,5,6], [4,56,2,8],[3,5,6,7],[2,3,4,5]]
arr2=np.array(data2)
In [103]:
arr2
Out[103]:
In [104]:
arr2.ndim #Number of array dimensions.
Out[104]:
In [105]:
arr2.shape
Out[105]:
In [24]:
arr2.size
Out[24]:
除非显示说明，np.array 会尝试为新建的这个数组推断出一个较为合适的数据类型。数据类型保存在一个特殊的 dtype 对象中。比如说，在上面的两个examples中，我们有
In [106]:
data.dtype
Out[106]:
In [107]:
arr2.dtype
Out[107]:
除 np.array 之外，还有一些函数可以新建数组。比如，zeros 和 ones 分别可创建指定长度或形状的全 0 和 全 1 数组。empty 可创建一个没有任何具体值的数组。要用这些方法创建多维数组，只需要传入一个表示形状的元祖即可：
In [108]:
np.zeros(10)
Out[108]:
In [109]:
arr4 = np.zeros((3,6,3))
arr4
Out[109]:
In [110]:
arr4.ndim
Out[110]:
In [111]:
arr3 = np.empty((2,4,2))
arr3
Out[111]:
In [112]:
arr3.ndim
Out[112]:
In [113]:
arr5 = np.empty((2,3,4,2))
arr5
Out[113]:
警告 认为 np.emptry 会返回全 0 数组的想法是不安全的。很多情况下（如上所示），它返回的都是一些未初始化的垃圾值。
arange 是 Python 内置函数range 的数组版：
In [114]:
np.arange(15)
Out[114]:
In [115]:
np.arange(2)
Out[115]:
下表列出了一些数组创建函数。由于Numpy关注的是数值计算，因此，如果没有特别的制定，数据类型一般都是 float64。
| 函数 | 说明 |
|---|---|
| array | 将输入数据(列表、元祖、数字或者其他数据类型)转换为 ndarray。要么推断出 dtype，要么显示地指定dtype。默认直接复制输入数据 |
| asarray | 将输入转为 ndarray，如果输入本身就是一个ndarray就不进行复制 |
| arange | 类似于python内置的range,但是返回的是一个ndarray,而不是一个列表 |
| ones、ones_like | 根据指定的形状和dtype创建一个全1数组。ones_like以另一个数组为参数，并根据其形状和dtype创建一个全1数组 |
| zeros、zeros_like | 类似上述命令，只是改为全0数组 |
| empty、empty_like | 创建新数组，只分配内存空间但不填充任何值 |
| eye、identity | 创建一个正方的N * N 单位矩阵（对角线为1，其余为0） |
In [116]:
data1 = (1,2,3,4)
In [117]:
np.asarray(data1)
Out[117]:
In [118]:
np.array(data1)
Out[118]:
In [119]:
data2 = ([2,2])
np.asarray(data2)
Out[119]:
In [120]:
import numpy as np
np.arange(15)
Out[120]:
In [121]:
ones
In [122]:
np.ones(19)
Out[122]:
In [123]:
np.zeros(10)
Out[123]:
In [124]:
np.empty(4)
Out[124]:
In [125]:
np.eye(3)
Out[125]:
In [126]:
np.eye(4)
Out[126]:
In [127]:
np.identity(2)
Out[127]:
In [128]:
np.identity(3)
Out[128]:
Recently I jsut moved from Shanghai to Kyoto, hence I have stopped taking notes for almost two weeks. From now on, I will continue writing this notes. Let's note~
YWFANG @Kyoto University November, 2017
dtype()
dtype 是一个特殊的对象，它含有ndarray将一块内存解释为特定数据类型的所需信息：
In [129]:
import numpy as np
arr1 = np.array([1,2,3], dtype = np.float64)
arr2 = np.array([1,2,3], dtype = np.int32)
In [130]:
arr1.dtype
Out[130]:
In [131]:
arr2.dtype
Out[131]:
dtype 是 NumPy 强大的原因之一。在多数情况下，它们直接映射到相应的机器表示，这使得“读写磁盘上的二进制数据流”以及“集成低级语言，如fortran"等工作变得简单。
下表记录了NumPy所支持的全部数据类型：（记不住没有关系，刚开始记不住也很正常）
| 类型 | 类型代码 | 说明 |
|---|---|---|
| int8、unit8 | i1、u1 | 有符号和无符号的8位（1个字节）整型 |
| int16、unit16 | i2、u2 | 有符号和无符号的16位（2字节）整型 |
| int32、unit32 | i4、u4 | 。。。32位。。。 |
| int64、unit64 | i8、u8 | 。。。64位。。。 |
| float16 | f2 | 半精度浮点数 |
| flaot32 | f4或者f | 标准单精度浮点数，与C的float兼容 |
| float64 | f8或d | 标准双精度浮点数，与C的double和Python的float对象兼容 |
| float128 | f16或者g | 扩展精度浮点数 |
| complex64、complex128 | c8、c16 | 分别用两个32位、64位或128位浮点数表示的复数 |
| complex256 | c32 | 复数 |
| bool | ？ | 存储True 或Flase 值的布尔类型 |
| object | O | Python多象类型 |
| string_ | S | 固定长度的字符串类型（每个字符1个字节）。例如，要创建一个长度位10的字符串，应使用S10 |
| unicode | U | 固定长度的unicode类型（字节数由平台决定）。跟字符串定义方式一样（如U10） |
我们可以通过 ndarray 的 astype 方法显示地转换其dtype：
In [132]:
import numpy as np
arr = np.array([1,2,3,4,5], dtype='i2')
print(arr.dtype)
print(arr)
In [133]:
float_arr = arr.astype(np.float64)
float_arr.dtype
Out[133]:
In the above example, an integer array was converted into a floating array.
In the following example, I will show you how to convert a float array to an int array. You will see that, if I cast some floating point numbers to be of interger type, the decimal part will be truncated.
In [134]:
import numpy as np
arr = np.array([1.2, 2.3, 4.5, 53.4,3.2,4.2])
print(arr.dtype)
print(arr)
print(id(arr)) #memoery address of arr
print('\n')
#conversion to integer
int_arr = arr.astype(np.int32)
print(int_arr.dtype)
print(int_arr)
If you have an array of strings representing numbers, you can also use 'astype' to convert them into numberic form:
In [135]:
import numpy as np
num_strings_arr = np.array(['1.25', '-9.6', '42'], dtype = np.string_)
print(num_strings_arr)
print(num_strings_arr.dtype)
float_arr = num_strings_arr.astype(np.float64)
# num_strings_arr.astype(float)
print(float_arr.dtype)
print(float_arr)
# alternatively, we can use a lazy writing
float1_arr = num_strings_arr.astype(float)
print(float_arr.dtype)
print(float_arr)
In addition, we can use another array’s dtype attribute:
In [136]:
# in this example, we can see that the int_arry will converted into
# a floating array, in particular, the dtype of calibers was used
# during the conversion using astype(calibers.dtype)
import numpy as np
int_array = np.arange(10)
print(int_array, int_array.dtype)
calibers = np.array([.22, .20, .23,.45, .44], dtype=np.float64)
print(calibers , calibers.dtype)
int_array_new = int_array.astype(calibers.dtype)
print(int_array_new, int_array_new.dtype)
In [137]:
#when stating an array, we can use the short code in the table to assign
# the dtype of the array
# for example
import numpy as np
empty_array = np.empty(8, dtype='u4')
print(empty_array)
print('\n')
zero_array = np.zeros(12, dtype='u4')
print(zero_array, zero_array.dtype)
print('\n')
one_array = np.ones(9, dtype='f8')
print(one_array, one_array.dtype)
print(*one_array)
点数（比如float64和float32）只能表示近似的分数值。因此复杂计算中，由于可能积累的浮点错误，比较浮点数字大小时，只能在一定的小数位数以内有效。
In [138]:
import numpy as np
arr = np.array([[1., 2., 3.,],[3.,5.,6.]])
print(arr.shape)
print(arr)
In [139]:
arr*arr
Out[139]:
In [140]:
arr-arr
Out[140]:
In [141]:
arr+arr
Out[141]:
同样地，当数组与标量进行算数运算时，也会遍历到各个元素
In [142]:
1/arr
Out[142]:
In [143]:
arr*2
Out[143]:
不同大小的数组之间的运算叫做广播 broadcasting，我们之后还会在第12章进行深度的学习。
In [144]:
import numpy as np
arr = np.arange(10, dtype='i1')
print(arr)
print(arr.dtype)
In [145]:
print(arr[0],arr[5])
print(arr[0:2])
In [146]:
arr[5:8]=12
print(arr)
In [147]:
#作为对比，我们回顾下之前列表的一些操作
list1=[0,1,2,3,4,5,6,7,8,9]
print(list1[:])
print(list1[0:2])
list1[5] = 12
print(list1[:])
list1[5:8]=12 #这里是跟数组很不同的地方
#如果不使用一个iterable，这里并无法赋值
print(list1[:])
如上面例子中看到的那种，当我们将标量赋值给一个切片时（arr[5:8]=12)，该值会自动传播（也就是12章将降到的broadcasting）到整个选区。跟列表最重要的区别在于，数组切片是原始数组的视图。这意味着数据不会被复制，视图上任何的修改都会直接反映到源数组上。
In [148]:
import numpy as np
arr = np.arange(10)
print(arr)
arr_slice = arr[5:8]
arr_slice[1] = 12345
print(arr)
arr_slice[:]=123
print(arr)
由于python常用来处理大数据，这种通过操作数组视图就可以改变源数组的方式，可以避免对数据的反复复制所带来的性能和内存问题。
如果我们想要得到的是一个数组切片的副本，而不是视图，就需要显式地进行复制操作，例如
In [149]:
import numpy as np
arr = np.arange(10)
arr_slice = arr[5:8]
arr_slice[1] = 12345
arr1 = arr[5:8]
print(arr1)
arr2 = arr[5:8].copy()
print(arr2)
#in this example，arr1仍然是数组的视图，
#但是arr2已经是通过复制得到的副本了
arr[5:8]=78
print('arr1 = ', arr1)
print('arr2 = ', arr2)
对于高维数组，能做的事情更多。在一个二维数组中，各个索引位置上的元素不再是标量，而是一维数组：
In [150]:
import numpy as np
arr2d = np.array([[1,2,3],[4,5,6],[7,8,9]])
arr2d[2]
Out[150]:
因此可以对各个元素进行递归的访问，不过这样需要做的事情有点多。我们可以传入一个以逗号隔开的索引列表来选区单个元素。也就是说，下面两种方式是等价的：
In [151]:
arr2d[0][2]
Out[151]:
In [152]:
arr2d[0,2]
Out[152]:
下图说明了二维数组的索引方式
在多维数组中，如果省略了后面的索引，则返回对象会是一个维度低一点的ndarray（它含有高一级维度上的所有数据）。
这里中文版的作者特别说明了上面这句话。括号外面的“维度”是一维、二维、三维之类的意思，而括号外面的应该理解为“轴”。也就是说，这里指的是“返回的低维度数组含有原始高维度数组某条轴上的所有数据。
下面看个例子来理解：
In [153]:
import numpy as np
arr3d = np.array([[[1,2,3],[4,5,6]],[[7,8,9],[10,11,12]]])
print(arr3d)
In [154]:
arr3d[0] #它是一个 2*3 数组
Out[154]:
标量值和数值都可以赋值给 arr3d[0]:
In [155]:
arr3d[0] = 42
print(arr3d)
In [156]:
print(arr3d[0,1])
print(arr3d[1,0])
注意，上面所有选取数组子集的例子中，返回的数组都是视图。
In [157]:
import numpy as np
arr = np.arange(10)
print(arr)
arr[4]=54
print(arr[1:6])
高维度对象的花样更多，我们可以在一个或者多个轴上进行切片、也可以跟整数索引混合使用。
In [158]:
import numpy as np
arr2d = np.array([[2,3,4],[3,5,5],[3,5,5]])
print(arr2d)
In [159]:
arr2d[:2]
Out[159]:
上述我们可以看出，这里的切片是沿着第0轴（即第一个轴）切片的。换句话说，切片是沿着一个轴向选取元素的。我们可以单次传入多个切片，就像传入多个索引那样：
In [160]:
arr2d[:2, :2]
Out[160]:
像上述这样的切片方式，只能得到相同维数的数组视图。我们还可以将整数索引与切片混合使用，从而得到低纬度的切片：
In [161]:
arr2d[2,:2]
Out[161]:
In [162]:
arr2d[:,:1] #这里，我们实现了对高维轴进行了切片
Out[162]:
自然地，对切片表达式的赋值操作也会被扩散到整个选区：
In [163]:
arr2d[:,:1] = 0
print(arr2d)
来看这样一个例子，假设我们有一个用于存储数据的数组以及一个存储姓名的数组（含有重复项）。在这里，我将使用 numpy.random 中的randn函数生成一些正态分布的随机数据。
In [164]:
%reset
In [165]:
import numpy as np
from numpy.random import randn
names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
#please make a comparison, if you use
# names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'], dtype='S4')
print(names, names.dtype)
type(names)
print('\n')
data = randn(7,4)
print(data, data.dtype, data.shape)
type(data)
Out[165]:
假设 names 数组中的每个名字都对应 data数组中的一行，而我们想要选出对应于名字‘Bob'的所有行。跟算数运算一样，数组的比较运算（如==）也是矢量化的。因此，对于names和字符串"Bob"的比较运算将会产生一个boolean array
In [166]:
names == 'Will'
Out[166]:
这个Boolean array可以用于数组索引，This boolean array can be passed when indexing the array:
In [167]:
data[names =='Will']
Out[167]:
当利用布尔型数组进行索引时候，必须注意布尔型数组的长度需要与被索引的轴长度一致。此外，还可以将布尔型数组跟切片、整数（或者整数序列，稍后对此进行详细的介绍）混合使用:
In [168]:
data[names =='Will', 2:]
Out[168]:
In [169]:
data[names =='Will', 2]
Out[169]:
要选择除了will以外的其他值，既可以使用不等于符号(!=)，也可以通过符号（-）对条件进行否定
In [170]:
names != 'Will'
Out[170]:
In [171]:
print(data[names != 'Will'])
In [172]:
data[-(names == 'Will')]
#this '-' was discarded in python3, alternatively we
# use '~'
In [173]:
data[~(names == 'Bob')]
# in python2, it should be
# data[-(names == 'Bob')]
Out[173]:
如果我们要选取这三个名字中的两个进行组合来应用多个布尔条件，需要使用&（和）、|（或）之类的布尔运算符：（注意，python关键字and和or在布尔型数组中是无效的）
In [174]:
mask = (names =='Bob') | (names == 'Will')
mask
Out[174]:
In [175]:
data[mask]
Out[175]:
值得注意的是，通过布尔索引选取数组中的数据，将总是创建数据的副本，即使返回一模一样的数组也是如此。
通过布尔型数组设置值是一种常用的方法。为了将data中的所有负数变为0，我们只需要
In [176]:
data[data<0] = 0
data
Out[176]:
通过一维布尔数组设置整行或列的值也很简单：
In [177]:
data[names != 'Will'] = 7
data
Out[177]:
In [178]:
#Suppose we had an 8 × 4 array:
import numpy as np
arr1 = np.zeros((8,4))
print(arr1)
print('\n')
for i in range(8):
arr1[i] = i+1
print(arr1)
为了以特定顺序选取行子集，只需传入一个用于指定顺序的整数列表或ndarray即可：
In [179]:
arr1[[4,3,0,6]]
Out[179]:
上面的代码的，我们用一个列表[4,3,0,6]就选出了arra1中的第4，3，0，6的子集。
如果我们使用负数进行索引，则选择的顺序将是从末尾到开头。
注意-0和0是一样的，还是开头的第一行作为0. 这是值得注意的地方。
In [180]:
arr1[[-4,-3,-1,-6,-0]]
Out[180]:
一次传入多个索引数组会会比较特别。它返回的是一个一维数组，其中的元素对应各个索引元组：
In [181]:
# 在第12章，我们会展开讲讲reshape，在这个例子中，我们只是使用 reshape
import numpy as np
arr = np.arange(32).reshape((8,4))
print(arr)
print('\n')
arr_select = arr[[1,5,7,2],[0,3,1,2]]
print(arr_select)
从上述代码的结果看不难看出，得出来的结果是[1,0] [5,3] [7,1] 和 [2,2]
那么怎么选取矩阵的行列子集呢？下面，我们只需要稍微改动下代码即可实现：（这部分最好再读几遍原书，字句不好理解）
In [182]:
import numpy as np
arr = np.arange(32).reshape((8,4))
print(arr)
print('\n')
arr_select = arr[[1,5,7,2]][:, [0,3,1,2]]
#1 5 7 2 选取行
#0 3 2 1 选取列
print(arr_select)
此外，还可以使用 np.ix_函数来实现上述的功能，它可以将两个一维整数数组转换为一个用于选取方形区域的索引器：
In [183]:
import numpy as np
arr = np.arange(32).reshape((8,4))
print(arr)
print('\n')
arr_select = arr[np.ix_([1,5,7,2],[0,3,1,2])]
print(arr_select)
It should be mentioned that, 花式索引与切片不一样，它总是将数据复制到新数组中。
In [184]:
import numpy as np
arr = np.arange(15).reshape((3,5))
print(arr)
In [185]:
print(arr.T)
print('\n')
print(arr)
当我们进行矩阵预算时候，进行需要用到转置操作。例如，要用 np.dot计算矩阵内积X$^T$X：
In [186]:
import numpy
from numpy.random import randn
arr = randn(6,3)
print(arr, '\n')
np.dot(arr.T, arr)
Out[186]:
对于更高维的数组，transpose 时需要得到一个由轴编号组成的元祖才能对这些轴进行转置（这个可能不好理解，得多阅读几次）：
In [187]:
#这里我简单举个例子
import numpy as np
arr = np.arange(16).reshape((2,2,4))
print(arr)
arr_transpose = arr.transpose((1))
从上面几个例子，我们可以看出，对于简单的低维矩阵，使用.T就可以实现转置，毕竟只是进行轴对换而已；但是对于高维数组，就显得麻烦好多。ndarray还有一个swapaxes方法，它需要接受一对轴编号：(注意swapaxes也是返回源数据的视图，并不会进行任何复制操作。)
In [188]:
import numpy as np
arr = np.arange(18).reshape(3,3,2)
print(arr, '\n')
arr_axes1 = arr.swapaxes(0,1)
print(arr_axes1)
print('\n')
arr_axes2 = arr.swapaxes(1,2)
print(arr_axes2)
In [ ]:
In [189]:
import numpy as np
arr = np.arange(10)
In [190]:
print(arr, '\n')
print(np.sqrt(arr))
In [191]:
print(arr,'\n')
np.exp(arr) #the results are e^N (N = 0, 1, 2,...)
Out[191]:
上述这些都是一元（unary）ufunc。另外一些（如add或maximum）接受2个数组（因此也叫二元binary ufunc），并返回一个结果数组：
In [192]:
import numpy as np
from numpy.random import randn
x = randn(8)
print(x,'\n')
y = randn(8)
print(y,'\n')
max_number = np.maximum(x,y)
print(max_number,'\n')
此外，有一小部分的ufunc，它们可以返回多个数组。mof就是一个例子，它是Python内置函数 divmod的矢量化版本，用于分离浮点数组的小数和整数部分。通过下面的例子，我们会发现，mof其实得到的是几个数组组成的tuple
In [193]:
import numpy as np
from numpy.random import randn
arr = randn(7)*5
print(arr,'\n')
arr_1 = np.modf(arr)
print(arr_1)
print(type(arr_1))
print(arr_1[1])
下表中列出了一些一元和二元ufunc
| 函数 | 说明 |
|---|---|
| abs, fabs | 计算整数、浮点数和负数的绝对值。对于复数值，可以使用更快的fabs |
| sqrt | 计算各元素的平方根。相当于 arr ** 0.5 |
| square | 计算各元素的平方。相当于是 arr ** 2 |
| exp | 计算各元素的e指数，即 e$^x$ |
| log,log10,log2,log1p | 分别对应自然对数（以e为底），底数是10的log，底数是2的log，以及log(1+x) |
| sign | 计算各元素的正负号：1代表整数，0代表零，-1代表负数 |
| ceil | 计算各元素的ceiling值，即大于等于该值的最小整数 |
| floor | 计算各元素的floor值，即小于等于该值的最大整数 |
| rint | 将各元素之四舍五入到最接近的整数，保留dtype |
| modf | 将数组的小数和整数部分以两个独立的数组形式返回 |
| isnan | 返回一个表示“哪些值是NaN（这不是一个数字）”的boolean数组 |
| isfinite、isinf | 分别返回一个表示“哪些元素是有穷的（非inf，非NaN）” 或者 “哪些元素是无穷的”的布尔型数组 |
| cos、cosh、sin、sinh、tan，tanh | 普通型和双曲型三角函数 |
| arccos、arccosh、arcsin、arcsinh、arctan、arctanh | 反三角函数 |
| logical_not | 计算各个元素not x的真值。相当于-arr |
| 函数 | 说明 | |
|---|---|---|
| add | 将数组中对应的元素相加 | |
| substract | 从第一个数组中减去第二个数组中的元素 | |
| multiply | 数组元素相乘 | |
| divide、floor_divide | 除法或向下圆整除法（丢弃余数） | |
| power | 对第一个数组中的元素A，根据第二个数组中的相应好元素B，计算A$^B$ | |
| maximum, fmax | 元素级的最大值计算。fmax将忽略NaN | |
| minimum、fmin | 元素级的最小值计算。fmin将忽略NaN | |
| mod | 元素级的求模计算（除法的余数） | |
| copysign | 将第二个数组中的值的符号复制给第一个数组中的值 | |
| greater、greater_equal、less、less_equal、equal、not_equal | 执行元素级的比较运算，最终产生boolean型数组。相当于中缀运算>, >=, <, <=, ==, != | |
| logical_and、logical_or、logical_xor | 执行元素级的真值逻辑运算。相当于中缀运算符 '&'，'$ | $'，'^' |
In [194]:
import numpy as np
from numpy.random import randn
new = randn(10)
In [195]:
new
Out[195]:
In [196]:
np.sign(new)
Out[196]:
In [ ]:
In [197]:
import numpy as np
new = np.arange(10)
new = new+0.1
print(new,'\n')
np.ceil(new)
Out[197]:
In [ ]:
In [198]:
import numpy as np
from numpy.random import randn
new = randn(10)
print(new,'\n')
print('rint function:', np.rint(new))
print('isnan function: ', np.isnan(new))
print('isfinite function', np.isfinite(new))
print('isinf function: ', np.isinf(new))
print('logical_not function: ', np.logical_not(new))
In [199]:
#Revieing some knowledge I have learnt
import numpy as np
arr1 = np.arange(16,dtype='i4').reshape(2,2,4)
arr2 = np.arange(10,dtype='float')
print(arr1)
print('\n')
print(arr2)
print('\n')
arr3=arr1.copy()
arr3[1]=23
print(arr3.astype(arr2.dtype))
In [200]:
sum(arr1,arr3)
Out[200]:
In [201]:
print('mean value = ', arr1.mean(), '\n' 'max value is ',
arr1.max(), '\n' 'std root = ', arr1.std(), '\n'
'The sum of all the elements = ', arr1.cumsum(),
'\n' 'The multipy of all the elements = ', arr1.cumprod())
In [ ]:
NumPy数组的矢量化在很大程度上简化了数据处理方式。一般而言，矢量化运算要比等价的纯python方式快1-2个数量级，尤其是在数值计算处理过程中这个优势更加的明显。在后面的第12章节中，我们将了解到广播，它是一种针对矢量化计算的强大手段。
假设我们想要在一组值（网格型）上计算sqrt(x^2+y^2)。我们当然可以选择用loop的方式来计算，但是我们在这里使用数组的方法。
np.meshgrid 函数接受两个一维数组，并产生两个二维矩阵（对英语两个数组中所有的(x,y)对）：
In [202]:
import numpy as np
points = np.arange(-1,1,0.5) # 产生4个间隔同为0.5的点。
print(points[:10],'\n')
xs, ys = np.meshgrid(points, points)
print('xs is \n',xs,'\n')
print('transposed xs is \n', xs.T)
print('ys is \n', ys, '\n')
现在，我们来计算xs二次方与ys二次方的和：
In [203]:
z = np.sqrt(xs**2 + ys**2)
print('z = \n', z)
我们试着将上述这个z函数画出来
In [204]:
%matplotlib inline
import matplotlib.pyplot as plt
#Here, the matplotlib function 'imshow' was used
# to create an image plot from a 2D array of function values
plt.imshow(z, cmap=plt.cm.gray);
plt.colorbar()
Out[204]:
上面，我们只使用了很“疏”的点，接下来，我们尝试使用很密集的点，这样有利于我们可视化sqrt(x^2+y^2)这个函数。
In [205]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
arr1=np.arange(-1,1,0.001)
print(arr1)
xs1,ys1=np.meshgrid(arr1,arr1)
#print(xs)
z1 = np.sqrt(xs1**2+ys1**2)
print(z1)
plt.imshow(z1, cmap=plt.cm.gray)
plt.colorbar()
Out[205]:
In [206]:
import numpy as np
xarr = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
yarr = np.array([2.1, 2.2, 2.3, 2.4, 2.5])
cond = np.array([True, False, True, True, False])
假设我们想要根据 cond 中的值来决定我们是选取 xarr 还是 yarr 的值。当 cond 中的值为 True 时，我们选取 xarr 中的值，否则选用 yarr 中的数值。
python中列表推导式的写法如下所示：
In [215]:
result = [(x if c else y)
for x, y, c in zip(xarr, yarr, cond)]
print(result)
It has multiple problems here. First, it will not be fast for large arrages (because all the work is being done in interpreted python code,即纯python处理)；second, it will not work with multidimensional array,即无法处理多维数组。
如果我们使用 np.where，we can wirte this code very concisely：
In [213]:
result_where = np.where(cond, xarr, yarr)
print(result_where)
np.where的第二个和第三个参数不必是数组，它们可以是标量。在数据分析工作中，where 通常用于根据另一个数组而产生一个新的数组。假设有一个由随机数据组成的矩阵，我们想将所有正的值替换为2，所有负值改为-2。那么我们可以写为：
In [221]:
from numpy.random import randn
import numpy as np
arr_a = randn(10)
print(arr_a)
arr_b = np.where(arr_a <0, -2, 2)
print(arr_a)
print(arr_b)
如果我们只需要把负的值改为 -3， 那么我们可以用
In [223]:
arr_c = np.where(arr_a < 0, -3, arr_a)
print(arr_c)
Highlight： 我们可以使用where表现更加复杂的逻辑。想象这样一个例子：有两个boolean array，分别叫做cond1和conda2，希望使用四种不同的布尔值组合实现不同的赋值操作.
如果我们不用where，那么这个pseudo code 的逻辑大概如下
虽然不是那么容易看出来，我们可以使用 where 的嵌套来实现上述的pseudocode逻辑
np.where(conda1 & conda2, 0, np.where(conda1, 1, np.where(conda2, 2, 3)))
在这个特殊的例子中，我们还可以利用“布尔值在计算过程中被当作0或者1处理”这个事实，将上述result的结果改写成
In [ ]:
result = 1*（cond1 - cond2) + 2 *(cond2 & -cond1) + 3*-(cond1|cond2)
#这种写法我觉得并不是太推荐，在2017年的新版中，原作者写删除了这部分的讨论
现在我们来应用下上面的嵌套np.where
In [13]:
import numpy as np
from numpy.random import randn
x1 = randn(10)
y1 = randn(10)
cond1 = np.where(x1<0, True, False)
cond2 = np.where(y1>0, True, False)
result=np.where(cond1 & cond2, 0,
np.where(cond1, 1,
np.where(cond2, 2, 3)))
print(result)
In [22]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
arr = np.random.randn(5,4)
print(arr)
plt.imshow(arr, cmap=plt.cm.gray)
plt.colorbar()
Out[22]:
In [30]:
arr[0,2]
Out[30]:
In [31]:
arr[0,3]
Out[31]:
上面代码中，我产生了一些 enormally distributed random data，并且用imshow function 把这个二维数组给画了出来。我们可以使用 aggregate statistics 做一些计算. （其实我在前面已经用到过了这些 array 实例方法。
In [23]:
arr.mean() #这里的使用方法就是作为 array instance method
Out[23]:
In [24]:
np.mean(arr) # 这里的使用方法就是 top-level Numpy function
Out[24]:
In [25]:
arr.sum()
Out[25]:
mean 和 sum 这类的函数可以接受一个 axis 参数 （用于计算该轴向上的统计值），最终结果是一个相对于原数组少了一维的数组：
In [28]:
arr.mean(axis=1) # compute mean across the columns
Out[28]:
In [32]:
arr.mean(axis = 0 ) # compute mean down the rows
#如果对axis = 0 和 axis = 1 两个结果的列数有不解，
# 可以回顾一下前面的二维数组的索引方式
# 特别是那个 NumPy 数组元素索引的图，那个图上描述了 axis1 和 axis0 的相对朝向。
Out[32]:
其他如 ‘cumsum’， ‘cumprod’ 这类函数方法并不聚合，而是产生一个由中间结果组成的数组：
English: Other methods like cumsum and cumprod donot aggregate, instead producing an array of the intermediate results:
In [39]:
arr1 = np.array([0,1,2,3,4,5,6,7])
arr1.cumsum()
Out[39]:
In multidimensional arrays, accumulation functions like cumsum return an array of the same size, but with the partial aggregates computed along the indicated axis according to each lower dimensional slice:
In [36]:
arr2 = np.array([[0,1,2],[3,4,5],[6,7,8]])
In [37]:
arr2.cumsum(0)
Out[37]:
在上述方法中，布尔值会被强制转换为 1 （True） 和 0 （False）。因此，sum 经常被用来对Boolean数组中的True值计算：
In [2]:
import numpy as np
from numpy.random import randn
arr = randn(100)
print(arr)
(arr>0).sum() # 正值的个数
Out[2]:
另外还有两个方法 any 和 all，它们对 Boolean array 很有用。any用于测试数组中是否存在一个或多个True，而all则检查数组中所有值是否都是True。
In [3]:
bools = np.array([False, False, True, False])
In [4]:
bools.any() # any of them if True, then the return result is True
Out[4]:
In [5]:
bools.all() # all of them should be True, otherwise the return result is False
Out[5]:
跟Python内置的列表一样，NumPy 数组也可以通过 sort 方法就地排序
In [8]:
import numpy as np
from numpy.random import randn
arr_a = randn(8)
print(arr_a)
In [10]:
arr_a.sort() # 注意，它将直接改变数组本身
print(arr_a)
对于多维数组，只要我们给定确定的轴编号，它就会沿着特定的轴进行排序。我们这里拿一个二维数组举例
In [43]:
import numpy as np
arr_b = randn(4,5)
print(arr_b)
arr_c = arr_b.copy()
print('\n')
print(arr_c)
In [44]:
arr_b.sort(1) # 这里我们沿着 axis = 1 方向进行排序，我们发现每个一位数组中的元素都被排序了
print(arr_b)
In [45]:
arr_c[2].sort() #这里我们只选择了编号为2的那个一维数组进行排序
print(arr_c)
The top-level method 'np.sort' returns a sorted copy of an array instead of modifying the array in-place. 这个需要我们区分 np.sort 和数组实例 sort 的地方。
In [48]:
np.sort(arr_c)
print(arr_c, '\n')
print(np.sort(arr_c))
数组 sort 的应用之一，就是确定数组的分位数(quantile)。
A quick-and-dirty way to compute the quantiles of an array is to sort it, and select the value at a particular rank.
In [57]:
short_arr = randn(10)
short_arr.sort()
print(short_arr, '\n', len(short_arr))
In [52]:
short_arr[int(0.1*len(large_arr))] #(处于0.1分位数位置上)
Out[52]:
In [53]:
short_arr[int(0*len(large_arr))] #（处于0分位上）
Out[53]:
上面我们只是使用了很小的数组，我们一眼就可以看出各分位数上的数值；当数组变得很大时候，才能凸显出 sort 的便捷。例如：
In [56]:
large = randn(2000)
large.sort()
large[int(0.5*len(large))]
Out[56]:
关于 NumPy 排序方法以及诸如间接排序之类的高级技术，我们在第12章还会详细的讨论，在 Pandas 中也有一些特别的排数技术。
In [58]:
import numpy as np
names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
np.unique(names)
Out[58]:
In [59]:
ints = np.array([3,3,41,4424,523,523,22,22,43]
)
np.unique(ints)
Out[59]:
我们可以拿着与 np.unique 等价的纯python代码来比较一下（Contrast no.unique with the pure Python alternative:)
In [60]:
sorted(set(names))
Out[60]:
In [61]:
sorted(set(ints))
Out[61]:
Anotehr function, np.in1d, tests membership of the values in one array in another, returning a boolean array.
另一个函数np.in1d用于测试一个数组的值在另一个数组中的成员资格，返回一个Boolean array
In [64]:
values = np.array([6,623,43,22,3])
np.in1d(values,[6,43,22])
Out[64]:
In [67]:
np.unique(values)
Out[67]:
In [68]:
np.intersect1d([3,6],[3,22,43])
Out[68]:
In [70]:
np.union1d([3,6],[3,22,43])
Out[70]:
In [71]:
np.in1d([3,6],[3,22,43])
Out[71]:
In [72]:
np.in1d([3,6],[3,6,22])
Out[72]:
In [79]:
np.setdiff1d([3,22,6],[6])
Out[79]:
In [81]:
np.setxor1d([3,22,6],[6])
Out[81]:
In [82]:
import numpy as np
arr_c = np.arange(10)
np.save('./chapter04/some_array',arr_c)
ru如果文件路径末尾没有扩展名 .npy，那么这个扩展名会被自动补全。然后就可以通过 np.load 读取磁盘上的数组：
In [86]:
np.load('./chapter04/some_array.npy') # 注意，需要指明文件后缀名。
Out[86]:
通过 np.savez 可以将多个数组保存到一个uncompressed npz文件中（注意原书和中文翻译的第一版都把这个npz说成了是压缩文件，这个是错误的，但是原作者第二版，即利用python 3的版本已经更正了，我也查阅了 NumPy 的文档，np.savez保存的并不是压缩文件，如果要压缩文件，可以使用 np.savez_compressed），将数组以关键字参数的形式传入即可：
In [87]:
np.savez('./chapter04/array_archive.npz', a = arr_c, b = arr_c)
In [101]:
np.savez_compressed('./chapter04/array_compressed.npz', a1 = arr_c,
b1 = arr_c)
When loading an .npz file, we get back a dict-like object （我们得到的是一个类似字典的对象）that laods the individual arrays lazily (该对象会对各个数组进行延迟加载)
In [88]:
arch = np.load('./chapter04/array_archive.npz')
In [89]:
arch['b']
Out[89]:
In [104]:
arch = np.load('./chapter04/array_compressed.npz')
In [105]:
arch['b1']
Out[105]:
In [90]:
!cat ./chapter04/array_ex.txt #for windows system, use !type
该文件可以被加载到一个二维数组中，如下所示：
In [93]:
arr = np.loadtxt('chapter04/array_ex.txt', delimiter = ',')
arr
Out[93]:
In [94]:
print(arr)
np.savetxt 执行的是相反的操作：将数组写到以某种分隔符分开的文本文件中去。 genfromtxt 跟 loadtxt 差不多，只不过它面向的是结构化数组和缺失数据处理。在12章中，我们还会仔细讨论结构化数组的知识。
In [96]:
np.savetxt('./chapter04/array_ex-savetxt.txt', arr)
In [97]:
!cat chapter04/array_ex-savetxt.txt
Linear algebra, like matrix multiplication, decompisitions, determinants, and other square matrix math, is an important part of any array library. Unlike MATLAB, multiplying two two-dimensional arrays with * is an element-wise product instead a matrix dot product. Thus, there is a function 'dot', both an array method and a function in the numpy namespace, for matrix multiplication:
In [106]:
x = np.array([[1., 2., 3.],[4., 5., 6.]])
y = np.array([[6.,23.,], [-1, 7], [8, 9]])
In [107]:
x
Out[107]:
In [108]:
y
Out[108]:
In [109]:
x.dot(y)
Out[109]:
x.dot(y) is equivalent to np.dot(x,y)
In [111]:
import numpy as np
np.dot(x,y)
Out[111]:
A matrix product between a 2D array and a suitably sized 1D array result in a 1D array:
In [114]:
np.dot(x, np.ones(3))
Out[114]:
In [6]:
import numpy as np
x1 = np.array([[2,2],[3,3]])
y1 = np.array([[1,1],[1,1]])
dotvalue1=x1.dot(y1)
dotvalue2=np.dot(x1,y1)
print('dotvalue1 = \n', dotvalue1, '\n' 'dotvalue2 = \n', dotvalue2)
numpy.linalg 中有一组标准的矩阵分解运算以及诸如求逆行列式之类的东西。它们跟 matlab 和 R 等语言所使用的是相同的行业标准级 Fortran 库，如 BLAS、LAPACK、Intel MKL （可能有，这个取决于所使用的 NumPy 版本）等：
In [9]:
from numpy.linalg import inv, qr
from numpy.random import randn
In [14]:
X = randn(5,5)
In [15]:
mat = X.T.dot(X)
In [16]:
inv(mat)
Out[16]:
In [17]:
mat.dot(inv(mat))
Out[17]:
In [23]:
q, r = qr(mat)
# QR decomposition is decomposition of a matrix A into a product
# A = Q R
# where Q is an orthogonal matrix, and R is upper triangular matrix
In [24]:
q
Out[24]:
In [26]:
q.dot(q.T) # we can see that the result is an unit matrix
Out[26]:
In [27]:
r
Out[27]:
这里，我对常用 numpy.linalg 函数进行一些案例对说明，原书利用了一个表格，但是我自己为了看这本书也重复写了几个表格了，记忆情况并不佳，可能还是一个函数一个例子对这种方法更加容易让人记忆深刻一些。（2017/12/25）
| 函数 | 说明 |
|---|---|
| diag | 以一维数组的形式返回方阵的对角线（或非对角线）元素 |
In [38]:
import numpy as np
from numpy.random import randn
arr_d = np.arange(10)
arr_e= randn(16).reshape(4,4)
print(arr_d, '\n', arr_e)
In [32]:
np.diag(arr_d) # 对于1D array，它将转化为一个对角化的矩阵
Out[32]:
In [41]:
np.diag(arr_e) #当diag作用在一个2D以上数组时，则返回对角线上的元素。
Out[41]:
| 函数 | 说明 |
|---|---|
| dot | matrix multiplication, 矩阵乘法 |
这已经在前面举过例子，这里略了。
| 函数 | 说明 |
|---|---|
| trace | 计算对角线元素的和 |
In [51]:
np.trace(arr_e)
Out[51]:
In [55]:
np.trace(np.diag(arr_d))
Out[55]:
In [72]:
arr_f = np.arange(64).reshape(4,4,4)
print(arr_f
)
In [82]:
np.trace(arr_f[1])
#可以看到对于多维数组，我们还可以对其中低维度的求trace
Out[82]:
| 函数 | 说明 |
|---|---|
| det | 计算矩阵的行列式 |
In [58]:
from numpy.linalg import det
det(arr_e)
Out[58]:
| 函数 | 说明 |
|---|---|
| eig | 计算方阵的本征值和本征向量 |
In [90]:
import numpy as np
from numpy.linalg import eig
eig(arr_e)
# computes the eigenvalues and eigenvectors of a square matrix
#关于这个函数的数学含义，请参考线性代数相关的书籍
Out[90]:
| 函数 | 说明 |
|---|---|
| inv | 计算方阵的逆 |
In [111]:
from numpy.linalg import inv
# if 'a' is a matrix object,
# the return value is a matrix as well
a = np.array([[1., 2.], [3., 4.]])
ainv = inv(a)
print(a, '\n', ainv)
In [123]:
# inverses of several matrices can be computed at once:
b = np.array([
[
[1.,2.],[3., 4.]
],
[
[1., 3.],[3., 5.]
]
])
binv = inv(b)
binv
Out[123]:
| 函数 | 说明 |
|---|---|
| pinv | 计算矩阵的Moore-Penrose伪逆 |
Compute the Moore-Penrose pseudo-inverse of a matrix: The pseudo-inverse of a matrix $A^+$, is defined as "the matrix that 'sloves' [the least-squares problem ] Ax = b," i.e., if $\bar{x}$ is said solution, then $A^+$ is that matrix such that $\bar{x}$ = $A^+$b
For more information, please refere to linear algebra books
In [128]:
from numpy.linalg import pinv
from numpy.random import randn
c = randn(9,6)
bpinv = pinv(c)
bpinv
Out[128]:
| 函数 | 说明 |
|---|---|
| qr | copute the QR decompisition |
上面提过了，此处略
| 函数 | 说明 |
|---|---|
| svd | 计算奇异值分解， compute the singular value decomposition SVD |
In [134]:
import numpy as np
from numpy.random import randn
a = randn(9,6) + 1j*randn(9,6)
a
Out[134]:
In [135]:
# Reconstruction based on full SVD
# factors the matrix a as u * np.diag(s) * V,
# where u and v are unitary and s is a 1D array of a's
# singular values
U, s, V = np.linalg.svd(a, full_matrices = True)
U.shape, s.shape, V.shape
Out[135]:
In [139]:
s
Out[139]:
| 函数 | 说明 |
|---|---|
| solve | 计算方程组 Ax = b， 其中 A 为一个方阵 |
Note：The solutions are computed using LAPACK routine_gesv. 'a' must be square and of full-rank, i.e., all rows (or, equivalently, columns) must be linearly independent; if either is not true, use lstsq for the least-squares best 'solutions' of the system/equation.
In [141]:
# Solve the systems of equatons 3* x0 + x1 = 9 and x0 + 2*x1 =8
from numpy.linalg import solve
a = np.array([[3,1],[1,2]])
b = np.array([9,8])
x = solve(a,b)
In [142]:
a
Out[142]:
In [143]:
b
Out[143]:
In [144]:
x
Out[144]:
In [145]:
# check that the solution is correct:
np.allclose(np.dot(a, x), b)
Out[145]:
| 函数 | 说明 |
|---|---|
| lstsq | 计算方程组 $ax = b$ 的最小二乘解 |
numpy.linalg.lstsq(a, b, rcond=-1)
Return the equation $a$ $x$ = $b$ by computing a vector $x$ that minimizes the Eouliden 2-norm $ ||b - ax ||^2$. The equation may be under-, over-, or well-determined (i.e. the numer of linearly independent rows of $a$ can be less less, greater, or less than its number of linearly independent columns). If $a$ is square and of full rank, then $x$ (but for round-off error) is the "exact" solution of the equation. (reivsed on the content from scipy.org)
In [147]:
x = np.array([0,1,2,3])
y = np.array([-1, 0.2, 0.9, 2.1])
通过查看上面x和y的点，我们可以发现这个曲线的大概斜率在1左右，而在y轴上的cut off在-1左右。
我们可以重新写一下上面这个线性方程：$y$ = A p, 此处 A = [ [x, 1] ] 并且 p = [[m], [c]]。现在我们使用 lstsq 去解 p
In [152]:
A = np.vstack([x, np.ones(len(x))]).T
# np.vstack 可以按顺序把array堆叠在一起
# 此处是把 x 和 np.ones(4) 堆在了一起
In [153]:
A #即 A = [[x, 1]]
Out[153]:
In [155]:
from numpy.linalg import lstsq
m, c = lstsq(A, y)[0]
print(m, c)
In [173]:
# 我们把数据和拟合的线可以画出来
####basic settings started
import matplotlib.style
import matplotlib as mpl
mpl.style.use('classic')
mpl.rcParams['figure.facecolor'] = '1'
#if choose the grey backgroud, use 0.75
mpl.rcParams['figure.figsize'] = [6.4,4.8]
mpl.rcParams['lines.linewidth'] = 1.5
mpl.rcParams['legend.fancybox'] = True
#####basic settings finished
%matplotlib inline
# plot inline jupyter
import matplotlib.pyplot as plt
# plot orginal data (i.e. four points)
plt.plot(x, y, 'o', label = 'Original data', ms =14)
# plot the fitted line using red line style and linewidth = 2
plt.plot(x, m*x + c, 'r', lw=2, label = 'Fitted line')
# plot the legend
plt.legend()
# plot grid
plt.grid()
plt.show()
In [194]:
#因为上面用到了 numpy.stack,那么我就顺便再举一个例子来说明 vstack 的用法
# 与 vstack 相反的操作是 vsplit
import numpy as np
a = np.array([1,2,3])
b = np.array([4,5,6])
ab = np.vstack((a, b))
m, n= np.vsplit(ab, 2) # 把ab分成2个，分别存在在m和n中
print(ab)
print(m)
numpy.random 模块对 PYthon 内置对 random 进行了补充，增加了一些搞笑生成随机样本对函数。例如，我们可以用normal来得到一个标准正态分布对 4 * 4 样本数组：
In [9]:
import numpy as np
samples = np.random.normal(size=(4,4))
In [10]:
samples
Out[10]:
与此对比地，在python对内置random函数中，一次只能生成一个样本值。下面我们就来对比下这两种方法对区别，我们将会看到 numpy中对模块有更优越对效率：
In [11]:
from random import normalvariate
N = 1000000
In [12]:
%timeit samples = [normalvariate(0,1) for _ in range(N)]
#中文翻译版本中这行代码是错的，翻译者写成了 xrange
In [14]:
%timeit np.random.normal(size=N)
下表列出了 numpy.random 中的部分函数。在下一节中，我门将给出一些利用这些函数一次性生成大量样本值的案例。
| 函数 | 说明 |
|---|---|
| seed | 确定随机数生成器的种子 |
| permutation | 返回一个序列的随机排列或返回一个随机排列的范围 |
| shuffle | 对一个序列就地随机排列 |
| rand | 产生均匀分布的样本值 |
| randint | 从给定的上下限范围内随机选取整数 |
| randn | 产生正态分布（平均值为0，标准差为1）的样本值，类似于matlab接口 |
| binomial | 产生二项分布的样本值 |
| normal | 产生正态（高斯）分布的样本值 |
| beta | 产生Beta分布的样本值 |
| chisquare | 产生卡方分布的样本值 |
| gamma | 产生Gamma分布的样本值 |
| uniform | 产生在[0,1)中均匀分布的样本值 |
In [41]:
####basic settings started
import matplotlib.style
import matplotlib as mpl
mpl.style.use('classic')
mpl.rcParams['figure.facecolor'] = '1'
#if choose the grey backgroud, use 0.75
mpl.rcParams['figure.figsize'] = [6.4,4.8]
mpl.rcParams['lines.linewidth'] = 1.5
mpl.rcParams['legend.fancybox'] = True
#####basic settings finished
%matplotlib inline
# plot inline jupyter
import matplotlib.pyplot as plt
import numpy as np
position = 0 # 初始化位置
walk = []
steps = 1000
for _ in range(steps):
stepwidth = 1 if np.random.randint(0,2) else -1
position += stepwidth
walk.append(position)
#print(walk)
#plot this trajectory
plt.plot(walk[:1000])
Out[41]:
注意我上面的代码跟原书上的区别，主要在于我并不是python自身的random standard library。我使用的是numpy.random，这两个是有区别的，这主要在于
numpy.random.randint(a, b)，返回的值是a ~ (b-1)之间的整数值（包括a 和 b-1）；
而python自带的random.randint(a,b) 返回的值是 a ～ b之间的整数值（包括a和b）
In [42]:
import random
position = 0 # 初始化位置
walk = []
steps = 1000
for _ in range(steps):
stepwidth = 1 if random.randint(0,1) else -1
position += stepwidth
walk.append(position)
#print(walk)
#plot this trajectory
plt.plot(walk[:1000])
Out[42]:
上面的walk数值，其实就是随机数的累计和。不过上面的方式中，我门都是走一步然后产生一个随机数，其实我们可以用numpy.random.randint一次性地产生N个随机数，这里以N=1000为例
In [1]:
nsteps = 1000
In [2]:
import numpy as np
draws = np.random.randint(0, 2, size = nsteps)
In [4]:
steps = np.where(draws > 0, 1, -1)
In [5]:
walk = steps.cumsum()
In [6]:
walk.min()
Out[6]:
In [7]:
walk.max()
Out[7]:
A more complicated statics is the 'first crossing time', the step at which the random walk reaches a particular value. Here we might want to know how long it look the random walk to get at least 10 steps aways from the origin 0 in either direction. np.abs(walk) >= 10 gives us a boolean array indicating where the walk has reached or exceeded 10, but we want the index of the first 10 or -10.
In [9]:
(np.abs(walk) >= 10).argmax()
#Note that using argmax here is not always efficient because
#it always makes a full scan of the array. In this special case,
#once a True is observed we know it to be the maxi‐ mum value.
Out[9]:
如果希望模拟多个随机漫步过程，只需要对上面对代码做一点微调。我们只需要给numpy.random 传入一个二元元祖即可产生一个二维数组，然后我们就能一次性计算5000个随机漫步过程（一行一个）的累计和了。
In [1]:
import numpy as np
In [2]:
nwalks = 5000
In [3]:
nsteps = 1000
In [4]:
draws = np.random.randint(0, 2, size=(nwalks, nsteps)) # 0 or 1
In [6]:
steps = np.where(draws > 0, 1, -1)
In [12]:
walks = steps.cumsum(1)
walks
Out[12]:
In [13]:
walks.max()
Out[13]:
In [14]:
walks.min()
Out[14]:
得到这些数据后，我们可以来计算出30或者-30的最小穿越时间。这里得要稍微动一下脑子，因为不是5000个过程都到达了30。我们可以用any方法来对此进行检查
In [15]:
hits30 = (np.abs(walks) >= 30).any(1)
In [16]:
hits30
Out[16]:
In [18]:
hits30.sum() #到达30或者-30的数量
Out[18]:
然后我们再利用这个boolean array选出哪些穿越了30（绝对值）的随机漫步（行），并调用argmax在轴1上获取穿越时间：
In [20]:
crossing_times= (np.abs(walks[hits30]) >= 30).argmax(1)
In [21]:
crossing_times.mean()
Out[21]:
这里请尝试其他分布方式得到漫步数据。只需要使用不同的随机数生成函数即可。例如，normal 用于生成指定均值和标准差的正态分布数据
In [23]:
steps = np.random.normal(loc=0, scale=0.25,
size=(nwalks,nsteps))
In [24]:
steps
Out[24]:
In [ ]:
In [ ]:
In [1]:
#python list
x = [1,2,3,4]
y = [5,6,7,8]
In [2]:
x*2
Out[2]:
In [3]:
x+10
In [4]:
x+y
Out[4]:
In [5]:
#numpy arrays
import numpy as np
ax = np.array([1,2,3,4])
ay = np.array([5,6,7,8])
In [6]:
ax*2
Out[6]:
In [7]:
ax+10
Out[7]:
In [8]:
ax+ay
Out[8]:
从上面可以看出，numpy数组操作时候是对被作用对所有元素的，这一事实使得数组计算都变得简单和快速。比如我们可以快速地计算多项式：
In [9]:
def f(x):
return 3*x**2 + 2*x +7
f(ax)
Out[9]:
NumPy提供了一些通用函数集合，他们也能对数组进行直接对操作。这些通用函数可以作为math模块中所对应函数对替代。示例如下：
In [10]:
np.sqrt(ax)
Out[10]:
In [11]:
np.cos(ax)
Out[11]:
使用NumPy中的通用函数，其效率要比对数组进行迭代然后使用math模块中的函数每次只处理一个元素快上数倍。因此，只要有可能就应该直接使用这些通用函数。
在底层，NumPy 数组的内层分配方式和C和Fortran是一样的。他们在大块的连续内存中存储。正因如此，NumPy才能创建比通常Python列表大许多的数组。例如，如果像创建10000 * 10000的二维浮点函数，这对numpy而言是很轻松的事情：
In [12]:
grid = np.zeros(shape=(10000,10000), dtype=float)
In [13]:
grid
Out[13]:
所有的通用操作仍然可以同时施加于所有的元素之上：
In [14]:
grid+10
Out[14]:
In [16]:
np.sin(grid+10)
Out[16]:
关于NumPy，一个特别值得提起的方面就是NumPy扩展了python列表的索引功能——尤其是针对多维数组时更是如此。现在我们来构建一个简单的二维数组然后做一些简单的experiment
In [60]:
import numpy as np
x=list(range(1,5))
y=list(range(5,9))
z=list(range(9,13))
a = np.array(x)
b = np.array(y)
c = np.array(z)
array1 = np.concatenate((a, b, c), axis=0)
array2 = np.stack((a, b, c), axis=0)
In [54]:
array1
Out[54]:
In [61]:
array2
Out[61]:
In [63]:
#select row 1
array2[1]
Out[63]:
In [67]:
#select column 1
array2[:,1]
Out[67]:
In [68]:
array2[1:3,1:3]
Out[68]:
In [69]:
array2[1:3,1:3] += 10
In [70]:
array2
Out[70]:
In [71]:
#broadcast a row vector across an operation on all rows
array2+[100,101,102,103]
Out[71]:
In [72]:
array2
Out[72]:
In [73]:
#conditional assigan on an array
np.where(a < 10, a, 10) #a<10 is the condition, if ture, return a. I have introduced np.where before in this chapter
Out[73]:
In [74]:
import numpy as np
m = np.matrix([[1,-2,3],[0,4,5],[7,8,-9]])
m
Out[74]:
In [75]:
#Return transpose 转置矩阵
m.T
Out[75]:
In [76]:
#return inverse 逆矩阵
m.I
Out[76]:
In [77]:
# create a vector and multiply
v = np.matrix([[2],[3],[4]])
v
Out[77]:
In [78]:
m*v
Out[78]:
更多的操作可以在numpy.linalg子模块中找到，例如：
In [79]:
import numpy as np
import numpy.linalg as nlg
#Determinant
nlg.det(m)
Out[79]:
In [80]:
#Eigenvalues
nlg.eigvals(m)
Out[80]:
In [81]:
#Solve for x in mx = v
x = nlg.solve(m,v)
x
Out[81]:
In [82]:
m*x
Out[82]:
In [83]:
v
Out[83]:
In [ ]: