NumPy基础：数组和矢量计算

0 为什么NumPy？——快

import numpy as np

my_arr = np.arange(1000000)
my_list = list(range(1000000))

# 序列分别乘以2,计算10次的时间，numpy表现明显更优
%time for i in range(10): my_arr2 = my_arr * 2
%time for i in range(10): my_list2 = [x * 2 for x in my_list]

1 Numpy的ndarray：一种多维数组对象

注： “数组”、“NumPy数组”、“ndarray” 指的都是 ndarray对象

import numpy as np

1.1 创建ndarray

data1 = [6, 7.5, 8, 0, 1]
arr1 = np.array(data1)
arr1, arr1.dtype

data2 = [[1, 2, 3, 4], [5, 6, 7, 8]]
arr2 = np.array(data2)
arr2, arr2.ndim, arr2.shape, arr2.dtype

np.random.randn(2, 3)

np.zeros(10)

np.ones((3,6))

np.empty((2,3,2))  #empty返回的是一些未经初始化的垃圾值

np.arange(15)  #arange是Python内置函数range的数组版

np.eye(5)  #生成N×N的单位阵，等同于np.identity()

np.full((2,3), 5)  #用5填满数组

arr_zeros = np.zeros((3,4))
np.ones_like(arr_zeros)  #类似还有zeros_like、empty_like、full_like

1.2 ndarray的数据类型

# 指定ndarray的数据类型
arr1 = np.array([1, 2, 3], dtype=np.float64)
arr2 = np.array([1, 2, 3], dtype=np.int32)

arr = np.array([3.7, -1.2, -2.6, 0.5, 12.9, 10.1]) #自动生成最匹配的数据类型
arr.dtype

arr.astype(np.int32)  # astype转换数据类型。
# 调用astype会创建一个新的数据（数据的备份），所以arr本身的数据类型没改变

numerical_strings = np.array(['1.25', '-9.6', '42'], dtype=np.string_)
numerical_strings.astype(arr1.dtype)

np.empty(4, dtype='u4')  # dtype类型代码

1.3 NumPy数组的运算

注：大小相等的数组之间的任何运算都会应用到元素级

arr1 = np.array([[1., 2., 3.], [4., 5., 6.]])
arr2 = np.array([[0., 4., 1.],[7., 2., 12.]])
# 元素级
# arr1 * arr2
# arr1 - arr2
# 1/arr1
# arr1 ** 0.5
# arr1 > arr2

1.4 基本的索引和切片

和列表不同，数组切片是原始数组的视图，这意味着数据不会被复制，视图上的任何修改都会直接反映到源数组上。如果想到得到ndarray切片的副本而非视图，就需要明确进行复制操作，如 arr[:].copy()

arr = np.arange(10)
arr[:] = 12
arr   # >> array([12, 12, 12, 12, 12])

# 没copy是视图
arr = np.arange(5)
arr_slice = arr[2:4]
arr_slice[1] = 9
arr   # >> array([0, 1, 2, 9, 4])

# copy是备份
arr = np.arange(5)
arr_slice = arr[2:4].copy()
arr_slice[1] = 9
arr  # >> array([0, 1, 2, 3, 4])

arr2d = np.array([[1,2, 3],[4, 5, 6],[7, 8, 9]])
# arr2d
# arr2d[2]
# arr2d[0][2]
# arr2d[0, 2]

1.5 切片索引

1.6 布尔型索引

names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
names == 'Bob'  # 布尔型数组

# 布尔型数组可用于数组索引
data = np.random.randn(7,4)
data[names == 'Bob']
data[names == 'Bob', 2:]
data[names != 'Bob']

# python关键字and和or在布尔型数组中无效，要使用&与|
cond = (names == 'Bob') | (names == 'Will')
data[cond]

data[data < 0] = 0
data

data[names != 'Joe'] = 7
data

1.7 花式索引

注：花式索引和切片不一样，它总是将数据复制到新数组中

