Numpy（四）模块

1、linalg模块

     线性代数是数学的一个重要分支。numpy.linalg模块包含线性代数的函数。使用这个模块，我们可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等。

1.1计算逆矩阵

import numpy as npa=np.mat('1 0;0 2')print a#逆矩阵print a.Iprint np.linalg.inv(a)#原矩阵*逆矩阵=单位矩阵print a.I*aprint np.linalg.inv(a)*a

1.2求解线性方程组

numpy.linalg中的函数solve可以求解形如 Ax = b 的线性方程组，其中 A 为矩阵，b 为一维或二维的数组，x 是未知变量

import numpy as np
# (1) 创建矩阵A和数组b
A=np.mat('1 2 -1;0 2 -8;-4 5 9')
b=np.array([0,8,-9])
print 'A=',A
print 'b=',b
# (2) 调用solve函数求解线性方程组：
x=np.linalg.solve(A,b)
print 'x=',x
#(3) 使用dot函数检查求得的解是否正确
print "check=",np.dot(A,x)

1.3特征值和特征向量

     特征值（eigenvalue）即方程 Ax = ax 的根，是一个标量。其中，A 是一个二维矩阵，x 是一个一维向量。特征向量（eigenvector）是关于特征值的向量。在numpy.linalg模块中，eigvals函数可以计算矩阵的特征值，而eig函数可以返回一个包含特征值和对应的特征向量的元组。

import numpy as np
# (1) 创建一个矩阵
A=np.mat('3 -2;1 0')
print 'A=',A
# (2) 调用eigvals函数求解特征值
print 'Eigenvalue',np.linalg.eigvals(A)
# (3) 使用eig函数求解特征值和特征向量。该函数将返回一个元组，按列排放着特征值和对应的特征向量，其中第一列为特征值，第二列为特征向量。
eigenvalues,eigenvectors=np.linalg.eig(A)
print "First tuple of eig", eigenvalues
print "Second tuple of eig\n", eigenvectors
# (4) 使用dot函数验证求得的解是否正确。分别计算等式 Ax = ax 的左半部分和右半部分，检查是否相等
for i in range(len(eigenvalues)):
    print "Left", np.dot(A, eigenvectors[:,i])
    print "Right", eigenvalues[i] * eigenvectors[:,i]
    print

1.4 奇异值分解
        SVD（Singular Value Decomposition，奇异值分解）是一种因子分解运算，将一个矩阵分解为3个矩阵的乘积。奇异值分解是前面讨论过的特征值分解的一种推广。在numpy.linalg模块中的svd函数可以对矩阵进行奇异值分解。该函数返回3个矩阵——U、Sigma和V，其中U和V是正交矩阵，Sigma包含输入矩阵的奇异值。

       星号表示厄米共轭（Hermitian conjugate）或共轭转置（conjugate transpose）。

import numpy as np
# (1) 首先，创建一个矩阵
A=np.mat('4 11 14;8 7 -2')
print 'A=',A
# (2) 使用svd函数分解矩阵
U,Sigma,V=np.linalg.svd(A,full_matrices=0)
print "U=",U
print "Sigma=",Sigma
print "V=",V
# (3) 不过，我们并没有真正得到中间的奇异值矩阵——得到的只是其对角线上的值，而非对角线上的值均为0。我们可以使用diag函数生成完整的奇异值矩阵。将分解出的3个矩阵相乘
print 'Produce\n',U*np.diag(Sigma)*V

1.5 广义逆矩阵
    使用numpy.linalg模块中的pinv函数进行求解,计算广义逆矩阵需要用到奇异值分解。inv函数只接受方阵作为输入矩阵，而pinv函数则没有这个限制。

import numpy as np
# (1) 首先，创建一个矩阵
A=np.mat('4 11 14;8 7 -2')
print 'A=',A
# (2) 使用pinv函数计算广义逆矩阵
preudoinv=np.linalg.pinv(A)
print 'pseudo inverse',preudoinv

# (3) 将原矩阵和得到的广义逆矩阵相乘
print 'Check=',A*preudoinv

1.6 行列式

       行列式（determinant）是与方阵相关的一个标量值，在数学中得到广泛应用。对于一个n×n的实数矩阵，行列式描述的是一个线性变换对“有向体积”所造成的影响。行列式的值为正表示保持了空间的定向（顺时针或逆时针），为负则表示颠倒了空间的定向。numpy.linalg模块中的det函数可以计算矩阵的行列式。

import numpy as np
# (1) 首先，创建一个矩阵
A=np.mat('3 4;5 6')
print "A=",A
# (2) 使用det函数计算行列式
print 'Determinant',np.linalg.det(A)

