《Python数据分析基础教程：NumPy学习指南（第2版）》笔记6：第三章 常用函数2——中位数、方差、日期、展平

本章将介绍NumPy的常用函数。具体来说，我们将以分析历史股价为例，介绍怎样从文件中载入数据，以及怎样使用NumPy的基本数学和统计分析函数。这里还将学习读写文件的方法，并尝试函数式编程和NumPy线性代数运算。

第三章 常用函数

3.9 统计分析

股票交易者对于收盘价的预测很感兴趣。常识告诉我们，这个价格应该接近于某种均值。算数平均值和加权平均值都是在数值分布中寻找中心点的方法。然而，它们对于异常值（outlier）既不鲁棒也不敏感。举例来说，如果我们有一个高达100万美元的收盘价，这将影响到我们的计算结果。

3.10 动手实践：简单统计分析

我们可以用一些阈值来除去异常值，但其实有更好的方法，那就是中位数。将各个变量值按大小顺序排列起来，形成一个数列，居于数列中间位置的那个数即为中位数。例如，我们有1、 2、3、 4、 5这5个数值，那么中位数就是中间的数字3。下面是计算中位数的步骤。

• (1) 计算收盘价的中位数。创建一个新的Python脚本文件，命名为simplestats.py。你已经知道如何从CSV文件中读取数据到数组中了，因此只需要复制一行代码并确保只获取收盘价数据即可，如下所示：

c=np.loadtxt('data.csv', delimiter=',', usecols=(6,), unpack=True)

• (2) median函数将帮助我们找到中位数。我们调用它并立即打印出结果。添加下面这行代码：

print("median =", np.median(c))

这段代码的输出内容如下：

median = 352.055

• (3) 既然这是我们首次使用median函数，我们来检查一下结果是否正确。这可不是因为我们多疑！当然，我们可以将整个数据文件浏览一遍并人工找到正确的答案，但那样太无趣了。我们将对价格数组进行排序，并输出排序后居于中间位置的值，这也就是模拟了寻找中位数的算法。 msort函数可以帮我们完成第一步。我们将调用这个函数，获得排序后的数组，并输出结果。

sorted_close = np.msort(c)
print("sorted =", sorted_close)

这段代码的输出内容如下：

sorted = [336.1  338.61 339.32 342.62 342.88 343.44 344.32 345.03 346.5  346.67
 348.16 349.31 350.56 351.88 351.99 352.12 352.47 353.21 354.54 355.2
 355.36 355.76 356.85 358.16 358.3  359.18 359.56 359.9  360.   363.13]

太好了，代码生效了！现在，我们来获取位于中间的那个数字：

N = len(c)
print("middle =", sorted_close[int((N - 1)/2)])

输出如下：

middle = 351.99

• (4) 咦，这个值和median函数给出的值不一样，这是怎么回事？经过仔细观察我们发现，median函数返回的结果甚至根本没有在我们的数据文件里出现过。这就更奇怪了！在给NumPy团队提交bug报告之前，我们先来看下文档。原来这个谜团很容易解开。原因就在于我们的简单算法模拟只对长度为奇数的数组奏效。对于长度为偶数的数组，中位数的值应该等于中间那两个数的平均值。因此，输入如下代码：

print("average middle =", (sorted_close[int(N /2)] + sorted_close[int((N - 1) / 2)]) / 2)

输出结果如下：

average middle = 352.055

成功了！

• (5) 另外一个我们关心的统计量就是方差。方差能够体现变量变化的程度。在我们的例子中，方差还可以告诉我们投资风险的大小。那些股价变动过于剧烈的股票一定会给持有者制造麻烦。
  在NumPy中，计算方差只需要一行代码，看下面：

print "variance =", np.var(c)

将给出如下结果：

variance = 50.126517888888884

• (6) 既然我们不相信NumPy的函数，那就再次根据文档中方差的定义来复核一下结果。注意，这里方差的定义可能与你在统计学的书中看到的不一致，但这个定义在统计学上更为通用。方差是指各个数据与所有数据算术平均数的离差平方和除以数据个数所得到的值。
  一些书里面告诉我们，应该用数据个数减1去除离差平方和¹。

print("variance from definition =", np.mean((c - c.mean())**2))

输出结果如下：

variance from definition = 50.126517888888884

• simplestats.py完整代码如下：