arr = np.arange(12).reshape((3,4))
arr 
arr[[2,1]]
arr[[2,1], [1,3]] 
arr[[2,1]][:, [1,3]]

1.8 数组转置和轴对换

arr = np.arange(15).reshape((3,5))
arr.T

arr = np.arange(16).reshape((2, 2, 4))
arr.T   # 轴顺序 (0,1,2) -> (2,1,0)
arr.transpose((1, 0, 2)) # 轴顺序 (0,1,2) -> (1，0，2)
arr.swapaxes(0, 1)  # 同上，swapaxes也是返回源数据的视图，不进行复制

2 通用函数：快速的元素级数组函数

注：通用函数（ufunc）是一种对ndarray中的数据执行元素级运算的函数

3 利用数组进行数据处理

import matplotlib.pyplot as plt

points = np.arange(-5, 5, 0.01)  #在[-5, 5)之间以0.01为步长的1000个点
x, y = np.meshgrid(points, points)
z = np.sqrt(xs ** 2 + ys ** 2)  #计算横纵坐标的 sqrt(平方和)

plt.imshow(z, cmap=plt.cm.RdBu);
plt.colorbar
plt.title("Image plot of $\sqrt{x^2 + y^2}$ for a grid of values")

3.1 将条件逻辑表述为数组运算

xarr = np.array([1, 2, 3, 4, 5])
yarr = np.array([11, 12, 13, 14, 15])
cond = np.array([True, False, True, True, False])

# 纯python方法
result = [(x if c else y) for x, y, c in zip(xarr, yarr, cond)]
result

# np.where方法
result = np.where(cond, xarr, yarr)
result

3.2 数学和统计方法

sum、mean、std、var、max、min、argmax、argmin、cumsum、cumprod

arr = np.random.randn(2, 4)
arr.mean()
np.mean(arr)
arr.mean(axis=1)  # 指定的轴就是运算的方向
arr.sum()
arr.sum(axis=0)

arr = np.array([0, 1, 2, 3, 4, 5, 6, 7])
arr.cumsum()  # 累加函数

arr = np.arange(9).reshape(3, 3)
arr.cumsum(axis=0)
arr.cumprod(axis=1)  #累乘函数

3.3 用于布尔型数组的方法

arr = np.random.randn(100)
(arr>0).sum()  #计算正项的个数

bools = np.array([False, False, True, False])
bools.any()  # 检查数组中是否有True，所有非零元素都会被当作True
bools.all()  # 检查数组中是否都是True

3.4 排序

arr = np.random.randn(6)
arr.sort()

arr = np.random.randn(5, 3)
arr.sort(axis=1)

large_arr = np.random.randn(1000)
large_arr.sort()
large_arr[int(0.05 * len(large_arr))]  # 百分位数

3.5 唯一化以及其他的几何逻辑

NumPy中的集合函数：

• unique(x)：计算x中的唯一元素，并返回有序结果
• intersect1d(x,y)：计算x和y的交集，并返回有序结果
• union1d：计算x和y的并集，并返回有序结果
• in1d(x,y)：得到一个表示“x的元素是否包含于y”的布尔型数组
• setdiff1d(x,y)：计算x和y的差集，即在x中不在y中
• setxor1d(x,y)：集合的对称差

names = np.array(['Bob', 'Joe', 'Will', 'Bob', 'Will', 'Joe', 'Joe'])
%time np.unique(names)  # np.unique唯一化并排序
# 纯python方法
%time sorted(set(names))  # set唯一化，sorted排序

values = np.array([6, 0, 0, 3, 2, 5, 6])
np.in1d(values, [2, 3, 6])  #x的元素是否包含于y

4 用于数组的文件输入输出

arr = np.arange(10)
np.save('some_array', arr)  #如果不加扩展名.npy，则会被自动加上
np.load('some_array.npy')

