Python数据分析基础教程：Numpy学习指南

第二章 Numpy基础

2.6 改变数组维度

ravel()、flatten() 将多维数组展平

b.transpose() 矩阵转置，等同于b.T，一维数组不变

reshape() 改变数组维度

2.8 组合数组

hstack((a, b)) 水平组合，等同于 concatenate((a, b), axis=1)
vstack((a, b)) 垂直组合，等同于 concatenate((a, b), axis=0)

column_stack((a, b)) 列组合，二维等同于hstack
row_stack((a, b)) 行组合，二维等同与vstack

2.10 分割数组

In: a
Out:
array([[0, 1, 2],
      [3, 4, 5],
      [6, 7, 8]])
In:hsplit(a, 3)    #水平分割，等同于 split(a,3,axis=1)
Out:
[array([[0],
        [3],
        [6]]),
array ([[1],
        [4],
        [7]]),
array ([[2],
        [5],
        [8]])]

vsplit(a,3) 垂直分割 ，等同于 split(a,3,axis=0)

2.11 数组属性

ndim 数组维度，或数组轴的个数

size 数组元素总数

itemsize 数组元素在内存中所占字节数

nbytes 数组所占存储空间 = itemsize * size

b = array([1.j + 1, 2.j + 3]) 虚数
b.real 复数数组的实部 b.imag 虚部

flat属性将返回一个numpy.flatiter对象，可以让我们像遍历一维数组一样去遍历任意的多维数组。

In: b = arange(4).reshape(2,2)
In: b
Out:
array([[0, 1],
      [2, 3]])
In: f = b.flat
In: f
Out: <numpy.flatiter object at 0x103013e00>
In: for item in f: print item
.....:

2.12 数组转换

tolist函数将NumPy数组转换成Python列表。

In: b
Out: array([ 1.+1.j, 3.+2.j])
In: b.tolist()
Out: [(1+1j), (3+2j)]

astype函数可以在转换数组时指定数据类型。

In: b
Out: array([ 1.+1.j, 3.+2.j])
In: b.astype(int)
/usr/local/bin/ipython:1: ComplexWarning: Casting complex values to real discards the imaginary part #虚部丢失，转为b.astype('complex') 则不会发生错误。

#!/usr/bin/python
Out: array([1, 3])

3.2 读写文件

savetxt

import numpy as np
i2 = np.eye(2)
np.savetxt("eye.txt", i2)

3.4 读入CSV文件

# AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800

c,v=np.loadtxt('data.csv', delimiter=',', usecols=(6,7), unpack=True) #index从0开始

3.6.1 算术平均值

np.mean(c) = np.average(c)

3.6.2 加权平均值

t = np.arange(len(c))
np.average(c, weights=t)

3.8 极值

np.min(c)
np.max(c)

np.ptp(c) 最大值与最小值的差值

3.10 统计分析

np.median(c) 中位数
np.msort(c) 升序排序
np.var(c) 方差

3.12 分析股票收益率

np.diff(c) 可以返回一个由相邻数组元素的差
值构成的数组

returns = np.diff( arr ) / arr[ : -1]  #diff返回的数组比收盘价数组少一个元素

np.std(c) 标准差

对数收益率

logreturns = np.diff( np.log(c) )  #应检查输入数组以确保其不含有零和负数

where 可以根据指定的条件返回所有满足条件的数
组元素的索引值。
posretindices = np.where(returns > 0)

np.sqrt(1./252.) 平方根，浮点数

3.14 分析日期数据

# AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800

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

def datestr2num(s):
    return datetime.datetime.strptime(s, "%d-%m-%Y").date().weekday()
# 星期一 0
# 星期二 1
# 星期三 2
# 星期四 3
# 星期五 4
# 星期六 5
# 星期日 6

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

averages = np.zeros(5)
for i in range(5):
    indices = np.where(dates == i)
    prices = np.take(close, indices)  #按数组的元素运算,产生一个数组作为输出。

>>> a = [4, 3, 5, 7, 6, 8]
>>> indices = [0, 1, 4]
>>> np.take(a, indices)
array([4, 3, 6])

np.argmax(c) #返回的是数组中最大元素的索引值
np.argmin(c)

3.16 汇总数据

# AAPL,28-01-2011, ,344.17,344.4,333.53,336.1,21144800

#得到第一个星期一和最后一个星期五
first_monday = np.ravel(np.where(dates == 0))[0]
last_friday = np.ravel(np.where(dates == 4))[-1]

#创建一个数组，用于存储三周内每一天的索引值
weeks_indices = np.arange(first_monday, last_friday + 1)
#按照每个子数组5个元素，用split函数切分数组
weeks_indices = np.split(weeks_indices, 5)
#output
[array([1, 2, 3, 4, 5]), array([ 6, 7, 8, 9, 10]), array([11,12, 13, 14, 15])]

weeksummary = np.apply_along_axis(summarize, 1, weeks_indices,open, high, low, close)

def summarize(a, o, h, l, c): #open, high, low, close
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)

np.savetxt("weeksummary.csv", weeksummary, delimiter=",", fmt="%s") #指定了文件名、需要保存的数组名、分隔符（在这个例子中为英文标点逗号）以及存储浮点数的格式。

格式字符串以一个百分号开始。接下来是一个可选的标志字符：-表示结果左对齐，0表示左端补0，+表示输出符号（正号+或负号-）。第三部分为可选的输出宽度参数，表示输出的最小位数。第四部分是精度格式符，以”.”开头，后面跟一个表示精度的整数。最后是一个类型指定字符，在例子中指定为字符串类型。

numpy.apply_along_axis(func1d, axis, arr, *args, **kwargs)

>>> def my_func(a):
...     """Average first and last element of a 1-D array"""
...     return (a[0] + a[-1]) * 0.5
>>> b = np.array([[1,2,3], [4,5,6], [7,8,9]])
>>> np.apply_along_axis(my_func, 0, b)  #沿着X轴运动，取列切片
array([ 4.,  5.,  6.])
>>> np.apply_along_axis(my_func, 1, b)  #沿着y轴运动，取行切片
array([ 2.,  5.,  8.])


>>> b = np.array([[8,1,7], [4,3,9], [5,2,6]])
>>> np.apply_along_axis(sorted, 1, b)
array([[1, 7, 8],
       [3, 4, 9],
       [2, 5, 6]])

3.20 计算简单移动平均线

(1) 使用ones函数创建一个长度为N的元素均初始化为1的数组，然后对整个数组除以N，即可得到权重。如下所示：

N = int(sys.argv[1])
weights = np.ones(N) / N
print "Weights", weights

在N = 5时，输出结果如下：

Weights [ 0.2 0.2 0.2 0.2 0.2]  #权重相等

(2) 使用这些权重值，调用convolve函数：