import numpy as np

c = np.loadtxt('data.csv', delimiter=',', usecols=(6,), unpack=True)

print("median =", np.median(c))
sorted_close = np.msort(c)
print("sorted =", sorted_close)

N = len(c)
print("middle =", sorted_close[int((N - 1)/2)])
print("average middle =", (sorted_close[int(N /2)] + sorted_close[int((N - 1) / 2)]) / 2)

print("variance =", np.var(c))
print("variance from definition =", np.mean((c - c.mean())**2))

3.11 股票收益率

在学术文献中，收盘价的分析常常是基于股票收益率和对数收益率的。简单收益率是指相邻两个价格之间的变化率，而对数收益率是指所有价格取对数后两两之间的差值。我们在高中学习过对数的知识，“a”的对数减去“b”的对数就等于“a除以b”的对数。因此，对数收益率也可以用来衡量价格的变化率。注意，由于收益率是一个比值，例如我们用美元除以美元（也可以是其他货币单位），因此它是无量纲的。总之，投资者最感兴趣的是收益率的方差或标准差，因为这代表着投资风险的大小。

3.12 动手实践：分析股票收益率

按照如下步骤分析股票收益率。

• (1) 首先，我们来计算简单收益率。 NumPy中的diff函数可以返回一个由相邻数组元素的差值构成的数组。这有点类似于微积分中的微分。为了计算收益率，我们还需要用差值除以前一天的价格。不过这里要注意， diff返回的数组比收盘价数组少一个元素。经过仔细考虑，我们使用如下代码：

returns = np.diff( c ) / c[ : -1]

注意，我们没有用收盘价数组中的最后一个值做除数。接下来，用std函数计算标准差：

print("Standard deviation =", np.std(returns))

输出结果如下：

Standard deviation = 0.012922134436826306

• (2) 对数收益率计算起来甚至更简单一些。我们先用log函数得到每一个收盘价的对数，再对结果使用diff函数即可。

logreturns = np.diff( np.log(c) )

一般情况下，我们应检查输入数组以确保其不含有零和负数。否则，将得到一个错误提示。
不过在我们的例子中，股价总为正值，所以可以将检查省略掉。

• (3) 我们很可能对哪些交易日的收益率为正值非常感兴趣。在完成了前面的步骤之后，我们只需要用where函数就可以做到这一点。 where函数可以根据指定的条件返回所有满足条件的数组元素的索引值。输入如下代码：

posretindices = np.where(returns > 0)
print("Indices with positive returns", posretindices)

即可输出该数组中所有正值元素的索引。

Indices with positive returns (array([ 0,  1,  4,  5,  6,  7,  9, 10, 11, 12, 16, 17, 18, 19, 21, 22, 23,
       25, 28], dtype=int64),)

• (4) 在投资学中，波动率（volatility）是对价格变动的一种度量。历史波动率可以根据历史价格数据计算得出。计算历史波动率（如年波动率或月波动率）时，需要用到对数收益率。**年波动率等于对数收益率的标准差除以其均值，再除以交易日倒数的平方根，通常交易日取252天。**我们用std和mean函数来计算，代码如下所示：

annual_volatility = np.std(logreturns)/np.mean(logreturns)
annual_volatility = annual_volatility / np.sqrt(1./252.)
print("Yearly volatility",annual_volatility)

年波动率输出为：

Yearly volatility 129.27478991115132

• (5) 请注意sqrt函数中的除法运算。在Python中，整数的除法和浮点数的除法运算机制不同，我们必须使用浮点数才能得到正确的结果。与计算年波动率的方法类似，计算月波动率如下：

print("Monthly volatility", annual_volatility * np.sqrt(1./12.))

月波动率输出为：

Monthly volatility 37.318417377317765

分析股票收益率示例完整代码return.py文件代码如下：

import numpy as np

c = np.loadtxt('data.csv', delimiter=',', usecols=(6,), unpack=True)

returns = np.diff( c ) / c[ : -1]
print("Standard deviation =", np.std(returns))

logreturns = np.diff( np.log(c) )
posretindices = np.where(returns > 0)
print("Indices with positive returns", posretindices)

annual_volatility = np.std(logreturns)/np.mean(logreturns)
annual_volatility = annual_volatility / np.sqrt(1./252.)
print("Yearly volatility", annual_volatility)
print("Monthly volatility", annual_volatility * np.sqrt(1./12.))