2、fft模块

  2.1 快速傅里叶变换

      FFT（Fast Fourier Transform，快速傅里叶变换）是一种高效的计算DFT（Discrete FourierTransform，离散傅里叶变换）的算法。FFT算法比根据定义直接计算更快，计算复杂度为O(NlogN) 。DFT在信号处理、图像处理、求解偏微分方程等方面都有应用。在NumPy中，有一个名为fft的模块提供了快速傅里叶变换的功能。在这个模块中，许多函数都是成对存在的，也就是说许多函数存在对应的逆操作函数。例如，fft和ifft函数就是其中的一对。

import numpy as np
from matplotlib.pyplot import plot, show
# (1) 创建一个包含30个点的余弦波信号
x=np.linspace(0,2*np.pi,30)
wave=np.cos(x)
# (2) 使用fft函数对余弦波信号进行傅里叶变换。
transformed=np.fft.fft(wave)
# (3) 对变换后的结果应用ifft函数，应该可以近似地还原初始信号。
print np.all(np.abs(np.fft.ifft(transformed)-wave)<10**-9)
# (4) 使用Matplotlib绘制变换后的信号。
plot(transformed)
show()

2.2移频
      numpy.linalg模块中的fftshift函数可以将FFT输出中的直流分量移动到频谱的中央。 ifftshift函数则是其逆操作。

import numpy as np
from matplotlib.pyplot import plot, show
# (1) 创建一个包含30个点的余弦波信号
x=np.linspace(0,2*np.pi,30)
wave=np.cos(x)
# (2) 使用fft函数对余弦波信号进行傅里叶变换。
transformed=np.fft.fft(wave)
# (3) 使用fftshift函数进行移频操作。
shifted=np.fft.fftshift(transformed)
# (4) 用ifftshift函数进行逆操作，这将还原移频操作前的信号。
print np.all((np.fft.ifftshift(shifted) - transformed) < 10 ** -9)
# print np.all((np.fft.ifftshift(shifted) - transformed) < 10 ** -9)
plot(transformed,lw=2)
plot(shifted,lw=3)
show()

3、随机数模块

3.1 随机数

      随机数在蒙特卡罗方法（Monto Carlo method）、随机积分等很多方面都有应用。真随机数的产生很困难，因此在实际应用中我们通常使用伪随机数。在大部分应用场景下，伪随机数已经足够随机，当然一些特殊应用除外。有关随机数的函数可以在NumPy的random模块中找到。随机数发生器的核心算法是基于马特赛特旋转演算法（Mersenne Twister algorithm）的。随机数可以从离散分布或连续分布中产生。分布函数有一个可选的参数size，用于指定需要产生的随机数的数量。该参数允许设置为一个整数或元组，生成的随机数将填满指定形状的数组。支持的离散分布包括几何分布、超几何分布和二项分布等。

import numpy as np
from matplotlib.pyplot import plot, show
# (1) 初始化一个全0的数组来存放剩余资本。以参数10000调用binomial函数，
cash=np.zeros(10000)
cash[0]=10000
outcome=np.random.binomial(9,0.5,size=len(cash))
# (2) 模拟每一轮抛硬币的结果并更新cash数组。打印出outcome的最小值和最大值，以检查输出中是否有任何异常值：
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')
print outcome.min(),outcome.max()
# (3) 使用Matplotlib绘制cash数组：
plot(np.arange(len(cash)),cash)
show()

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

import numpy as np
from matplotlib.pyplot import plot, show
# (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)
# (3)使用Matplotlib绘制points数组。
plot(np.arange(len(points)),points)
show()

4、连续分布和离散分布

4.1 连续分布

       连续分布可以用PDF（Probability Density Function，概率密度函数）来描述。随机变量落在某一区间内的概率等于概率密度函数在该区间的曲线下方的面积。NumPy的random模块中有一系列连续分布的函数——beta、chisquare、exponential、f、gamma、gumbel、laplace、lognormal、logistic、multivariate_normal、noncentral_chisquare、noncentral_f、normal等。

import numpy as np
import matplotlib.pyplot as plt
import math
# (1) 使用NumPy random模块中的normal函数产生指定数量的随机数。
N=10000
normal_values = np.random.normal(size=N)

# (2) 绘制分布直方图和理论上的概率密度函数（均值为0、方差为1的正态分布）曲线。我们将使用Matplotlib进行绘图。
n, bins, patches = plt.hist(normal_values, int(np.sqrt(N)), normed=1, facecolor='blue', alpha=0.5)
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()

4.2 对数正态分布

      对数正态分布（lognormal distribution） 是自然对数服从正态分布的任意随机变量的概率分布。NumPy random模块中的lognormal函数模拟了这个分布。

import numpy as np
import matplotlib.pyplot as plt
# (1) 使用NumPy random模块中的normal函数产生随机数。
N=10000
lognormal_values=np.random.lognormal(size=N)
# (2) 绘制分布直方图和理论上的概率密度函数（均值为0、方差为1）。我们将使用Matplotlib进行绘图。
dummy, bins, dummys = plt.hist(lognormal_values, int(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()