4.2 通用函数:快速的逐元素数组函数
通用函数,也可以称为ufunc,是一种在ndarray数据中进行逐元素操作的函数。某些简单函数接受一个或多个标量数值,并产生一个或多个标量结果,而通用函数就是对这些简单函数的向量化封装。
有很多ufunc是简单的主元素转换,比如sqrt或exp函数:
In [100]: arr = np.arange(10)
In [101]: arr
Out[101]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [102]: np.sqrt(arr)
Out[102]:
array([0. , 1. , 1.41421356, 1.73205081, 2. ,
2.23606798, 2.44948974, 2.64575131, 2.82842712, 3. ])
In [106]: np.exp(arr)
Out[106]:
array([1.00000000e+00, 2.71828183e+00, 7.38905610e+00, 2.00855369e+01,
5.45981500e+01, 1.48413159e+02, 4.03428793e+02, 1.09663316e+03,
2.98095799e+03, 8.10308393e+03])
那么sqrt和exp到底是啥意思呢……
sqrt就是平方根,就是√ ̄。
exp就是e的X次方,0次方是1,1次方是2.7183,2次方是7.389。
这个e就是数学中的一个常数,是一个无限不循环小数。具体的建议去看百度百科中对于e的解释:
地址如下:自然常数
这些是所谓的一元通用函数。还有一些通用函数,比如add或maximum则会接收两个数组并返回一个数组作为结果,因此成为二元通用函数:
In [119]: y = np.random.random(10)
In [120]: x = np.random.random(10)
In [121]: x
Out[121]:
array([0.84011849, 0.87387061, 0.76205089, 0.69165302, 0.34540561,
0.30167627, 0.46833912, 0.84250118, 0.59886878, 0.64750627])
In [122]: y
Out[122]:
array([0.72860809, 0.35024632, 0.62468672, 0.94319557, 0.32032339,
0.56525057, 0.12750191, 0.46299197, 0.18831764, 0.20538675])
In [123]: np.maximum(x,y)
Out[123]:
array([0.84011849, 0.87387061, 0.76205089, 0.94319557, 0.34540561,
0.56525057, 0.46833912, 0.84250118, 0.59886878, 0.64750627])
这里numpy.maximum逐个元素地将X和y中元素的最大值计算出来。
maximum:对比两组数值(X,Y),输出数值较大的数值(0.84 > 0.72,所以生成0.84)。
也有一些通用函数返回多个数组。比如modf,是Python内建函数divmod的向量化版本。它返回了一个浮点值数组的小数部分和整数部分。
In [126]: arr = np.random.random(10) * 5
In [127]: arr
Out[127]:
array([1.514357 , 2.67348413, 2.38685845, 0.3327048 , 1.249205 ,
2.68107619, 3.7438636 , 2.57037747, 2.54322009, 1.21793862])
In [128]: a,b = np.modf(arr)
In [129]: a
Out[129]:
array([0.514357 , 0.67348413, 0.38685845, 0.3327048 , 0.249205 ,
0.68107619, 0.7438636 , 0.57037747, 0.54322009, 0.21793862])
In [130]: b
Out[130]: array([1., 2., 2., 0., 1., 2., 3., 2., 2., 1.])
modf比较好理解,就是一个数字,比如3.1415926,它由整数和小数点后面的数字组成。modf就是把整数和小数拆开,变成3和0.1415926。
通用函数接收一个可选参数out,允许对数组按位置操作:
In [132]: arr
Out[132]:
array([1.514357 , 2.67348413, 2.38685845, 0.3327048 , 1.249205 ,
2.68107619, 3.7438636 , 2.57037747, 2.54322009, 1.21793862])
In [133]: np.sqrt(arr)
Out[133]:
array([1.23059213, 1.63507924, 1.5449461 , 0.57680568, 1.1176784 ,
1.63739921, 1.93490661, 1.60323968, 1.59474766, 1.10360256])
In [134]: np.sqrt(arr,arr)
Out[134]:
array([1.23059213, 1.63507924, 1.5449461 , 0.57680568, 1.1176784 ,
1.63739921, 1.93490661, 1.60323968, 1.59474766, 1.10360256])
In [135]: arr
Out[135]:
array([1.23059213, 1.63507924, 1.5449461 , 0.57680568, 1.1176784 ,
1.63739921, 1.93490661, 1.60323968, 1.59474766, 1.10360256])
这个看着有点云里雾里的,书上的例子比较难懂,我换一个方式吧:
In [135]: arr
Out[135]:
array([1.23059213, 1.63507924, 1.5449461 , 0.57680568, 1.1176784 ,
1.63739921, 1.93490661, 1.60323968, 1.59474766, 1.10360256])
In [138]: np.sqrt(0,arr)
Out[138]: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
In [139]: arr
Out[139]: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
意思就是通用函数增加一个可选参数,然后直接赋值给原来的变量。np.sqrt(0,arr)中,我们把0赋给arr,然后再开平方根,0的平方根是0。就是这个意思。
4.3 使用数组进行面向数组编程
使用Numpy数组可以使你利用简单的数组表达式完成多种数据操作任务,而无需写些大量循环。这种利用数组表达式来替代显式循环的方法,称之为向量化。通常,向量化的数组操作会比纯Python的等价实现在速度上快一到两个数量级(甚至更多),这对所有种类的数值计算产生了最大的影响。
作为一个简单的示例,假设我们想要对一些网格数据计算函数sqrt(x^2 +y^2)的值。np.meshgrid函数接收两个一维数组,并根据两个数组的所有(x,y)对生成一个二维矩阵:
In [149]: po = np.arange(-5,5,0.01)
In [150]: xs,ys = np.meshgrid(po,po)
In [152]: ys
Out[152]:
array([[-5. , -5. , -5. , ..., -5. , -5. , -5. ],
[-4.99, -4.99, -4.99, ..., -4.99, -4.99, -4.99],
[-4.98, -4.98, -4.98, ..., -4.98, -4.98, -4.98],
...,
[ 4.97, 4.97, 4.97, ..., 4.97, 4.97, 4.97],
[ 4.98, 4.98, 4.98, ..., 4.98, 4.98, 4.98],
[ 4.99, 4.99, 4.99, ..., 4.99, 4.99, 4.99]])
现在,你可以用和两个坐标值相同的表达式来使用函数:
In [157]: z
Out[157]:
array([[7.07106781, 7.06400028, 7.05693985, ..., 7.04988652, 7.05693985,
7.06400028],
[7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.04985815,
7.05692568],
[7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.04278354,
7.04985815],
...,
[7.04988652, 7.04279774, 7.03571603, ..., 7.0286414 , 7.03571603,
7.04279774],
[7.05693985, 7.04985815, 7.04278354, ..., 7.03571603, 7.04278354,
7.04985815],
[7.06400028, 7.05692568, 7.04985815, ..., 7.04279774, 7.04985815,
7.05692568]])
作为第9章的前瞻,我使用matplotlib来生成这个二维数组的可视化:
import numpy as np
import matplotlib.pyplot as plt
po = np.arange(5,-5,0.01)
ys,xs = np.meshgrid(po,po)
z = np.sqrt(xs**2 + ys**2)
plt.imshow(z,cmap=plt.cm.gray);plt.colorbar()
plt.title("Image plot of $\sqrt{x^2 + y^2}$ for a grid of values")
这块按照书上写的代码会报错……并生成不了图,也不知道为什么。百度上也没有相关的资料,B站up主致敬大神直接就略过这块了。
难受……那就略过吧。
4.3.1 将条件逻辑作为数组操作
numpy.where函数是三元表达式 x if condition else y的向量化版本。假设我们有一个布尔值数组和两个数值数组:
In [165]: xarr = np.array([1.1,1.2,1.3,1.4,1.5])
In [166]: yarr = np.array([2.1,2.2,2.3,2.4,2.5])
In [167]: cond = np.array([True,False,True,True,False])
假设cond中的元素为True时,我们取xarr中的对应元素值,否则yarr中的元素。我们可以通过列表推导式来完成,像下面这样:
In [171]: re = [(x if c else y) # 选择x的值,如果是c(=True)则选择y。
...: for x,y,c in zip(xarr,yarr,cond)]
In [172]: re
Out[172]: [1.1, 2.2, 1.3, 1.4, 2.5]
这会产生多个问题。首先,如果数组很大的话,速度会变很慢(因为所有的工作都是通过解释器解释Python代码完成)。其次,当数组是多维时,就无法奏效了。而使用np.where时,就可以非常简单地完成:
In [173]: re1 = np.where(cond,xarr,yarr)
In [174]: re1
Out[174]: array([1.1, 2.2, 1.3, 1.4, 2.5])
np.where的第二个和第三个参数并不需要是数组,它们可以是标量。where在数据分析中的一个典型用法就是根据一个数组来生成一个新的数组。假设你有一个随机生成的矩阵数据,并且你想将其中的正值都替换为2,将所有的负值替换为-2,使用np.where会很容易实现。
In [177]: arr
Out[177]:
array([[-1.37021712, -0.87138269, -0.15010701, 0.22736407],
[-2.38806998, 0.11452806, -0.01859951, -0.3149819 ],
[ 0.81163732, 0.84179709, -0.28342479, 0.23882525],
[-0.03122251, 0.50378108, 0.83066658, -0.4822985 ]])
In [178]: arr > 0
Out[178]:
array([[False, False, False, True],
[False, True, False, False],
[ True, True, False, True],
[False, True, True, False]])
In [179]: np.where(arr>0,2,-2)
Out[179]:
array([[-2, -2, -2, 2],
[-2, 2, -2, -2],
[ 2, 2, -2, 2],
[-2, 2, 2, -2]])
你可以使用np.where将标量和数组联合,例如,可以像下面的代码那样arr中的所有正值替换为常数2:
In [187]: np.where(arr>0,2,arr)
Out[187]:
array([[-1.37021712, -0.87138269, -0.15010701, 2. ],
[-2.38806998, 2. , -0.01859951, -0.3149819 ],
[ 2. , 2. , -0.28342479, 2. ],
[-0.03122251, 2. , 2. , -0.4822985 ]])
4.3.2 数学和统计方法
许多关于计算整个数组统计值或关于轴向数据的科学函数,可以作为数组类型的方法被调用。你可以使用聚合函数(通常也叫缩减函数),比如sum、mean和std(标准差),既可以直接调用数组实例的方法,也可以使用顶层的numpy函数。
此处我生成了一个正态分布的随机数,并计算了部分聚合统计数据:
In [188]: arr = np.random.randn(5,4)
In [189]: arr
Out[189]:
array([[-1.29373731, 0.03123039, -0.96616799, -0.8496834 ],
[ 0.68013113, 0.63368681, 1.10114311, -1.02762889],
[ 0.65165446, 0.46103685, 1.55264834, 0.87102779],
[ 1.91937296, -1.42439436, 0.30368359, 1.00357199],
[ 0.33076767, -0.9469198 , 0.93020695, 1.07798312]])
In [190]: arr.mean() # 平均值
Out[190]: 0.2519806710822715
In [191]: np.mean(arr)
Out[191]: 0.2519806710822715
In [192]: arr.sum() # 元素的累和
Out[192]: 5.03961342164543
像mean、sum等函数可以接收一个可选参数axis,这个参数可以用于计算给定轴向上的统计值,形成一个下降一维度的数组:
In [193]: arr.mean(axis=1)
Out[193]: array([-0.76958958, 0.34683304, 0.88409186, 0.45055855, 0.34800948])
In [194]: arr.sum(axis=0)
Out[194]: array([ 2.28818891, -1.2453601 , 2.921514 , 1.07527061])
arr.mean(1)表示“计算每一列的平均值”,而arr.sum(0)表示“行轴向的和”(跟上一篇文章讲的一样,1就是列,0就是行)。
其他的方法,例如cumsum和cumprod并不会聚合,它们会产生一个中间结果:
In [195]: arr = np.array([[0,1,2],[3,4,5],[6,7,8]])
In [196]: arr
Out[196]:
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
In [197]: arr.cumsum(axis=0)
Out[197]:
array([[ 0, 1, 2],
[ 3, 5, 7],
[ 9, 12, 15]])
# 0+3 = 3, 3+3 = 6 1+4 = 5,5+7=12 2+5=7,7+8=15
4.3.3 布尔值数组的方法
在前面介绍的方法,布尔值会被强制为1(True)和0(False)。因此,sum通常可以用于计算布尔值数组中的True的个数:
In [202]: arr = np.random.randn(100)
In [203]: (arr>0).sum() # 正值的个数
Out[203]: 49
对于布尔值数组,有两个非常有用的方法any和all。any检查数组中是否至少有一个true,而all检查是否每个值都是true:
In [205]: bools = np.array([False,False,True,False])
In [206]: bools.any()
Out[206]: True
In [207]: bools.all()
Out[207]: False
这些方法也可适用于非布尔值数组,所有的非0元素都会按True处理。
4.3.4 排序
和Python的内建列表类型相似,numpy数组可以使用sort方法按位置排序:
In [208]: arr = np.random.randn(6)
In [209]: arr
Out[209]:
array([ 0.71143173, 0.31392501, -1.50800015, 0.47079999, -1.17964639,
-0.51921076])
In [210]: arr.sort()
In [211]: arr
Out[211]:
array([-1.50800015, -1.17964639, -0.51921076, 0.31392501, 0.47079999,
0.71143173])
你可以在多维数组中根据传递的axis值,沿着轴向对每个一维数据段进行排序:
In [214]: arr
Out[214]:
array([[ 0.53034875, -0.14207949, -1.0790517 ],
[-1.70692011, 0.72902593, -1.43168883],
[ 0.43739608, 1.66885615, 1.68066729],
[-1.73041806, -0.70178695, -0.08415239],
[-0.68712738, -0.01258072, 0.12001392]])
In [215]: arr.sort(1)
In [216]: arr
Out[216]:
array([[-1.0790517 , -0.14207949, 0.53034875],
[-1.70692011, -1.43168883, 0.72902593],
[ 0.43739608, 1.66885615, 1.68066729],
[-1.73041806, -0.70178695, -0.08415239],
[-0.68712738, -0.01258072, 0.12001392]])
# 每一行的排序
顶层的np.sort方法返回的是已经排序好的数组拷贝,而不是对原数组按照位置排序。下面的例子计算的是一个数组的分位数,并选出分位数所对应的值,这是一种应急的方式:
In [217]: large = np.random.randn(1000)
In [218]: large.sort()
In [219]: large[int(0.05*len(large))]
Out[219]: -1.6398203401470943
这里需要注意,randn(1000)的意思是创建1000个符合正态分布的数值,然后进行排序,在第50个数字的数是-1.63982034……
4.3.5 唯一值与其他集合逻辑
Numpy包含一些针对一维ndarray的基础集合操作。常用的一个方法是np.unique,返回的是数组中唯一值排序后形成的数组:
In [224]: name = np.array(['quanxiaosheng','piaozhiyan','piaoxiaomin','piaozhiyan','jinzhenxi','jinxuanya','jinzhenxi'])
In [225]: np.unique(name)
Out[225]:
array(['jinxuanya', 'jinzhenxi', 'piaoxiaomin', 'piaozhiyan',
'quanxiaosheng'], dtype='<U13')
In [226]: int = np.array([1,3,2,5,3,2,1,5,6,7,8,5,6])
In [227]: np.unique(int)
Out[227]: array([1, 2, 3, 5, 6, 7, 8])
将np.unique和纯Python实现相比较:
In [228]: sorted(set(name))
Out[228]: ['jinxuanya', 'jinzhenxi', 'piaoxiaomin', 'piaozhiyan', 'quanxiaosheng']
另一个函数,np.in1d,可以检查一个数组中的值是否在另外一个数组中,并返回一个布尔值数组:
In [229]: var = np.array([6,0,0,3,2,5,6])
In [230]: np.in1d(var,[2,3,6])
Out[230]: array([ True, False, False, True, True, False, True])
这里就是拿var里面的每个数字去后面的数组中轮询,得出布尔值。6在236里有,就是true,0在236里没有,就是false。
4.4 使用数组进行文件输入和输出
Numpy可以在硬盘中将数据以文本或二进制文件的形式进行存入硬盘或由硬盘载入。在本节,我将只讨论Numpy的内建二进制格式,因为大部分用户更倾向于使用pandas或其他工具来载入文本或表格型数据。
np.save和np.load是高效存取硬盘数据的两大工具函数。数组在默认情况下是以未压缩的格式进行存储的,后缀名是.npy。
In [231]: arr = np.array(range(10))
In [232]: np.save('some_array',arr)
如果文件存放路径中没写.npy时,后缀名会被自动加上。硬盘上的数组可以使用np.load进行载入:
In [236]: np.load('some_array.npy')
Out[236]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
你可以使用np.savez并将数组作为参数传递给该函数,用于在未压缩文件中保存多个数组:
np.savez('array.npz',a=arr,b=arr)
当载入一个.npy文件的时候,你会获得一个字典型的对象,并通过该对象很方便地载入单个数组:
In [240]: arch = np.load('array.npz')
In [241]: arch['b']
Out[241]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
如果你的数据已经压缩好了,你可能会想要使用numpy.savez_compressed将数据存入已经压缩的文件:
np.savez_compressed('array_1.npz',a=arr,b=arr)
这里额外补充一点内容:
当我们将内容打包进npz之后,我们是没有办法直接打开npz的,只能通过’key’也就是上述代码中的[‘b’]去查看。但是如果我们不知道这个npz到底有哪些文件,也不知道这个npz里面的key是什么该怎么办呢?
这里补充一个方法:
.files
In [257]: np.load('array_1.npz').files
Out[257]: ['a', 'b']
4.5 线性代数
线性代数,比如矩阵乘法、分解、行列式等方阵数学,是所有数组类库的重要组成部分(胭惜雨:幸好不是全部……)。和matlab等其他语言相比,numpy的线性代数中所不同的是*是矩阵的逐元素乘积,而不是矩阵的点乘积。因此numpy的数组方法和numpy命名空间中都有一个函数dot,用于矩阵的操作:
关于矩阵内积,今天我苦思冥想,终于有点思路了,还拿昨天的数据来说吧。
In [84]: arr
Out[84]:
array([[ 0, 1, 2, 3, 4],
[ 5, 6, 7, 8, 9],
[10, 11, 12, 13, 14]])
In [85]: arr.T
Out[85]:
array([[ 0, 5, 10],
[ 1, 6, 11],
[ 2, 7, 12],
[ 3, 8, 13],
[ 4, 9, 14]])
In [86]: old = arr.copy()
In [87]: np.dot(arr.T,arr)
Out[87]:
array([[125, 140, 155, 170, 185],
[140, 158, 176, 194, 212],
[155, 176, 197, 218, 239],
[170, 194, 218, 242, 266],
[185, 212, 239, 266, 293]])
125 = 0*0 + 5*5 + 10*10
140 = 1*0 + 6*5 + 11*10
第一行内容以此类推。那么到第四行呢?
170 = 0*3 + 5*8 + 10*13
如果我们要找212这个值,应该从哪找呢?
ok,差不多明白了吧,接下来我们进入今天的数据:
In [260]: x = np.array([[1.,2.,3.],[4.,5.,6.]])
In [261]: y = np.array([[6.,23.],[-1,7],[8,9]])
In [262]: x
Out[262]:
array([[1., 2., 3.],
[4., 5., 6.]])
In [263]: y
Out[263]:
array([[ 6., 23.],
[-1., 7.],
[ 8., 9.]])
In [264]: x.dot(y)
Out[264]:
array([[ 28., 64.],
[ 67., 181.]])
x.dot(y)等价于np.dot(x,y):
In [264]: x.dot(y)
Out[264]:
array([[ 28., 64.],
[ 67., 181.]])
In [265]: np.dot(x,y)
Out[265]:
array([[ 28., 64.],
[ 67., 181.]])
一个二维数组和一个长度合适的一位数组之间的矩阵乘积,其结果是一个一维数组:
In [266]: np.dot(x,np.ones(3)) #ones(3) = array([1., 1., 1.])
Out[266]: array([ 6., 15.])
特殊符号@也作为中缀操作符,用于点成矩阵操作:
In [269]: x @ np.ones(3)
Out[269]: array([ 6., 15.])
numpy.linalg 拥有一个矩阵分解的标注函数集,以及其他常用函数,例如求逆和行列式求解。这些函数都是通过在matlab和r等其他语言使用的相同的行业标准线性代数库实现的,例如BLAS、LAPACK或英特尔专有的MKL(数学核心库)(是否使用MKL取决于使用numpy的版本):
In [271]: from numpy.linalg import inv,qr
In [272]: x =np.random.randn(5,5)
In [273]: mat = x.T.dot(x)
In [274]: inv(mat)
Out[274]:
array([[ 0.41146636, 0.14102918, -0.38213338, 0.09299996, -0.68646348],
[ 0.14102918, 0.23375719, -0.18007064, 0.16723144, -0.6318615 ],
[-0.38213338, -0.18007064, 1.19473588, -0.31378151, 1.96712614],
[ 0.09299996, 0.16723144, -0.31378151, 0.55961514, -1.10198253],
[-0.68646348, -0.6318615 , 1.96712614, -1.10198253, 5.20331582]])
In [275]: mat.dot(inv(mat))
Out[275]:
array([[ 1.00000000e+00, 3.93334552e-17, -9.03268796e-17,
7.35113472e-18, 5.43524235e-17],
[-6.80240673e-17, 1.00000000e+00, -3.19504841e-16,
1.10155711e-16, -2.34266807e-16],
[-5.28357836e-17, -3.39184656e-17, 1.00000000e+00,
5.16667005e-17, 7.51371078e-17],
[ 1.36509468e-17, -4.02066164e-17, 3.78861318e-16,
1.00000000e+00, 8.14882255e-16],
[-2.36362053e-17, -3.41065208e-17, 3.02067825e-16,
-8.25166167e-17, 1.00000000e+00]])
In [276]: q,r = qr(mat)
In [277]: r
Out[277]:
array([[-4.73157356, 5.61240256, -2.43122059, -0.85451614, 0.7993408 ],
[ 0. , -6.09867896, 0.82343319, 0.33031748, -1.01204728],
[ 0. , 0. , -2.65419264, 2.02744257, 1.50751986],
[ 0. , 0. , 0. , -2.88378432, -0.62317305],
[ 0. , 0. , 0. , 0. , 0.17400141]])
表达式x.T.dot(X)计算的是X和它的转置矩阵X.T的点乘积。
4.6 伪随机数生成
numpy.random模块填补了Python内建的random模块的不足,可以高效地生成多种概率分布下的完整样本值数组。例如,你可以使用normal来获得一个4X4的正态分布样本数组:
In [278]: sam = np.random.normal(size=(4,4))
In [279]: sam
Out[279]:
array([[-0.03016345, 1.15460783, -0.21663213, -0.25145808],
[-0.38961357, 0.075285 , 0.93884707, 1.00314845],
[ 1.13388719, 1.26916511, -0.9402673 , 0.66699758],
[ 0.27428168, 1.78101786, -0.02000926, 0.06556419]])
然而Python内建的random模块一次只能生成一个值。你可以从下面的示例中看到,numpy.random在生成大型样本时比纯Python的方法快了一个数量级:
from random import normalvariate
In [292]: N = 100000
In [293]: %timeit sam = [normalvariate(0,1) for _ in range(N)]
458 ms ± 45.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
In [294]: %timeit np.random.normal(size=N)
16.7 ms ± 638 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
我们可以称这些为伪随机数,因为它们是由具有确定性行为的算法根据随机数生成器中的随机数种子生成的。你可以通过np.random.seed 更改numpy的随机数种子:
np.random.seed(1234)
numpy.random中的数据生成函数使用一个全局随机数种子。为了避免全局状态,你可以使用numpy.random.RandomState创建一个随机数生成器,使数据独立于其他的随机数状态:
In [298]: rng = np.random.RandomState(1234)
In [299]: rng.randn(10)
Out[299]:
array([ 0.47143516, -1.19097569, 1.43270697, -0.3126519 , -0.72058873,
0.88716294, 0.85958841, -0.6365235 , 0.01569637, -2.24268495])
补充:
gamma :从伽马分布中抽取样本
uniform:从均匀[0,1)分布中抽取样本
4.7 示例:随机漫步
随机漫步的模拟提供了一个使用数组操作的说明性应用。首先,让我们考虑一个简单的随机漫步,从0开始,步进为1和-1,且两种不仅发生的概率相等。
import matplotlib.pyplot as plt
import random
pos = 0
walk = [pos]
steps = 1000
for i in range(steps):
step = 1 if random.randint(0, 1) else -1
pos += step
walk.append(pos)
plt.plot(walk[:100])
plt.show()
你可能会观察到walk只是对随机步进的累积,并且可以通过一个数组表达式实现。因此,我使用np.random模块一次性抽取1000次投掷硬币的结果,每次投掷的结果为1或-1,然后计算累积值:
In [331]: ns = 1000
In [332]: draws = np.random.randint(0,2,size = ns)
In [333]: set = np.where(draws>0,1,-1)
In [334]: walk = set.cumsum()
由此我们开始从漫步轨道上提取一些统计数据,比如最大值、最小值等:
In [335]: walk.min()
Out[335]: -9
In [336]: walk.max()
Out[336]: 60
更复杂的统计是第一次穿越时间,即随机漫步的某一步达到了某个特定值。这里假设我们想要知道漫步中是何时连续朝某个方向连续走了十步。np.abs(walk)>=10给我们一个布尔值数组,用于表明漫步是否连续在同一个方向走了十步,但是我们想要的是第一次走了十步或者-10步的位置。我们可以使用argmax来计算,该函数可以返回布尔值数组中最大值的第一个位置:
In [337]: (np.abs(walk)>= 10).argmax()
Out[337]: 297
请注意,这里使用argmax效率并不高,因为它总是完整地扫描整个数组。在这个特殊的示例中,一旦True被发现,我们就知道最大值了。
4.7.1 一次性模拟多次随机漫步
如果你的目标是模拟多次随机漫步,比如说5000步。你可以稍微地修改下之前的代码来生成所有的随机步。如果传入一个2个元素的元组,numpy.random中的函数可以生成一个二维的抽取数组,并且我们可以一次性地跨越行算出全部5000个随机步的累积和:
import numpy as np
nwalks = 5000
nsteps = 1000
draws = np.random.randint(0,2,size=(nwalks,nsteps))
steps = np.where(draws>0,1,-1)
walks = steps.cumsum(1)
print(walks)
[[ 1 0 -1 ... -58 -57 -58]
[ 1 2 1 ... -86 -85 -84]
[ 1 0 1 ... 38 39 38]
...
[ -1 -2 -3 ... 42 41 40]
[ 1 0 -1 ... -10 -9 -8]
[ 1 2 1 ... 38 39 38]]
现在我们可以计算出这些随机步的最大值和最小值了:
walks.max()
walks.min()
让我们在这些随机步中计算出30 或-30的最小穿越时间。这有点棘手,因为不是所有的5000个都达到了30.我们可以使用any方法来检查:
In [343]: hit30 = (np.abs(walks)>= 30).any(1) #abs()就是取绝对值
In [344]: hit30
Out[344]: array([ True, True, True, ..., True, False, True])
In [345]: hit30.sum()
Out[345]: 3368
我们可以使用布尔值数组来选出绝对步数超过30的步所在的行,并使用argmax从轴向1上获取穿越时间:
In [346]: corssing_times = (np.abs(walks[hit30])>30).argmax(1)
In [348]: corssing_times.mean()
Out[348]: 495.1543942992874
利用其它分布而不是等概率的掷硬币试验来随机漫步也是很容易的。你只需要使用一个不同的随机数生成函数,比如normal,再根据特定的均值和标准差即可生成正态分布下的随机步:
In [350]: steps = np.random.normal(loc=0,scale=0.25,size=(nwalks,nsteps))
本章后半段的一部分内容在之前《Python编程》中有描述,不过更进了一步,有了一个拔高的过程。
不过受限于某些原因,今天就到这里了,明天见~
胭惜雨
2021年03月24日