3.13 日期分析

你是否有时候会有星期一焦虑症和星期五狂躁症？想知道股票市场是否受上述现象的影响？我认为这值得深入研究。

3.14 动手实践：分析日期数据

首先，我们要读入收盘价数据。随后，根据星期几来切分收盘价数据，并分别计算平均价格。

最后，我们将找出一周中哪一天的平均收盘价最高，哪一天的最低。在我们动手之前，有一个善意的提醒：你可能希望利用分析结果在某一天买股票或卖股票，然而我们这里的数据量不足以做出可靠的决策，请先咨询专业的统计分析师再做决定！

程序员不喜欢日期，因为处理日期总是很烦琐。 NumPy是面向浮点数运算的，因此需要对日期做一些专门的处理。请自行尝试如下代码，单独编写脚本文件或使用本书附带的代码文件：

import numpy as np
from datetime import datetime

dates, close=np.loadtxt('data.csv', delimiter=',', usecols=(1,6), unpack=True)

执行以上代码后会得到一个错误提示：

ValueError: invalid literal for float(): 28-01-2011

按如下步骤处理日期。

• (1) 显然， NumPy尝试把日期转换成浮点数。我们需要做的是显式地告诉NumPy怎样来转换日期，而这需要用到loadtxt函数中的一个特定的参数。这个参数就是converters，它是一本数据列和转换函数之间进行映射的字典。

为此，我们必须写出转换函数：

def datestr2num(s):
    return datetime.strptime(s.decode('ascii'), "%d-%m-%Y").date().weekday()

我们将日期作为字符串传给datestr2num函数，如“28-01-2011”。这个字符串首先会按照指定的形式"%d-%m-%Y"转换成一个datetime对象。补充一点，这是由Python标准库提供的功能，与NumPy无关。随后， datetime对象被转换为date对象。最后，调用weekday方法返回一个数字。如同你在注释中看到的，这个数字可以是0到6的整数， 0代表星期一， 6代表星期天。当然，具体的数字并不重要，只是用作标识。

注意：data.csv文件第二列数据为日期格式字符串，loadtxt()默认读取数据需要为二进制编码格式，而返回的值为字节字符串bytes，所以需要把它转化二进制格式的string，故需要对字符串解码，使用函数decode(‘asii’)，变成string格式。否则，抛出TypeError: strptime() argument 1 must be str, not bytes异常。

• (2) 接下来，我们将日期转换函数挂接上去，这样就可以读入数据了。

dates, close=np.loadtxt('data.csv', delimiter=',', usecols=(1,6), converters={1: datestr2num}, unpack=True)
print("Dates =", dates)

输出结果如下：

Dates = [4. 0. 1. 2. 3. 4. 0. 1. 2. 3. 4. 0. 1. 2. 3. 4. 1. 2. 3. 4. 0. 1. 2. 3.
 4. 0. 1. 2. 3. 4.]

如你所见，没有出现星期六和星期天。股票交易在周末是休市的。

• (3) 我们来创建一个包含5个元素的数组，分别代表一周的5个工作日。数组元素将初始化为0。

averages = np.zeros(5)

这个数组将用于保存各工作日的平均收盘价。

• (4) 我们已经知道， where函数会根据指定的条件返回所有满足条件的数组元素的索引值。
  **take函数可以按照这些索引值从数组中取出相应的元素。**我们将用take函数来获取各个工作日的收盘价。在下面的循环体中，我们将遍历0到4的日期标识，或者说是遍历星期一到星期五，然后用where函数得到各工作日的索引值并存储在indices数组中。在用take函数获取这些索引值相应的元素值。最后，我们对每个工作日计算出平均值存放在averages数组中。代码如下：

averages = np.zeros(5)

for i in range(5):
    indices = np.where(dates == i) 
    prices = np.take(close, indices)
    avg = np.mean(prices)
    print("Day", i, "prices", prices, "Average", avg)
    averages[i] = avg

输出结果如下：

Day 0 prices [[339.32 351.88 359.18 353.21 355.36]] Average 351.7900000000001
Day 1 prices [[345.03 355.2  359.9  338.61 349.31 355.76]] Average 350.63500000000005
Day 2 prices [[344.32 358.16 363.13 342.62 352.12 352.47]] Average 352.1366666666666
Day 3 prices [[343.44 354.54 358.3  342.88 359.56 346.67]] Average 350.8983333333333
Day 4 prices [[336.1  346.5  356.85 350.56 348.16 360.   351.99]] Average 350.0228571428571