#savez可以将多个数组保存到一个未压缩文件中
np.savez('array_archive.npz', a=arr, b=arr)
arch = np.load('array_archive.npz')
arch['b']

#要将数据压缩，可以使用savez_compressed
np.savez_compressed('array_compressed.npz', a=arr, b=arr)

5 线性代数

常用的numpy.linalg函数：

• diag：对角线，或将一维数组转换成方阵
• dot：矩阵乘法
• trace：迹，对角线元素的和
• det：行列式
• eig：特征值和特征向量
• inv：方阵的逆
• pinv：计算矩阵的Moore-Penrose违逆
• qr：计算QR分解
• svd：计算奇异值分解（SVD）
• solve：解线性方程组Ax+b=b
• lstsq：计算Ax=b的最小二乘解

from numpy.linalg import inv, qr, svd

X = np.random.randn(5, 5)
inv(X)  #逆
q, r = qr(X)  # QR分解
s, v, d = svd(X)  # svd 分解

# 对于np.ndarray类型的A和B
A * B   # 元素级乘法
A @ B   # 矩阵乘法
np.multiply(A, B) # 元素级乘法
np.dot(A, B)  # 矩阵乘法

# 对于np.matrix类型的A和B (矩阵类型)
A * B   # 矩阵乘法
A @ B   # 矩阵乘法
np.multiply(A, B) # 元素级乘法
np.dot(A, B)  # 矩阵乘法

6 伪随机数生成

部分numpy.random函数：

• seed：确定随机数生成器的种子
• permutation：返回一个序列的随机排列或返回一个随机排列的范围
• shuffle：对一个序列就地随机排列
• rand：产生均匀分布的样本值
• randint：从给定的上下限范围内随机选取整数
• randn：产生标准正态分布的样本值
• binomial：产生二项分布的样本值
• normal：产生正态分布的样本值
• bata：产生Beta分布的样本值
• chisquare：产生卡方分布的样本值
• gamma：产生Gamma分布的样本值
• uniform：产生在[0,1)中均匀分布的样本值

# np.random
np.random.rand()     # [0,1)之间的随机数
np.random.rand(3,2)  # 形状为3*2的[0,1)之间的随机数
np.random.randn()    # 服从标准正态分布的随机数
np.random.randn(3,2) # 形状为3*2的服从标准正态分布的随机数
np.random.randint(2,5)  # [2,5)之间的随机整数
np.random.randint(2,size=5)  # [0,2)之间的5个随机整数
np.random.choice(np.arange(4),
                 size=(3,3),
                 p=[0.1,0.2,0.3,0.4], 
                 replace=True) # p指定了每个元素被选择的概率, repalce=True表示可以重复
np.random.seed(1)
np.random.rand()  # >> 0.417022004702574
np.random.seed(1)
np.random.rand()  # >> 0.417022004702574, 同一个seed生成的随机数相同

7 示例：随机漫步

随机漫步

# 纯python方法
import random
import matplotlib.pyplot as plt
position = 0
walk = [position]
steps = 1000
for i in range(steps):
    step = 1 if random.randint(0, 1) else -1
    position += step
    walk.append(position)
%time plt.plot(walk[:100])

# 用numpy库
import numpy as np
import matplotlib.pyplot as plt

position = 0
nsteps = 1000
step = np.random.choice([1,-1], size=nsteps)
walk = step.cumsum()

%time plt.plot(walk[:100])

(np.abs(walk)>=10).argmax()  #首达时

一次模拟多个随机漫步

#一次性模拟5000个随机漫步
import numpy as np
import matplotlib.pyplot as plt

nwalks = 5000
nsteps = 1000
steps = np.random.choice([1,-1], size=(nwalks, nsteps))
walks = steps.cumsum(1)

walks
walks.max()
walks.min()
hits30 = (np.abs(walks) >= 30).any(1)
hits30.sum()
crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)  #所有绝对值达到30的随机漫步的首达时