12-11第四章3 NumPy基础:数组和矢量计算

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

Numpy能从磁盘直接存储和加载数据,不论是文本格式还是二进制模式。这里我们只考虑Numpy的二进制模式,因为大多数用户更喜欢用pandas或其他工具来加载text或tabular数据。

np.save和np.load。数组会以未压缩的原始二进制模式被保存,后缀为.npy:

In [213]: arr = np.arange(10)

In [214]: np.save('some_array', arr)

如果文件路径末尾没有扩展名.npy,则该扩展名会被自动加上。然后就可以通过np.load读取磁盘上的数组:

In [215]: np.load('some_array.npy')
Out[215]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

用np.savez能保存多个数组,还可以指定数组对应的关键字,不过是未压缩的npz格式:

In [216]: np.savez('array_archive.npz', a=arr, b=arr)

加载.npz文件时,你会得到一个类似字典的对象,该对象会对各个数组进行延迟加载:

In [217]: arch = np.load('array_archive.npz')

In [218]: arch['b']
Out[218]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

如果要将数据压缩,可以使用numpy.savez_compressed:

In [219]: np.savez_compressed('arrays_compressed.npz', a=arr, b=arr)

4.5 线性代数

线性代数(如矩阵乘法、矩阵分解、行列式以及其他方阵数学等)是任何数组库的重要组成部分。不像某些语言(如MATLAB),通过*对两个二维数组相乘得到的是一个元素级的积,而不是一个矩阵点积。因此,NumPy提供了一个用于矩阵乘法的dot函数(既是一个数组方法也是numpy命名空间中的一个函数):

In [223]: x = np.array([[1., 2., 3.], [4., 5., 6.]])

In [224]: y = np.array([[6., 23.], [-1, 7], [8, 9]])

In [225]: x
Out[225]: 
array([[ 1.,  2.,  3.],
       [ 4.,  5.,  6.]])

In [226]: y
Out[226]: 
array([[  6.,  23.],
       [ -1.,   7.],
       [  8.,   9.]])

In [227]: x.dot(y)
Out[227]: 
array([[  28.,   64.],
       [  67.,  181.]])

x.dot(y)等价于np.dot(x, y):

In [228]: np.dot(x, y)
Out[228]: 
array([[  28.,   64.],
       [  67.,  181.]])

一个二维数组跟一个大小合适的一维数组的矩阵点积运算之后将会得到一个一维数组:

In [229]: np.dot(x, np.ones(3))# 3个是1的矩阵
#ones(shape[, dtype, order] 依据给定形状和类型(shape[, dtype, order])返回一个新的元素全部为1的数组。
Out[229]: array([  6.,  15.])

@符(类似Python 3.5)也可以用作中缀运算符,进行矩阵乘法:

In [230]: x @ np.ones(3)
Out[230]: array([  6.,  15.])

np.linalg能用来做矩阵分解,以及比如转置和求秩之类的事情:

In [231]: from numpy.linalg import inv, qr

In [232]: X = np.random.randn(5, 5)

In [233]: mat = X.T.dot(X) #X.T.dot(X)计算的是X和X的转置的矩阵乘法。

In [234]: inv(mat)  #inv是求逆矩阵
Out[234]: 
array([[  933.1189,   871.8258, -1417.6902, -1460.4005,  1782.1391],
       [  871.8258,   815.3929, -1325.9965, -1365.9242,  1666.9347],
       [-1417.6902, -1325.9965,  2158.4424,  2222.0191, -2711.6822],
       [-1460.4005, -1365.9242,  2222.0191,  2289.0575, -2793.422 ],
       [ 1782.1391,  1666.9347, -2711.6822, -2793.422 ,  3409.5128]])

In [235]: mat.dot(inv(mat))
Out[235]: 
array([[ 1.,  0., -0., -0., -0.],
       [-0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.],
       [-0.,  0.,  0.,  1., -0.],
       [-0.,  0.,  0.,  0.,  1.]])

In [236]: q, r = qr(mat) 

In [237]: r
Out[237]: 
array([[-1.6914,  4.38  ,  0.1757,  0.4075, -0.7838],
       [ 0.    , -2.6436,  0.1939, -3.072 , -1.0702],
       [ 0.    ,  0.    , -0.8138,  1.5414,  0.6155],
       [ 0.    ,  0.    ,  0.    , -2.6445, -2.1669],
       [ 0.    ,  0.    ,  0.    ,  0.    ,  0.0002]])
linalg矩阵分解

4.6 伪随机数生成

numpy.random模块提供了很多生成随机数的函数,可以选择生成符合某种概率分布的随机数。比如我们可以用normal得到一个4 x 4的,符合标准正态分布的数组:

import numpy as np
In [238]: samples = np.random.normal(size=(4, 4)) #符合标准正态分布的数组

In [239]: samples
Out[239]: 
array([[ 0.5732,  0.1933,  0.4429,  1.2796],
       [ 0.575 ,  0.4339, -0.7658, -1.237 ],
       [-0.5367,  1.8545, -0.92  , -0.1082],
       [ 0.1525,  0.9435, -1.0953, -0.144 ]])

相对的,python内建的random模块一次只能生成一个样本。在生成大量样本方法,numpy.random是非常快的:

In [240]: from random import normalvariate
#.normalvariate(mu,sigma) 正态分布的随机数
In [241]: N = 1000000