• (5) 如果你愿意，还可以找出哪个工作日的平均收盘价是最高的，哪个是最低的。这很容易做到，用max和min函数即可，代码如下：

top = np.max(averages)
print( "Highest average", top)
print( "Top day of the week", np.argmax(averages))

bottom = np.min(averages)
print( "Lowest average", bottom)
print( "Bottom day of the week", np.argmin(averages))

输出结果如下：

Highest average 352.1366666666666
Top day of the week 2
Lowest average 350.0228571428571
Bottom day of the week 4

刚才做了些什么
argmin函数返回的是averages数组中最小元素的索引值，这里是4，也就是星期五。而
argmax函数返回的是averages数组中最大元素的索引值，这里是2，也就是星期三。

示例代码完整代码如下：

import numpy as np
from datetime import datetime

def datestr2num(s):
    return datetime.strptime(s.decode('ascii'), "%d-%m-%Y").date().weekday()

dates, close=np.loadtxt('data.csv', delimiter=',', usecols=(1,6), converters={1: datestr2num}, unpack=True)
print("Dates =", dates)

averages = np.zeros(5)

for i in range(5):
    indices = np.where(dates == i) 
    prices = np.take(close, indices)
    avg = np.mean(prices)
    print("Day", i, "prices", prices, "Average", avg)
    averages[i] = avg
    
top = np.max(averages)
print( "Highest average", top)
print( "Top day of the week", np.argmax(averages))

bottom = np.min(averages)
print( "Lowest average", bottom)
print( "Bottom day of the week", np.argmin(averages))

3.15 周汇总

在之前的“动手实践”教程中，我们用的是盘后数据。也就是说，这些数据是将一整天的交易数据汇总得到的。如果你对棉花市场感兴趣，并且有数十年的数据，你可能希望对数据做进一步的汇总和压缩。开始动手吧。我们来把苹果股票数据按周进行汇总。

3.16 动手实践：汇总数据

我们将要汇总整个交易周中从周一到周五的所有数据。数据覆盖的时间段内有一个节假日：2月21日是总统纪念日。这天是星期一，美国股市休市，因此在我们的示例数据中没有这一天的数据记录。数据中的第一天为星期五，处理起来不太方便。按照如下步骤来汇总数据。

• (1) 为了简单起见，我们只考虑前三周的数据，这样就避免了节假日造成的数据缺失。你可以稍后尝试对其进行拓展。

import numpy as np
from datetime import datetime

def datestr2num(s):
    return datetime.strptime(s.decode('ascii'), "%d-%m-%Y").date().weekday()

dates, open, high, low, close=np.loadtxt('data.csv', delimiter=',', usecols=(1, 3, 4, 5, 6), converters={1: datestr2num}, unpack=True)

close = close[:16]
dates = dates[:16]

• (2) 首先我们来找到示例数据中的第一个星期一。回忆一下，在Python中星期一对应的编码是0，这可以作为where函数的条件。接着，我们要取出数组中的首个元素，其索引值为0。但where函数返回的结果是一个多维数组，因此要用ravel函数将其展平。

# 找到第一个星期一
first_monday = np.ravel(np.where(dates == 0))[0]
print( "The first Monday index is", first_monday)

输出结果如下：

The first Monday index is 1

• (3) 下面要做的是找到示例数据的最后一个星期五，方法和找第一个星期一类似。星期五相对应的编码是4。此外，我们用-1作为索引值来定位数组的最后一个元素。

# 找到最后一个星期五
last_friday = np.ravel(np.where(dates == 4))[-1]
print( "The last Friday index is", last_friday)

输出结果如下：

The last Friday index is 15

接下来创建一个数组，用于存储三周内每一天的索引值。

weeks_indices = np.arange(first_monday, last_friday + 1)
print( "Weeks indices initial", weeks_indices)

输出结果如下：

Weeks indices initial [ 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15]

• (4) 按照每个子数组5个元素，用split函数切分数组：

weeks_indices = np.split(weeks_indices, 3)
print( "Weeks indices after split", weeks_indices)