c = np.loadtxt('data.csv', delimiter=',', usecols=(6,),unpack=True)
sma = np.convolve(weights, c)[N-1:-N+1]   #卷积是分析数学中一种重要的运算，定义为一个函数与经过翻转和平移的另一个函数的乘积的积分。

t = np.arange(N - 1, len(c))   #作图
plot(t, c[N-1:], lw=1.0)
plot(t, sma, lw=2.0)
show()

3.22 计算指数移动平均线

指数移动平均线（exponential moving average）。指数移动平均线使用的权重是指数衰减的。对历史上的数据点赋予的权重以指数速度减小，但永远不会到达0。

x = np.arange(5)
print "Exp", np.exp(x)
#output
Exp [ 1. 2.71828183 7.3890561 20.08553692 54.59815003]

Linspace 返回一个元素值在指定的范围内均匀分布的数组。

print "Linspace", np.linspace(-1, 0, 5) #起始值、终止值、可选的元素个数
#output 
Linspace [-1. -0.75 -0.5 -0.25 0. ]

(1)权重计算

N = int(sys.argv[1])
weights = np.exp(np.linspace(-1. , 0. , N))

(2)权重归一化处理

weights /= weights.sum()
print "Weights", weights
#output
Weights [ 0.11405072 0.14644403 0.18803785 0.24144538 0.31002201]

(3)计算及作图

c = np.loadtxt('data.csv', delimiter=',', usecols=(6,),unpack=True)
ema = np.convolve(weights, c)[N-1:-N+1]
t = np.arange(N - 1, len(c))
plot(t, c[N-1:], lw=1.0)
plot(t, ema, lw=2.0)
show()

3.26 用线性模型预测价格

(x, residuals, rank, s) = np.linalg.lstsq(A, b) #系数向量x、一个残差数组、A的秩以及A的奇异值
print x, residuals, rank, s
#计算下一个预测值
print np.dot(b, x)

3.28 绘制趋势线

>>> x = np.arange(6)
>>> x = x.reshape((2, 3))
>>> x
array([[0, 1, 2],
       [3, 4, 5]])
>>> np.ones_like(x)   #用1填充数组
array([[1, 1, 1],
       [1, 1, 1]])

类似函数
zeros_like
empty_like
zeros
ones
empty

3.30 数组的修剪和压缩

a = np.arange(5)
print "a =", a
print "Clipped", a.clip(1, 2) #将所有比给定最大值还大的元素全部设为给定的最大值，而所有比给定最小值还小的元素全部设为给定的最小值
#output
a = [0 1 2 3 4]
Clipped [1 1 2 2 2]

a = np.arange(4)
print a
print "Compressed", a.compress(a > 2) #返回一个根据给定条件筛选后的数组
#output
[0 1 2 3]
Compressed [3]

b = np.arange(1, 9)
print "b =", b
print "Factorial", b.prod() #输出数组元素阶乘结果
#output
b = [1 2 3 4 5 6 7 8]
Factorial 40320

print "Factorials", b.cumprod()
#output
Factorials [ 1 2 6 24 120 720 5040 40320] #数组元素遍历阶乘

4.2 股票相关性分析

covariance = np.cov(a，b)

详解协方差与协方差矩阵

获取对角元素
covariance.diagonal()

4.4 多项式拟合

bhp=np.loadtxt('BHP.csv', delimiter=',', usecols=(6,), unpack=True)
vale=np.loadtxt('VALE.csv', delimiter=',', usecols=(6,),unpack=True)
t = np.arange(len(bhp))
poly = np.polyfit(t, bhp - vale, int(sys.argv[1])) # sys.argv[1]为3,即用3阶多项式拟合数据
print "Polynomial fit", poly

#output
Polynomial fit [ 1.11655581e-03 -5.28581762e-02 5.80684638e-01 5.79791202e+01]
#预测下个值
print "Next value", np.polyval(poly, t[-1] + 1)

使用polyder函数对多项式函数求导(以求极值)

der = np.polyder(poly)
print "Derivative", der
#output
Derivative [ 0.00334967 -0.10571635 0.58068464]

求出导数函数的根，即找出原多项式函数的极值点

print "Extremas", np.roots(der)
#output
Extremas [ 24.47820054 7.08205278]

注：书中提示，3阶多项式拟合数据的结果并不好，可尝试更高阶的多项式拟合。

4.6 计算OBV（On-Balance Volume）净额成交量

diff函数可以计算数组中两个连续元素的差值，并返回一个由这些差值组成的数组。
change = np.diff(c)

sign函数可以返回数组中每个元素的正负符号，数组元素为负时返回-1，为正时返回1，否则返回0
np.sign(change)

使用piecewise（分段的）函数来获取数组元素的正负。使用合适的返回值和对应的条件调用该函数：

pieces = np.piecewise(change, [change < 0, change > 0], [-1, 1])
print "Pieces", pieces

检查一致性
np.array_equal(a, b)

np.vectorize 替代循环

>>> def myfunc(a, b):
...     "Return a-b if a>b, otherwise return a+b"
...     if a > b:
...         return a - b
...     else:
...         return a + b

>>> vfunc = np.vectorize(myfunc)
>>> vfunc([1, 2, 3, 4], 2)
array([3, 4, 1, 2])

The vectorize function is provided primarily for convenience, not for performance. The implementation is essentially a for loop.

4.10 使用hanning 函数平滑数据

(1) 调用hanning函数计算权重，生成一个长度为N的窗口（在这个示例中N取8）

N = int(sys.argv[1])
weights = np.hanning(N)
print "Weights", weights

#output
Weights [ 0. 0.1882551 0.61126047 0.95048443 0.95048443 0.61126047 0.1882551 0. ]

bhp = np.loadtxt('BHP.csv', delimiter=',', usecols=(6,),unpack=True)   #某股票数据
bhp_returns = np.diff(bhp) / bhp[ : -1] #股票收益率计算
smooth_bhp = np.convolve(weights/weights.sum(), bhp_returns) [N-1:-N+1]  #使用weights平滑股票收益率

#绘图
t = np.arange(N - 1, len(bhp_returns))
plot(t, bhp_returns[N-1:], lw=1.0)
plot(t, smooth_bhp, lw=2.0)
show()

两个多项式做差运算
poly_sub = np.polysub(a, b)

选择函数
numpy.select(condlist, choicelist, default=0)

>>> x = np.arange(10)
>>> condlist = [x<3, x>5]  
#输出两个array [true...false...],[false,...true]
>>> choicelist = [x, x**2]
>>> np.select(condlist, choicelist)  
array([ 0,  1,  2,  0,  0,  0, 36, 49, 64, 81])

trim_zeros函数可以去掉一维数组中开头和末尾为0的元素：
np.trim_zeros(a)

5.2 创建矩阵（略）

5.4 从已有矩阵创建新矩阵（略）

5.6 创建通用函数（略）

5.7 通用函数的方法（略）

5.8 在add上调用通用函数的方法（略）

5.10 数组的除法运算

