深入 NumPy 模块
# 来源：NumPy Biginner's Guide 2e ch6
矩阵的逆
import numpy as np
A = np.mat("0 1 2;1 0 3;4 -3 8")
print "A\n", A
'''
A
[[ 0 1 2]
[ 1 0 3]
[ 4 -3 8]]
'''
# 求解矩阵的逆，不可逆会报错
inverse = np.linalg.inv(A)
print "inverse of A\n", inverse
'''
inverse of A
[[-4.5 7. -1.5]
[-2. 4. -1. ]
[ 1.5 -2. 0.5]]
'''
print "Check\n", A * inverse
'''
Check
[[ 1. 0. 0.]
[ 0. 1. 0.]
[ 0. 0. 1.]]
'''
求解线性方程组
import numpy as np
A = np.mat("1 -2 1;0 2 -8;-4 5 9")
print "A\n", A
'''
A
[[ 1 -2 1]
[ 0 2 -8]
[-4 5 9]]
'''
b = np.array([0, 8, -9])
print "b\n", b
'''
b
[ 0 8 -9]
'''
# solve 求出向量 x，使 Ax = b
# 内部使用 np.dot(A.I, b) 来计算
# 所以 A 不可逆时报错
x = np.linalg.solve(A, b)
print "Solution", x
'''
Solution [ 29. 16. 3.]
'''
print "Check\n", np.dot(A , x)
'''
Check
[[ 0. 8. -9.]]
'''
特征值和特征向量
# 如果 Ax = λx，并且 x 非零
# 则 λ 是 A 的特征值，x 是 A 关于 λ 的特征向量
import numpy as np
A = np.mat("3 -2;1 0")
print "A\n", A
'''
A
[[ 3 -2]
[ 1 0]]
'''
# eigvals 仅仅求解特征值
print "Eigenvalues", np.linalg.eigvals(A)
# Eigenvalues [ 2. 1.]
# eig 求解特征值和特征向量
# 特征向量是数组的列
eigenvalues, eigenvectors = np.linalg.eig(A)
print "First tuple of eig", eigenvalues
print "Second tuple of eig\n", eigenvectors
'''
First tuple of eig [ 2. 1.]
Second tuple of eig
[[ 0.89442719 0.70710678]
[ 0.4472136 0.70710678]]
'''
for i in range(len(eigenvalues)):
print "Left", np.dot(A, eigenvectors[:,i])
print "Right", eigenvalues[i] * eigenvectors[:,i]
print
'''
Left [[ 1.78885438]
[ 0.89442719]]
Right [[ 1.78885438]
[ 0.89442719]]
Left [[ 0.70710678]
[ 0.70710678]]
Right [[ 0.70710678]
[ 0.70710678]]
'''
奇异值分解
# 奇异值是 A * A.T 特征值的算术平方根
# 若 A 是 mxn 阶矩阵
# 奇异值分解将 A 分解成 USV
# U 是 mxm 阶。S 是 mxm 阶对角矩阵，由奇异值构成
# V 是 mxn 阶
import numpy as np
A = np.mat("4 11 14;8 7 -2")
print "A\n", A
'''
A
[[ 4 11 14]
[ 8 7 -2]]
'''
# linalg.svd 求解奇异值
U, Sigma, V = np.linalg.svd(A, full_matrices=False)
print "U"
print U
'''
U
[[-0.9486833 -0.31622777]
[-0.31622777 0.9486833 ]]
'''
print "Sigma"
print Sigma
'''
Sigma
[ 18.97366596 9.48683298]
'''
print "V"
print V
'''
V
[[-0.33333333 -0.66666667 -0.66666667]
[ 0.66666667 0.33333333 -0.66666667]]
'''
print "Product\n", U * np.diag(Sigma) * V
'''
Product
[[ 4. 11. 14.]
[ 8. 7. -2.]]
'''
print np.linalg.eigvals(A * A.T) ** 0.5
'''
[ 18.97366596 9.48683298]
'''
广义逆矩阵
# 若存在唯一矩阵 M，满足
# AMA = A，MAM = M，AM 和 MA 均为对称矩阵
# 则 M 是 A 的广义逆矩阵
import numpy as np
A = np.mat("4 11 14;8 7 -2")
print "A\n", A
'''
A
[[ 4 11 14]
[ 8 7 -2]]
'''
# linalg.pinv 求解广义逆
pseudoinv = np.linalg.pinv(A)
print "Pseudo inverse\n", pseudoinv
'''
Pseudo inverse
[[-0.00555556 0.07222222]
[ 0.02222222 0.04444444]
[ 0.05555556 -0.05555556]]
'''
print "Check", A * pseudoinv
'''
Check [[ 1.00000000e+00 0.00000000e+00]
[ 8.32667268e-17 1.00000000e+00]]
'''
计算矩阵行列式
import numpy as np
A = np.mat("3 4;5 6")
print "A\n", A
'''
A
[[ 3. 4.]
[ 5. 6.]]
'''
# linalg.det 求解行列式
# 对于二阶矩阵
# det(A) = A[0][0] * A[1][1] - A[0][1] * A[1][0]
print "Determinant", np.linalg.det(A)
# Determinant -2.0
快速傅里叶变换（FFT）
import numpy as np
from matplotlib.pyplot import plot, show
# 创建 30 个点的余弦波
x = np.linspace(0, 2 * np.pi, 30)
wave = np.cos(x)
# 使用 fft 函数变换余弦波
transformed = np.fft.fft(wave)
# 验证是否能够还原数据
print np.all(np.abs(np.fft.ifft(transformed) - wave) < 10 ** -9)
# True
# 绘制变换后的信号
plot(transformed)
show()
调整频率
import numpy as np
from matplotlib.pyplot import plot, show
# 创建 30 个点的余弦波
x = np.linspace(0, 2 * np.pi, 30)
wave = np.cos(x)
# 使用 fft 函数变换余弦波
transformed = np.fft.fft(wave)
# 使用 fftshift 函数调整信号
shifted = np.fft.fftshift(transformed)
# 检查是否可以还原
print np.all(np.abs(np.fft.ifftshift(shifted) - transformed) < 10 ** -9)
# 绘制信号
plot(transformed, lw=2)
plot(shifted, lw=3)
show()
随机游走
import numpy as np
from matplotlib.pyplot import plot, show
# 创建大小为 10000 的余额数组
# 初始余额为 1000
cash = np.zeros(10000)
cash[0] = 1000
# 生成随机数，其概率满足二项分布 B(n = 9, p = 0.5)
# P(x = k) = C(n, k) * p ** k * (1 - p) ** (n - k)
# x: 0 ~ 9
outcome = np.random.binomial(9, 0.5, size=len(cash))
for i in range(1, len(cash)):
# 如果随机数小于 5，余额减一，否则加一
# p = 0.5 时，二项分布是对称的
# 所以小于 5 的概率应该是 0.5
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()
# 0 9
# 绘制余额变化情况
plot(np.arange(len(cash)), cash)
show()
模拟游戏节目
# 游戏节目中有一个桶，里面有 25 个正常球
# 和一个倒霉球
# 选手正确回答问题后都需要取出三个球
# 如果三个都是普通球，分数加一
# 否则分数减六
import numpy as np
from matplotlib.pyplot import plot, show
# 创建分数数组
points = np.zeros(100)
# hypergeometric 函数生成的随机数满足超几何分布
# 超几何分布描述了，袋子里有两类球，从袋子里取出几个球，取到 k 个一类球的概率。
# 参数分别为：普通球的个数，倒霉球的个数，取出的个数
# 随机数模拟了取出的球中普通球的个数
outcomes = np.random.hypergeometric(25, 1, 3, size=len(points))
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]
# 绘制分数变化情况
plot(np.arange(len(points)), points)
show()
绘制正态分布
import numpy as np
import matplotlib.pyplot as plt
N=10000
# 生成随机数，其概率满足正态分布 N(mu = 0, sigma = 1)
normal_values = np.random.normal(size=N)
# 绘制直方图，横轴是分组，纵轴是出现频率
# 第二个参数表示 100 个组（默认 10 个）
# bins 是组中值
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()
绘制对数正态分布
# 和上一段代码基本一样
import numpy as np
import matplotlib.pyplot as plt
N=10000
# 将这里改成 lognormal 来生成对数正态分布的随机数
lognormal_values = np.random.lognormal(size=N)
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(-(np.log(x) - mu)**2 / (2 * sigma**2))/ (x * sigma * np.sqrt(2 * np.pi))
plt.plot(x, pdf,lw=3)
plt.show()