输出结果如下：

Weeks indices after split [array([1, 2, 3, 4, 5], dtype=int64), array([ 6,  7,  8,  9, 10], dtype=int64), array([11, 12, 13, 14, 15], dtype=int64)]

• (5) 在NumPy中，数组的维度也被称作轴。现在我们来熟悉一下apply_along_axis函数。这个函数会调用另外一个由我们给出的函数，作用于每一个数组元素上。目前我们的数组中有3个元素，分别对应于示例数据中的3个星期，元素中的索引值对应于示例数据中的1天。在调用apply_along_axis时提供我们自定义的函数名summarize，并指定要作用的轴或维度的编号（如取1）、目标数组以及可变数量的summarize函数的参数。

weeksummary = np.apply_along_axis(summarize, 1, weeks_indices, open, high, low, close)
print( "Week summary", weeksummary)

(6) 编写summarize函数。该函数将为每一周的数据返回一个元组，包含这一周的开盘价、最高价、最低价和收盘价，类似于每天的盘后数据。

def summarize(a, o, h, l, c):
    monday_open = o[a[0]]
    week_high = np.max( np.take(h, a) )
    week_low = np.min( np.take(l, a) )
    friday_close = c[a[-1]]

    return("APPL", monday_open, week_high, week_low, friday_close)

注意，我们用take函数来根据索引值获取数组元素的值，并用max和min函数轻松计算出一周的最高股价和最低股价。一周的开盘价即为周一的开盘价，而一周的收盘价即为周五的收盘价。

Week summary [['APPL' '335.8' '346.7' '334.3' '346.5']
 ['APPL' '347.89' '360.0' '347.64' '356.85']
 ['APPL' '356.79' '364.9' '349.52' '350.56']]

(7) 使用NumPy中的savetxt函数，将数据保存至文件。

np.savetxt("weeksummary.csv", weeksummary, delimiter=",", fmt="%s")

如代码中所示，我们指定了文件名、需要保存的数组名、分隔符（在这个例子中为英文标点逗号）以及存储浮点数的格式。
格式字符串以一个百分号开始。接下来是一个可选的标志字符： -表示结果左对齐， 0表示左端补0， +表示输出符号（正号+或负号-）。第三部分为可选的输出宽度参数，表示输出的最小位数。第四部分是精度格式符，以"."开头，后面跟一个表示精度的整数。最后是一个类型指定字符，在我们的例子中指定为字符串类型。

weeksummary.csv文件内容如下：

APPL,335.8,346.7,334.3,346.5
APPL,347.89,360.0,347.64,356.85
APPL,356.79,364.9,349.52,350.56

示例代码完整代码如下：

import numpy as np
from datetime import datetime

def datestr2num(s):
    return datetime.strptime(s.decode('ascii'), "%d-%m-%Y").date().weekday()

dates, open, high, low, close=np.loadtxt('data.csv', delimiter=',', usecols=(1, 3, 4, 5, 6), converters={1: datestr2num}, unpack=True)

close = close[:16]
dates = dates[:16]

# get first Monday
first_monday = np.ravel(np.where(dates == 0))[0]
print( "The first Monday index is", first_monday)

# get last Friday
last_friday = np.ravel(np.where(dates == 4))[-1]
print( "The last Friday index is", last_friday)

weeks_indices = np.arange(first_monday, last_friday + 1)
print( "Weeks indices initial", weeks_indices)

weeks_indices = np.split(weeks_indices, 3)
print( "Weeks indices after split", weeks_indices)

def summarize(a, o, h, l, c):
    monday_open = o[a[0]]
    week_high = np.max( np.take(h, a) )
    week_low = np.min( np.take(l, a) )
    friday_close = c[a[-1]]

    return("APPL", monday_open, week_high, week_low, friday_close)

weeksummary = np.apply_along_axis(summarize, 1, weeks_indices, open, high, low, close)
print( "Week summary", weeksummary)

np.savetxt("weeksummary.csv", weeksummary, delimiter=",", fmt="%s")

---

1. 注意样本方差和总体方差在计算上的区别。总体方差是用数据个数去除离差平方和，而样本方差则是用样本数据个数减1去除离差平方和，其中样本数据个数减1（即n-1）称为自由度。之所以有这样的差别，是为了保证样本方差是一个无偏估计量。——译者注 ↩︎