《Python数据分析基础教程：Numpy学习指南》- 速记 - 第十章

10.2 保存和加载.mat 文件

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

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

10.4 分析随机数

from scipy import statsimport matplotlib.pyplot as pltgenerated = 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_yahoofrom datetime import dateimport numpy as npfrom scipy import signalimport matplotlib.pyplot as pltfrom matplotlib.dates import DateFormatterfrom matplotlib.dates import DayLocatorfrom 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 npfrom seipy import interpolateimport 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 miscimport numpy as npimport matplotlib.pyplot as pltfrom 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 重复音频片段（略）