In [242]: %timeit samples = [normalvariate(0, 1) for _ in range(N)]
1.77 s +- 126 ms per loop (mean +- std. dev. of 7 runs, 1 loop each)
#用魔术函数%timeit。会自动多次执行以产生一个非常精确的平均执行时间。
In [243]: %timeit np.random.normal(size=N)
61.7 ms +- 1.32 ms per loop (mean +- std. dev. of 7 runs, 10 loops each)

之所以称之为伪随机数,是因为随机数生成算法是根据seed来生成的。也就是说,只要seed设置一样,每次生成的随机数是相同的:

In [244]: np.random.seed(1234)

当然,这个seed是全局的,如果想要避免全局状态,可以用numpy.random.RandomState来创建一个独立的生成器:

In [245]: rng = np.random.RandomState(1234)
#RandomState 设置成种子
In [246]: rng.randn(10)
Out[246]: 
array([ 0.4714, -1.191 ,  1.4327, -0.3127, -0.7206,  0.8872,  0.8596,
       -0.6365,  0.0157, -2.2427])
numpy.random中的部分函数

4.7 示例:随机漫步

这个例子让我了解一个在实际任务中如何利用数组操作。首先一个最简单的随机漫步:从0开始,步幅为1和-1,以相同的概率出现。

下面是纯python的实现方法,1000步:

In [247]: import random
   .....: position = 0
   .....: walk = [position]
   .....: steps = 1000
   .....: for i in range(steps):
   .....:     step = 1 if random.randint(0, 1) else -1 
#在python中的random.randint(a,b)用于生成一个指定范围内的整数。
   .....:     position += step
   .....:     walk.append(position)
   .....:
import matplotlib.pyplot as plt
%matplotlib inline
In [249]: plt.plot(walk[:100])
简单的随机漫步

随机漫步其实就是一个简单的累加。而用np.random能更快:

In [251]: nsteps = 1000

In [252]: draws = np.random.randint(0, 2, size=nsteps)

In [253]: steps = np.where(draws > 0, 1, -1)

In [254]: walk = steps.cumsum()

有了这些数据之后,我们就可以沿着漫步路径做一些统计工作了,比如求取最大值和最小值:

In [255]: walk.min()
Out[255]: -3

In [256]: walk.max()
Out[256]: 31

一个更复杂的统计值是在哪一步random walk到达了一个指定值。我们想知道从0走出10步用了多久,不论是正方向还是负方向。np.abs(walk) >= 10给我们一个布尔数组告诉我们是否超过10,但我们想要第一次出现的10或-10。因此,我们利用argmax来计算,这个会返回布尔数组中最大值的索引(Ture是最大值):

In [257]: (np.abs(walk) >= 10).argmax()  # argmax(f(x))是使得 f(x)取得最大值所对应的变量x
Out[257]: 37

注意,使用argmax并不总是效率的,因为它总会搜索整个数组。在这里例子里,一旦True被找到了,我们就返回为最大值。

Simulating Many Random Walks at Once(一次模拟多个随机漫步)

假设我们一次要模拟5000个随机漫步。传入一个2-tuple,np.random会生成一个二维数组,然后我们沿着每行来计算累加,这样就能一次模拟5000个:

##在python中的random.randint(a,b)用于生成一个指定范围内的整数。
In [258]: nwalks = 5000

In [259]: nsteps = 1000

In [260]: draws = np.random.randint(0, 2, size=(nwalks, nsteps)) # 0 or 1

In [261]: steps = np.where(draws > 0, 1, -1)

In [262]: walks = steps.cumsum(1)

In [263]: walks
Out[263]: 
array([[  1,   0,   1, ...,   8,   7,   8],
       [  1,   0,  -1, ...,  34,  33,  32],
       [  1,   0,  -1, ...,   4,   5,   4],
       ..., 
       [  1,   2,   1, ...,  24,  25,  26],
       [  1,   2,   3, ...,  14,  13,  14],
       [ -1,  -2,  -3, ..., -24, -23, -22]])

现在,我们来计算所有随机漫步过程的最大值和最小值:

In [264]: walks.max()
Out[264]: 138

In [265]: walks.min()
Out[265]: -133

得到这些数据之后,我们来计算30或-30的最小穿越时间。这里稍微复杂些,因为不是5000个过程都到达了30。我们可以用any方法来对此进行检查:

In [266]: hits30 = (np.abs(walks) >= 30).any(1)

In [267]: hits30
Out[267]: array([False,  True, False, ..., False,  True, False], dtype=bool)

In [268]: hits30.sum() # Number that hit 30 or -30
Out[268]: 3410

然后我们利用这个布尔型数组选出那些穿越了30(绝对值)的随机漫步(行),并调用argmax在轴1上获取穿越时间:

In [269]: crossing_times = (np.abs(walks[hits30]) >= 30).argmax(1)

In [270]: crossing_times.mean()
Out[270]: 498.88973607038122

请尝试用其他分布方式得到漫步数据。只需使用不同的随机数生成函数即可,如normal用于生成指定均值和标准差的正态分布数据:

In [271]: steps = np.random.normal(loc=0, scale=0.25,
   .....:                          size=(nwalks, nsteps))
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 199,902评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,037评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 146,978评论 0 332
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,867评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,763评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,104评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,565评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,236评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,379评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,313评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,363评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,034评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,637评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,719评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,952评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,371评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,948评论 2 341

推荐阅读更多精彩内容