在NumPy中，基本算术运算符+、-和*隐式关联着通用函数add、subtract和multiply。
也就是说，当你对NumPy数组使用这些算术运算符时，对应的通用函数将自动被调用。除法包含的过程则较为复杂，在数组的除法运算中涉及三个通用函数divide、true_divide和floor_division，以及两个对应的运算符/和//。

a = np.array([2, 6, 5])
b = np.array([1, 2, 3])
print "Divide", np.divide(a, b), np.divide(b, a)
#output
Divide [2 3 1] [0 0 0]

print "True Divide", np.true_divide(a, b), np.true_divide(b, a)
#output
True Divide [2. 3. 1.66666667] [0.5 0.33333333 0.6 ]



print "Floor Divide", np.floor_divide(a, b), np.floor_divide(b, a) c = 3.14 * b   #floor_divide函数总是返回整数结果，相当于先调用divide函数再调用floor函数。
print "Floor Divide 2", np.floor_divide(c, b), np.floor_divide(b, c) #floor函数将对浮点数进行向下取整并返回整数。

#output
Floor Divide [2 3 1] [0 0 0]
Floor Divide 2 [ 3. 3. 3.] [ 0. 0. 0.]

默认情况下，使用/运算符相当于调用divide函数
运算符//对应于floor_divide函数

print "// operator", a//b, b//a

• 
  

5.12 模运算

remainder函数逐个返回两个数组中元素相除后的余数。如果第二个数字为0，则直接返回0。

a = np.arange(-4, 4)
print "Remainder", np.remainder(a, 2) #等同于a % 2
#output
Remainder [0 1 0 1 0 1 0 1]

mod函数与remainder函数的功能完全一致

fmod函数处理负数的方式与remainder、mod和%不同。所得余数的正负由被除数决定，与除数的正负无关。

print "Fmod", np.fmod(a, 2)
#output
Fmod [ 0 -1 0 -1 0 1 0 1]

5.14 计算斐波那契数列

rint函数对浮点数取整，但结果仍为浮点数类型

a = np.arange(1,3,0.5) #[1.,1.5,2.,2.5]
pr  int np.rint(a)  #[1.,2.,2.,2.]

5.16 绘制利萨茹曲线（略）

5.18 绘制方波（略）

5.20 绘制锯齿波和三角波（略）

5.22 玩转二进制位

三个运用位操作的小技巧——检查两个整数的符号是否一致，检查一个数是否为2的幂数，以及计算一个数被2的幂数整除后的余数。

XOR或者^操作符。XOR操作符又被称为不等运算符，因此当两个操作数的符号不一致时，XOR操作的结果为负数。

import numpy as np
x = np.arange(-9, 9)
y = -x
print "Sign different?", (x ^ y) < 0 
print "Sign different?", np.less(np.bitwise_xor(x, y), 0) #^操作符对应于bitwise_xor函数，<操作符对应于less函数

在二进制数中，2的幂数表示为一个1后面跟一串0的形式，例如10、100、1000等。而比2的幂数小1的数表示为一串二进制的1，例如11、111、1111（即十进制里的3、7、15）等。如果我们在2的幂数以及比它小1的数之间执行位与操作AND，那么应该得到0。

print "Power of 2?\n", x, "\n", (x & (x - 1)) == 0
print "Power of 2?\n", x, "\n", np.equal(np.bitwise_and(x, (x - 1)), 0) #&操作符对应于bitwise_and函数，==操作符对应于equal函数。

计算余数的技巧实际上只在模为2的幂数（如4、8、16等）时有效。二进制的位左移一位，则数值翻倍。在前一个小技巧中我们看到，将2的幂数减去1可以得到一串1组成的二进制数，如11、111、1111等。这为我们提供了掩码（mask），与这样的掩码做位与操作AND即可得到以2的幂数作为模的余数。

print "Modulus 4\n", x, "\n", x & ((1 << 2) - 1)
print "Modulus 4\n", x, "\n", np.bitwise_and(x, np.left_shift(1, 2) - 1) #<<操作符对应于left_shift函数

6.2 计算逆矩阵（略）

6.3 求解线性方程

A = np.mat("1 -2 1;0 2 -8;-4 5 9")
b = np.array([0, 8, -9])
x = np.linalg.solve(A, b)

6.6 求解特征值和特征向量（略）

6.8 分解矩阵（略）

6.10 计算广义逆矩阵（略）

6.12 计算矩阵的行列式（略）

6.14 计算傅里叶变换（略）

6.16 移频（略）

6.18 硬币赌博游戏

对一个硬币赌博游戏下8份赌注。每一轮抛9枚硬币，如果少于5枚硬币正面朝上，将损失8份赌注中的1份；否则，将赢得1份赌注。模拟一下赌博的过程，初始资本为1000份赌注。

import numpy as np
from matplotlib.pyplot import plot, show
cash = np.zeros(10000) #将在赌场中玩10000轮硬币赌博游戏。
cash[0] = 1000
outcome = np.random.binomial(9, 0.5, size=len(cash)) #二项式分布
for i in range(1, len(cash)):
    if outcome[i] < 5:
        cash[i] = cash[i - 1] - 1
    elif outcome[i] < 10:
        cash[i] = cash[i - 1] + 1
    else:
        raise AssertionError("Unexpected outcome " + outcome)

print outcome.min(), outcome.max()
plot(np.arange(len(cash)), cash)
show() #剩余资本呈随机游走（random walk）状态

6.20 模拟游戏秀节目

超几何分布（hypergeometric distribution）是一种离散概率分布，它描述的是一个罐子里有两种物件，无放回地从中抽取指定数量的物件后，抽出指定种类物件的数量。

设想有这样一个游戏秀节目，每当参赛者回答对一个问题，他们可以从一个罐子里摸出3个球并放回。罐子里有一个“倒霉球”，一旦这个球被摸出，参赛者会被扣去6分。而如果他们摸出的3个球全部来自其余的25个普通球，那么可以得到1分。因此，如果一共有100道问题被正确回答，得分情况会是怎样的呢？

(1) 使用hypergeometric函数初始化游戏的结果矩阵。该函数的第一个参数为罐中普通球的数量，第二个参数为“倒霉球”的数量，第三个参数为每次采样（摸球）的数量。

points = np.zeros(100)
outcomes = np.random.hypergeometric(25, 1, 3, size=len(points))

(2) 根据上一步产生的游戏结果计算相应的得分。

for i in range(len(points)):
    if outcomes[i] == 3:
        points[i] = points[i - 1] + 1
    elif outcomes[i] == 2:
        points[i] = points[i - 1] - 6
    else:
        print outcomes[i]

(3) 使用Matplotlib绘制points数组。

plot(np.arange(len(points)), points)
show()

6.22 正态分布

连续分布可以用PDF（Probability Density Function，概率密度函数）来描述。随机变量落在某一区间内的概率等于概率密度函数在该区间的曲线下方的面积。
(1) 使用NumPy random模块中的normal函数产生指定数量的随机数。
N=10000

normal_values = np.random.normal(size=N)

(2) 绘制分布直方图和理论上的概率密度函数（均值为0、方差为1的正态分布）曲线。使用Matplotlib进行绘图。

dummy, bins, dummy = plt.hist(normal_values, np.sqrt(N), normed=True, lw=1)
sigma = 1
mu = 0
plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) * np.exp( - (bins -mu)**2 / (2 *sigma**2) ),lw=2)
plt.show()

6.24 对数正态分布

对数正态分布（lognormal distribution） 是自然对数服从正态分布的任意随机变量的概率分布。
(1) 使用NumPy random模块中的normal函数产生随机数。

N=10000
lognormal_values = np.random.lognormal(size=N)

(2) 绘制分布直方图和理论上的概率密度函数（均值为0、方差为1）。我们将使用Matplotlib进行绘图。

dummy, bins, dummy = plt.hist(lognormal_values,np.sqrt(N), normed=True, lw=1)
sigma = 1
mu = 0
x = np.linspace(min(bins), max(bins), len(bins))
pdf = np.exp(-(numpy.log(x) - mu)**2 / (2 * sigma**2))/ (x *sigma * np.sqrt(2 * np.pi))
plt.plot(x, pdf,lw=3)
plt.show()

7.1 排序

NumPy提供了多种排序函数，如下所示：
 sort函数返回排序后的数组；
 lexsort函数根据键值的字典序进行排序；
 argsort函数返回输入数组排序后的下标；
 ndarray类的sort方法可对数组进行原地排序；
 msort函数沿着第一个轴排序；
 sort_complex函数对复数按照先实部后虚部的顺序进行排序。

7.2 按字典序排序

lexsort函数返回输入数组按字典序排序后的下标
Sort names: first by surname, then by name.

>>> surnames =    ('Hertz',    'Galilei', 'Hertz')
>>> first_names = ('Heinrich', 'Galileo', 'Gustav')
>>> ind = np.lexsort((first_names, surnames))
>>> ind
array([1, 2, 0])
>>>
>>> [surnames[i] + ", " + first_names[i] for i in ind]
['Galilei, Galileo', 'Hertz, Gustav', 'Hertz, Heinrich']

Sort two columns of numbers:

>>> a = [1,5,1,4,3,4,4] # First column
>>> b = [9,4,0,4,0,2,1] # Second column
>>> ind = np.lexsort((b,a)) # Sort by a, then by b
>>> print ind
[2 0 4 6 5 3 1]
>>>
>>> [(a[i],b[i]) for i in ind]
[(1, 0), (1, 9), (3, 0), (4, 1), (4, 2), (4, 4), (5, 4)]

7.4 对复数进行排序（略）

7.5 搜索

 argmax函数返回数组中最大值对应的下标。

>>> a = np.array([2, 4, 8])
>>> np.argmax(a)
2

 nanargmax函数提供相同的功能，但忽略NaN值。

>>> b = np.array([np.nan, 2, 4])
>>> np.nanargmax(b)
2

 argmin和nanargmin函数的功能类似，只不过换成了最小值。
 argwhere函数根据条件搜索非零的元素，并分组返回对应的下标。

>>> a = np.array([2, 4, 8])
>>> np.argwhere(a <= 4)
array([[0],[1]])

 searchsorted函数可以为指定的插入值寻找维持数组排序的索引位置。该函数使用二分搜索算法，计算复杂度为O(log(n))。
 extract函数返回满足指定条件的数组元素。

7.6 使用searchsorted 函数

searchsorted函数为7和-2返回了索引5和0。用这些索引作为插入位置，生成数组[-2, 0, 1, 2, 3, 4, 7]，这样就维持了数组的排序

import numpy as np
a = np.arange(5)
indices = np.searchsorted(a, [-2, 7])  #索引可以维持数组排序的插入位置
print "Indices", indices
#output
Indices [0 5]

print "The full array", np.insert(a, indices, [-2, 7])
#output
The full array [-2 0 1 2 3 4 7]

7.8 从数组中抽取元素

从一个数组中抽取偶数元素。

import numpy as np
a = np.arange(7)
condition = (a % 2) == 0 #生成选择偶数元素的条件变量
print "Even numbers", np.extract(condition, a)
print "Non zero", np.nonzero(a)

使用nonzero函数抽取数组中的非零元素
np.nonzero(a)

7.9 金融函数

 fv函数计算所谓的终值（future value），即基于一些假设给出的某个金融资产在未来某一时间点的价值。
 pv函数计算现值（present value），即金融资产当前的价值。
 npv函数返回的是净现值（net present value），即按折现率计算的净现金流之和。
 pmt函数根据本金和利率计算每期需支付的金额。
 irr函数计算内部收益率（internal rate of return）。内部收益率是是净现值为0时的有效利率，不考虑通胀因素。
 mirr函数计算修正后内部收益率（modified internal rate of return），是内部收益率的改进版本。
 nper函数计算定期付款的期数。
 rate函数计算利率（rate of interest）。

7.10 计算终值

终值是基于一些假设给出的某个金融资产在未来某一时间点的价值。终值决定于4个参数——利率、期数、每期支付金额以及现值。在本节的教程中，我们以利率3%、每季度支付金额10、存款周期5年以及现值1 000为参数计算终值。

使用正确的参数调用fv函数，计算终值：

print "Future value", np.fv(0.03/4, 5 * 4, -10, -1000)
#output
Future value 1376.09633204

7.12 计算现值（略）

7.14 计算净现值（略）

7.16 计算内部收益率（略）

7.18 计算分期付款

pmt函数可以根据利率和期数计算贷款每期所需支付的资金。
假设贷款100万，年利率为10%，要用30年时间还完贷款，那么每月必须支付多少资金呢？使用刚才提到的参数值，调用pmt函数。

print "Payment", np.pmt(0.10/12, 12 * 30, 1000000)
#output
Payment -8775.71570089

7.20 计算付款期数

NumPy中的nper函数可以计算分期付款所需的期数。所需的参数为贷款利率、固定的月供以及贷款额。
考虑贷款9000，年利率10%，每月固定还款为100的情形。
通过nper函数计算出付款期数。

print "Number of payments", np.nper(0.10/12, -100, 9000)
#output
Number of payments 167.047511801 #需要167个月

7.22 计算利率

rate函数根据给定的付款期数、每期付款资金、现值和终值计算利率。
使用7.20节中的数值进行逆向计算，由其他参数得出利率。
填入之前教程中的数值作为参数。

print "Interest rate", 12 * np.rate(167, -100, 9000, 0)
#计算出的利率约为10%。
Interest rate 0.0999756420664

7.23 窗函数

窗函数（window function）是信号处理领域常用的数学函数，相关应用包括谱分析和滤波器设计等。这些窗函数除在给定区间之外取值均为0。NumPy中有很多窗函数，如bartlett、blackman、hamming、hanning和kaiser。

7.24 绘制巴特利特窗

巴特利特窗（Bartlett window）是一种三角形平滑窗。

window = np.bartlett(42)
plot(window)
show()

7.25 布莱克曼窗

布莱克曼窗（Blackman window）形式上为三项余弦值的加和。blackman函数返回布莱克曼窗。该函数唯一的参数为输出点的数量。如果数量为0或小于0，则返回一个空数组。

7.26 使用布莱克曼窗平滑股价数据（略）

7.27 汉明窗

汉明窗（Hamming window）形式上是一个加权的余弦函数。NumPy中的hamming函数返回汉明窗。该函数唯一的参数为输出点的数量。如果数量为0或小于0，则返回一个空数组。

7.28 绘制汉明窗

window = np.hamming(42)
plot(window)
show()

7.29 凯泽窗

凯泽窗（Kaiser window）是以贝塞尔函数（Bessel function）定义的，该函数的第一个数为输出点的数量。如果数量为0或小于0，则返回一个空数组。第二个参数为β值。

window = np.kaiser(42, 14)
plot(window)
show()

7.31 专用数学函数（略）

7.32 绘制修正的贝塞尔函数（略）

7.34 绘制sinc 函数（略）

8.1 断言函数

numpy.testing包：
- assert_almost_equal 如果两个数字的近似程度没有达到指定精度，就抛出异常
- assert_approx_equal 如果两个数字的近似程度没有达到指定有效数字，就抛出异常
- assert_array_almost_equal 如果两个数组中元素的近似程度没有达到指定精度，就抛出异常
- assert_array_equal 如果两个数组对象不相同，就抛出异常
- assert_array_less 两个数组必须形状一致，并且第一个数组的元素严格小于第二个数组的元素，否则就抛出异常
- assert_equal 如果两个对象不相同，就抛出异常
- assert_raises 若用填写的参数调用函数没有抛出指定的异常，则测试不通过
- assert_warns 若没有抛出指定的警告，则测试不通过
- assert_string_equal 断言两个字符串变量完全相同
- assert_allclose 如果两个对象的近似程度超出了指定的容差限，就抛出异常

8.2 使用assert_almost_equal 断言近似相等

(1) 调用函数，指定较低的精度（小数点后7位）：

print "Decimal 6", np.testing.assert_almost_equal(0.123456789, 0.123456780,decimal=7)
#output
Decimal 6 None

(2) 调用函数，指定较高的精度（小数点后8位）：

print "Decimal 7", np.testing.assert_almost_equal(0.123456789, 0.123456780,decimal=8)

#output
Decimal 7
Traceback (most recent call last):
...
raiseAssertionError(msg)
AssertionError:
Arrays are not almost equal
ACTUAL: 0.123456789
DESIRED: 0.12345678

8.4 使用assert_approx_equal 断言近似相等

与assert_almost_equal 类似，以小数点后7位为分界点。

8.6 断言数组近似相等

assert_array_almost_equal函数将抛出异常。该函数首先检查两个数组的形状是否一致，然后逐一比较两个数组中的元素。
(1) 调用函数，指定较低的精度：

print "Decimal 8", np.testing.assert_array_almost_equal([0, 0.123456789], [0,
0.123456780], decimal=8)
#output
Decimal 8 None

print "Decimal 9", np.testing.assert_array_almost_equal([0, 0.123456789], [0,0.123456780], decimal=9)

#output
Decimal 9
Traceback (most recent call last):
...
assert_array_compare
raiseAssertionError(msg)
AssertionError:
Arrays are not almost equal
(mismateh 50.0%)
x: array([ 0. , 0.12345679])
y: array([ 0. , 0.12345678])

• 
  

8.8 比较数组相等

(1) 调用assert_allclose函数：

print "Pass", np.testing.assert_allclose([0, 0.123456789, np.nan], [0, 0.123456780,np.nan], rtol=1e-7, atol=0)
#output
Pass None

(2) 调用assert_array_equal函数：

print "Fail", np.testing.assert_array_equal([0, 0.123456789, np.nan], [0, 0.123456780,np.nan])
#output
Fail
Traceback (most recent call last):
...
assert_array_compare
raiseAssertionError(msg)
AssertionError:
Arrays are not equal
(mismatch 50.0%)
x: array([ 0. ,0.12345679, nan]
y: array([ 0. ,0.12345678, nan])

8.10 核对数组排序

两个数组必须形状一致并且第一个数组的元素严格小于第二个数组的元素，否则assert_array_less函数将抛出异常。

print "Pass", np.testing.assert_array_less([0, 0.123456789, np.nan], [1, 0.23456780,np.nan])
#output
Pass None

8.12 比较对象

如果两个对象不相同，assert_equal函数将抛出异常。这里的对象不一定为NumPy数组，也可以是Python中的列表、元组或字典。

(1) 调用assert_equal函数：
print "Equal?", np.testing.assert_equal((1, 2), (1, 3))

#output
Equal?
Traceback (most recent call last):
...
raiseAssertionError(msg)
AssertionError:
Items are not equal:
item=1
ACTUAL: 2
DESIRED: 3

8.14 比较字符串

assert_string_equal函数断言两个字符串变量完全相同。如果测试不通过，将会抛出异常并显示两个字符串之间的差异。该函数区分字符大小写。

print "Fail", np.testing.assert_string_equal("NumPy", "Numpy")

#output
Fail
Traceback (most recent call last):
...
raiseAssertionError(msg)
AssertionError: Differences in strings:
- NumPy? ^
+ Numpy?

8.15 浮点数比较

浮点数在计算机中是以不精确的方式表示的，这给比较浮点数带来了问题。NumPy中的assert_array_almost_equal_nulp和assert_array_max_ulp函数可以提供可靠的浮点数比较功能。ULP是Unit of Least Precision的缩写，即浮点数的最小精度单位。根据IEEE 754标准，四则运算的误差必须保持在半个ULP之内。你可以用刻度尺来做对比。公制刻度尺的刻度通常精确到毫米，而更高精度的部分只能估读，误差上界通常认为是最小刻度值的一半，即半毫米。
机器精度（machine epsilon）是指浮点运算中的相对舍入误差上界。机器精度等于ULP相对于1的值。NumPy中的finfo函数可以获取机器精度。Python标准库也可以给出机器精度值，并应该与NumPy给出的结果一致。

8.16 使用assert_array_almost_equal_nulp 比较浮点数（略）

8.18 设置maxulp 并比较浮点数（略）

8.20 编写单元测试

import numpy as np
import unittest

def factorial(n):
    if n == 0:
        return 1
    if n < 0:
        raise ValueError, "Unexpected negative value"
    return np.arange(1, n+1).cumprod()

class FactorialTest(unittest.TestCase):
    def test_factorial(self):
    # 计算3的阶乘，测试通过
    self.assertEqual(6, factorial(3)[-1])
    np.testing.assert_equal(np.array([1, 2, 6]), factorial(3))

    def test_zero(self):
    # 计算0的阶乘，测试通过
    self.assertEqual(1, factorial(0))

    def test_negative(self):
    # 计算负数的阶乘，测试不通过
    # 这里应抛出ValueError异常，但我们断言其抛出IndexError异常
    self.assertRaises(IndexError, factorial(-10))

if _name_ == '_main_':
unittest.main()

8.21 nose 和测试装饰器

nose也是一种Python框架，使得（单元）测试更加容易。nose可以帮助你组织测试代码。根据nose的文档，任何能够匹配testMatch正则表达式（默认为(?:^|[b_.-])[Tt]est）的Python源代码文件、文件夹或库都将被收集用于测试。nose充分利用了装饰器（decorator）。Python装饰器是有一定含义的对函数或方法的注解。numpy.testing模块中有很多装饰器。

8.22 使用测试装饰器（略）

8.24 执行文档字符串测试

文档字符串（docstring）是内嵌在Python代码中的类似交互式会话的字符串。这些字符串可以用于某些测试，也可以仅用于提供使用示例。numpy.testing模块中有一个函数可以运行这些测试。

docstringtest.py

import numpy as np
def factorial(n):
    """
    Test for the factorial of 3 that should pass.
    >>> factorial(3)
    6
    Test for the factorial of 0 that should fail.
    >>> factorial(0)
    1
    """
return np.arange(1, n+1).cumprod()[-1]

>>>from numpy.testing import rundocs
>>>rundocs('docstringtest.py')

9.2 绘制多项式函数

import numpy as np
import matplotlib.pyplot as plt
func = np.poly1d(np.array([1, 2, 3, 4]).astype(float))
x = np.linspace(-10, 10, 30)
y = func(x)
plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.show()

9.4 绘制多项式函数及其导函数

import numpy as np
import matplotlib.pyplot as plt
func = np.poly1d(np.array([1, 2, 3, 4]).astype(float))
func1 = func.deriv(m=1)
x = np.linspace(-10, 10, 30)
y = func(x)
y1 = func1(x)
plt.plot(x, y, 'ro', x, y1, 'g--')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

9.6 绘制多项式函数及其导函数

import numpy as np
import matplotlib.pyplot as plt

func = np.poly1d(np.array([1, 2, 3, 4]).astype(float))
x = np.linspace(-10, 10, 30)
y = func(x)
func1 = func.deriv(m=1) #一次导数
y1 = func1(x)
func2 = func.deriv(m=2) #二次导数
y2 = func2(x)
plt.subplot(311)
plt.plot(x, y, 'r-' )
plt.title("Polynomial")
plt.subplot(312)
plt.plot(x, y1, 'b^')
plt.title("First Derivative")
plt.subplot(313)
plt.plot(x, y2, 'go')
plt.title("Second Derivative")
plt.xlabel('x')
plt.ylabel('y')
plt.show()

9.8 绘制全年股票价格

from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator
from matplotlib.finance import quotes_historical_yahoo
from matplotlib.finance import candlestick
import sys
from datetime import date
import matplotlib.pyplot as plt

today = date.today()
start = (today.year - 1, today.month, today.day) #将当前的日期减去1年作为起始日期

#创建定位器（locator），这些来自matplotlib.dates包中的对象可以在x轴上定位月份和日期。
alldays = DayLocator()
months = MonthLocator()

#创建一个日期格式化器（date formatter）以格式化x轴上的日期。该格式化器将创建一个字符串，包含简写的月份和年份。
month_formatter = DateFormatter("%b %Y")

#股票名称
symbol = 'DISH'
if len(sys.argv) == 2:
    symbol = sys.argv[1]

#从雅虎财经频道下载股价数据。
quotes = quotes_historical_yahoo(symbol, start, today)
fig = plt.figure()
ax = fig.add_subplot(111)

#将x轴上的主定位器设置为月定位器。该定位器负责x轴上较粗的刻度。
ax.xaxis.set_major_locator(months)

#将x轴上的次定位器设置为日定位器。该定位器负责x轴上较细的刻度。
ax.xaxis.set_minor_locator(alldays)

#将x轴上的主格式化器设置为月格式化器。该格式化器负责x轴上较粗刻度的标签。
ax.xaxis.set_major_formatter(month_formatter)

#可以指定K线图的矩形宽度，现在先使用默认值。
candlestick(ax, quotes)

#将x轴上的标签格式化为日期。为了更好地适应x轴的长度，标签将被旋转。
fig.autofmt_xdate()
plt.show()

9.10 绘制股价分布直方图

from matplotlib.finance import quotes_historical_yahoo
import sys
from datetime import date
import matplotlib.pyplot as plt
import numpy as np
today = date.today()
start = (today.year - 1, today.month, today.day)
symbol = 'DISH'
if len(sys.argv) == 2:
    symbol = sys.argv[1]
quotes = quotes_historical_yahoo(symbol, start, today)
quotes = np.array(quotes)
close = quotes.T[4] #提取出收盘价数据
plt.hist(close, np.sqrt(len(close)))
plt.show()

9.12 绘制股票成交量

当数据的变化范围很大时，对数坐标图（logarithmic plot）很有用。Matplotlib中有semilogx函数（对x轴取对数）、semilogy函数（对y轴取对数）和loglog函数（同时对x轴和y轴取对数）。

from matplotlib.finance import quotes_historical_yahoo
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator
import sys
from datetime import date
import matplotlib.pyplot as plt
import numpy as np
today = date.today()
start = (today.year - 1, today.month, today.day)

symbol = 'DISH'
if len(sys.argv) == 2:
    symbol = sys.argv[1]

quotes = quotes_historical_yahoo(symbol, start, today)
quotes = np.array(quotes)
dates = quotes.T[0]
volume = quotes.T[5]
alldays = DayLocator()
months = MonthLocator()
month_formatter = DateFormatter("%b %Y")

fig = plt.figure()
ax = fig.add_subplot(111)
plt.semilogy(dates, volume)
ax.xaxis.set_major_locator(months)
ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_formatter(month_formatter)
fig.autofmt_xdate()
plt.show()

9.14 绘制股票收益率和成交量变化的散点图

from matplotlib.finance import quotes_historical_yahoo
import sys
from datetime import date
import matplotlib.pyplot as plt
import numpy as np
today = date.today()
start = (today.year - 1, today.month, today.day)

symbol = 'DISH'
if len(sys.argv) == 2:
    symbol = sys.argv[1]

quotes = quotes_historical_yahoo(symbol, start, today)
quotes = np.array(quotes)
close = quotes.T[4]
volume = quotes.T[5]

#计算股票收益率和成交量的变化值
ret = np.diff(close)/close[:-1] #out[n] = a[n+1] - a[n] 
volchange = np.diff(volume)/volume[:-1]

fig = plt.figure()
ax = fig.add_subplot(111)
#数据点的颜色与股票收益率相关联，数据点的大小与成交量的变化相关联 c,color s,size
ax.scatter(ret, volchange, c=ret * 100, s=volchange * 100, alpha=0.5)
ax.set_title('Close and volume returns')
ax.grid(True)
plt.show()

9.16 根据条件进行着色

fill_between函数使用指定的颜色填充图像中的区域。我们也可以选择alpha通道的取值。该函数的where参数可以指定着色的条件。

对股价图进行了着色，低于平均值的收盘价使用了一种颜色，高于平均值的收盘价使用了另外一种不同的颜色。

from matplotlib.finance import quotes_historical_yahoo
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator
import sys
from datetime import date
import matplotlib.pyplot as plt
import numpy as np
today = date.today()
start = (today.year - 1, today.month, today.day)
symbol = 'DISH'
if len(sys.argv) == 2:
    symbol = sys.argv[1]
quotes = quotes_historical_yahoo(symbol, start, today)
quotes = np.array(quotes)
dates = quotes.T[0]
close = quotes.T[4]
alldays = DayLocator()
months = MonthLocator()
month_formatter = DateFormatter("%b %Y")
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(dates, close)

#对收盘价下方的区域进行着色，依据低于或高于平均收盘价使用不同的颜色填充。
plt.fill_between(dates, close.min(), close, where=close>close.mean(), facecolor="green",alpha=0.4)
plt.fill_between(dates, close.min(), close, where=close<close.mean(),facecolor="red", alpha=0.4)

ax.xaxis.set_major_locator(months)
ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_formatter(month_formatter)
ax.grid(True)
fig.autofmt_xdate()
plt.show()

9.18 使用图例和注释

以用legend函数创建透明的图例，并由Matplotlib自动确定其摆放位置。可以用annotate函数在图像上精确地添加注释，并有很多可选的注释和箭头风格。

from matplotlib.finance import quotes_historical_yahoo
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator
import sys
from datetime import date
import matplotlib.pyplot as plt
import numpy as np

today = date.today()
start = (today.year - 1, today.month, today.day)

symbol = 'DISH'
if len(sys.argv) == 2:
    symbol = sys.argv[1]

quotes = quotes_historical_yahoo(symbol, start, today)
quotes = np.array(quotes)
dates = quotes.T[0]
close = quotes.T[4]

fig = plt.figure()
ax = fig.add_subplot(111)

#计算并绘制指数移动平均线(第3章),分别使用9、12和15作为周期数计算和绘制指数移动平均线。添加了一个图例，并使用注释将其中两条曲线的交点标注了出来。
其中两条曲线的交点标注了出来
其中两条曲线的交点标注了出来。示
方法。分别使用9、12和15作为周期数计算和绘制指数移动平均线。
emas = []
for i in range(9, 18, 3):
    weights = np.exp(np.linspace(-1., 0., i))
    weights / = weights.sum()
    ema = np.convolve(weights, close)[i-1:-i+1]
    idx = (i - 6)/3
    ax.plot(dates[i-1:], ema, lw=idx, label="EMA(%s)" % (i))
    data = np.column_stack((dates[i-1:], ema))
    emas.append(np.rec.fromrecords(data, names=["dates", "ema"]))

#找到两条指数移动平均曲线的交点。
first = emas[0]["ema"].flatten()
second = emas[1]["ema"].flatten()
bools = np.abs(first[-len(second):] - second)/second < 0.0001
xpoints = np.compress(bools, emas[1])

#将找到的交点用注释和箭头标注出来，并确保注释文本在交点的不远处。
for xpoint in xpoints:
    ax.annotate('x', xy=xpoint, textcoords='offset points',xytext=(-50, 30),arrowprops=dict(arrowstyle="->"))

#添加一个图例并由Matplotlib自动确定其摆放位置。
leg = ax.legend(loc='best', fancybox=True)

#设置alpha通道值，将图例透明化。
leg.get_frame().set_alpha(0.5)

alldays = DayLocator()
months = MonthLocator()
month_formatter = DateFormatter("%b %Y")
ax.plot(dates, close, lw=1.0, label="Close")
ax.xaxis.set_major_locator(months)
ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_formatter(month_formatter)
ax.grid(True)
fig.autofmt_xdate()
plt.show()

9.20 在三维空间中绘图

绘制 z = x**2 + y**2

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm

fig = plt.figure()
#使用3d关键字来指定图像的三维投影。
ax = fig.add_subplot(111, projection='3d')

#使用meshgrid函数创建一个二维的坐标网格。
u = np.linspace(-1, 1, 100)
x, y = np.meshgrid(u, u)
z = x ** 2 + y ** 2 #
ax.plot_surface(x, y, z, rstride=4, cstride=4, cmap=cm.YlGnBu_r) #cmap colormap
plt.show()

9.22 绘制色彩填充的等高线图

Matplotlib中的等高线3D绘图有两种风格——填充的(contourf)和非填充的(contour)。

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm
fig = plt.figure()
ax = fig.add_subplot(111)
u = np.linspace(-1, 1, 100)
x, y = np.meshgrid(u, u)
z = x ** 2 + y ** 2
ax.contourf(x, y, z)
plt.show()

9.24 制作动画

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

fig = plt.figure()
ax = fig.add_subplot(111)
N = 10
x = np.random.rand(N)
y = np.random.rand(N)
z = np.random.rand(N)

#用不同颜色的圆形、小圆点和三角形来绘制三个数据集中的数据点。
circles, triangles, dots = ax.plot(x, 'ro', y, 'g^', z, 'b.')
ax.set_ylim(0, 1)
plt.axis('off')

#定期调用以更新屏幕上的内容。将随机更新两个数据集中的y坐标值。
def update(data):
    circles.set_ydata(data[0])
    triangles.set_ydata(data[1])
    return circles, triangles

#使用NumPy生成随机数。  
def generated:
    while True: yield np.random.rand(2, N)

anim = animation.FuncAnimation(fig, update, generate, interval=150)
plt.show()

10.2 保存和加载.mat 文件

MATLAB以及其开源替代品Octave都是流行的数学工具。scipy.io包的函数可以在Python中加载或保存MATLAB和Octave的矩阵和数组。loadmat函数可以加载.mat文件。savemat函数可以将数组和指定的变量名字典保存为.mat文件。

import numpy as np
from scipy import io
a = np.arange(7)
io.savemat("a.mat", {"array": a})

10.4 分析随机数

from scipy import stats
import matplotlib.pyplot as plt

generated = stats.norm.rvs(size=900)  #使用scipy.stats包按正态分布生成随机数。
print "Mean", "Std", stats.norm.fit(generated)  #用正态分布去拟合生成的数据，得到其均值和标准差。

#偏度（skewness）描述的是概率分布的偏斜（非对称）程度。偏度检验有两个返回值，其中第二个返回值为p-value，即观察到的数据集服从正态分布的概率，取值范围为0~1。
print "Skewtest", "pvalue", stats.skewtest(generated) 

#output
#Skewtest pvalue (-0.62120640688766893, 0.5344638245033837)  
#该数据集有53%的概率服从正态分布。

#峰度（kurtosis）描述的是概率分布曲线的陡峭程度。
print "Kurtosistest","pvalue",stats.kurtosistest(generated)

#正态性检验（normality test）可以检查数据集服从正态分布的程度。我们来做一个正态性检验。该检验同样有两个返回值，其中第二个返回值为p-value。
print "Normaltest", "pvalue", stats.normaltest(generated)

#得到95%处的数值如下
print "95 percentile", stats.scoreatpercentile(generated, 95)

#将前一步反过来，可以从数值1出发找到对应的百分比
print "Percentile at 1", stats.percentileofscore(generated, 1)


plt.hist(generated)
plt.show()

10.6 比较股票对数收益率（略）

Jarque-Bera检验基于数据样本的偏度和峰度，评价给定数据服从未知均值和方差正态分布的假设是否成立。

10.8 检测QQQ 股价的线性趋势

以detrend函数作为滤波器，该函数可以对信号进行线性拟合，然后从原始输入数据中去除这个线性趋势。

from matplotlib.finance import quotes_historical_yahoo
from datetime import date
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
from matplotlib.dates import DayLocator
from matplotlib.dates import MonthLocator

#获取QQQ的收盘价和对应的日期数据
today = date.today()
start = (today.year - 1, today.month, today.day)
quotes = quotes_historical_yahoo("QQQ", start, today)
quotes = np.array(quotes)
dates = quotes.T[0]
qqq = quotes.T[4]

#去除信号中的线性趋势
y = signal.detrend(qqq)

#创建月定位器和日定位器
alldays = DayLocator()
months = MonthLocator()

#创建一个日期格式化器以格式化x轴上的日期。该格式化器将创建一个字符串，包含简写的月份和年份
month_formatter = DateFormatter("%b %Y")

fig = plt.figure()
ax = fig.add_subplot(111)

#绘制股价数据以及将去除趋势后的信号从原始数据中减去所得到的潜在趋势
plt.plot(dates, qqq, 'o', dates, qqq - y, '-')

#设置定位器和格式化器
ax.xaxis.set_minor_locator(alldays)
ax.xaxis.set_major_locator(months)
ax.xaxis.set_major_formatter(month_formatter)

#将x轴上的标签格式化为日期
fig.autofmt_xdate()
plt.show()

10.10 对去除趋势后的信号进行滤波处理（略）

现实世界中的信号往往具有周期性。傅里叶变换（Fourier transform）是处理这些信号的常用工具。傅里叶变换是一种从时域到频域的变换，也就是将周期信号线性分解为不同频率的正弦和余弦函数。
傅里叶变换的函数可以在scipy.fftpack模块中找到（NumPy也有自己的傅里叶工具包即numpy.fft）。这个模块包含快速傅里叶变换、微分算子和拟微分算子以及一些辅助函数。

10.12 拟合正弦波（略）

在scipy.optimize模块中提供了一些优化算法，最小二乘法函数leastsq就是其中之一。当调用这个函数时，我们需要提供一个残差（误差项）函数。这样，leastsq将最小化残差的平方和。
得到的解与我们使用的数学模型有关。我们还需要为算法提供一个起始点，这应该是一个最好的猜测——尽可能接近真实解。否则，程序执行800轮迭代后将停止。

10.14 计算高斯积分（略）

quad函数可以求单变量函数在两点之间的积分，这些点之间的距离可以是无穷小或无穷大。该函数使用最简单的数值积分方法即梯形法则（trapezoid rule）进行计算。

高斯积分（Gaussian integral）出现在误差函数（数学中记为erf）的定义中，但高斯积分本身的积分区间是无穷的，它的值等于pi的平方根。我们将使用quad函数计算它。

10.16 一维插值

插值（interpolation）即在数据集已知数据点之间“填补空白”。scipy.interpolate函数可以根据实验数据进行插值。interp1d类可以创建线性插值（linear interpolation）或三次插值（cubic interpolation）的函数。默认将创建线性插值函数，三次插值函数可以通过设置kind参数来创建。interp2d类的工作方式相同，只不过用于二维插值。

import numpy as np
from seipy import interpolate
import matplotlib.pyplot as plt

#创建数据点并添加噪音
x = np.linspaee(-18, 18, 36)
noise = 0.1 * np.random.random(len(x))
signal = np.sinc(x) + noise

#创建一个线性插值函数，并应用于有5倍数据点个数的输入数组
interpreted = interpolate.interpld(x, signal)
x2 = np.linspace(-18, 18, 180)
y = interpreted(x2)

#执行与前一步相同的操作，不过这里使用三次插值
cubic = interpolate.interpld(x, signal, kind="cubic")
y2 = cubic(x2)

plt.plot(x, signal, 'o', label="data")
plt.plot(x2, y, '-', label="linear")
plt.plot(x2, y2, '-', lw=2, label="cubic" )
plt.legend()
plt.show()

10.18 处理Lena 图像

使用scipy.ndimage包进行图像处理。该模块包含各种图像滤波器和工具函数。

from scipy import misc
import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage

#载入Lena图像，并使用灰度颜色表将其在子图中显示出来
image = misc.lena().astype(np.float32)
plt.subplot(221)
plt.title("Original Image")
img = plt.imshow(image, cmap=plt.cm.gray)
plt.axis("off")

#中值滤波器扫描信号的每一个数据点，并替换为相邻数据点的中值。对图像应用中值滤波器并显示在第二个子图中。

plt.subplot(222)
plt.title("Median Filter")
filtered = ndimage.median_filter(image, size=(42,42))
plt.imshow(filtered, cmap=plt.cm.gray)
plt.axis("off" )

#旋转图像并显示在第三个子图中
plt.subplot(223)
plt.title("Rotated")
rotated = ndimage.rotate(image, 90)
plt.imshow(rotated, cmap=plt.cm.gray)
plt.axis("off")

#Prewitt滤波器是基于图像强度的梯度计算。对图像应用Prewitt滤波器并显示在第四个子图中
plt.subplot(224)
plt.title("Prewitt Filter")
filtered = ndimage.prewitt(image)

plt.imshow(filtered, cmap=plt.cm.gray)
plt.axis("off")
plt.show()

10.20 重复音频片段（略）

11.2 动手安装pygame(略)

11.4 制作简单游戏（略）

11.6 使用NumPy 和Pygame 制作动画对象（略）

11.8 在Pygame 中使用Matplotlib（略）

11.10 访问屏幕像素（略）

11.12 数据点聚类（略）

11.14 绘制谢尔宾斯基地毯（略）

11.16 模拟生命（略）

全书完