
SciPy是一个非常著名的开源计算库, 它专门用于科学研究领域, 这套库是依托在NumPy基础之上构建出来的, 它还额外提供了许多功能模块, 其中包括对数据进行数值积分的操作能力、实现最优化求解的功能、进行统计分析的工具以及调用各类专用函数的手段。1、先把那些.mat格式的文件给保存下来, 然后再去把它们加载进来。并且, 它们的开源替代方案同样是当下非常流行的数学计算工具, 其中 scipy.io 这个软件包里的各种函数, 能够被用来在中端进行数据矩阵还有数组的加载操作或是保存操作, 具体来说, 那个特定的函数可以负责把包含数据的 .mat 格式的文件给加载进来, 同时还有一个专门的函数, 它可以把数组数据连同指定的变量名字典一起组合好, 然后把这些内容统一保存成为标准的 .mat 文件格式。a np.arange(7)io.savemat(a.mat, {array: a})2、统计使用存储在scipy.stats包里面的那些统计方面的函数, 去对已经生成好的数据进行必要的分析操作工作。(1) 借助于scipy.stats模块来实施随机数的生成操作, 这些随机数需要满足正态分布的要求。 stats.norm.rvs(size900)(2) 使用当前的通用口语习惯, 强行转换成书面语的表达方式是这样我们拿那个叫做正态分布的东西去套那个刚刚做出来的数据, 然后就能得出一个平均数跟一个标准差。print Mean, Std, stats.norm.fit()(3) 偏度所描述的内容, 是关于概率分布的偏斜情况, 或者说它描述的是概率分布的非对称程度。针对这种情况, 接下来我们将实施一个统计检验动作, 这个步骤被称为进行偏度检验。这个检验操作会产生两个返回结果, 这里需要特别关注的是其中的第二个返回值, 也就是所谓的p-value。简单来说, p-value代表的是这样一个概率值, 即根据当前观察到的数据集来判定其能够完全符合正态分布的可能性, 现在取值。其取值区间是在零到一的范围之内的。先输出一个空字符串, 再次输出一个空字符串, 然后再去调用那个统计结果的函数。(4) 这个峰度是用来描述概率分布曲线的陡峭程度的这样一个概念, 然后, 我们现在是要来进行一个叫做峰度检验的这个动作, 那么这个检验它是跟。偏度检验的情况是类似的, 不过现在这里讨论的对象是针对峰度展开的分析。print , , stats.()(5) 正态性检验这一项操作也就是通常所说的test, 这个操作是用来对数据集是不是服从正态分布的这个状况进行检查的, 接下来我们要来进行一个关于正态性的检验工作。对数据进行检验。这个检验的操作结果里面, 包含有两个返回值的信息, 并且其中的第二个返回值的含义, 就是p-value这一项指标的值。print , , stats.()(6) 通过调用 SciPy 这个软件包, 我们能够很方便地拿到数据所在的区段中, 处于某一个百分比所在位置的具体数值, 此处示例为输出显示包含“95 ”这一内容以及执行相关统计函数并传入参数 95 的结果。(7) 我们还可以反过来, 从前面提到的那个步骤倒着回去操作一下, 这样我们也可以从数值1这个位置开始进行推导, 进而找到它所对应的那个百分比是多少。print at 1, stats.(, 1)我们按正态分布生成了一个随机数据集并使用scipy.stats模块分析了该数据集from scipy import statsimport matplotlib.pyplot as pltgenerated stats.norm.rvs(size900)print Mean, Std, stats.norm.fit(generated)print Skewtest, pvalue, stats.skewtest(generated)print Kurtosistest, pvalue, stats.kurtosistest(generated)print Normaltest, pvalue, stats.normaltest(generated)print 95 percentile, stats.scoreatpercentile(generated, 95)print Percentile at 1, stats.percentileofscore(generated, 1)plt.hist(generated)plt.show()3、信号处理scipy.模块中, 是包含有滤波函数的, 还有B样条插值B-函数。4、趋势分析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 MonthLocatortoday 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()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)fig.autofmt_xdate()plt.show()5、傅里叶分析傅里叶变换的函数可以在scipy这个模块中找到, 要知道NumPy也有自己的傅里叶工具包, 也就是numpy.fft。这个模块里面包含了快速傅里叶变换、微分算子以及拟微分算子, 同时还囊括了一些辅助函数。用户应该会感到高兴, 因为scipy模块中的许多函数与它们对应的函数名字是一样的, 而且这些函数的功能也差不多相近。去除了一个信号的趋势部分, 并且使用了scipy模块对其应用了一个简易的滤波器。from matplotlib.finance import quotes_historical_yahoofrom datetime import dateimport numpy as npfrom scipy import signalimport matplotlib.pyplot as pltfrom scipy import fftpackfrom matplotlib.dates import DateFormatterfrom matplotlib.dates import DayLocatorfrom matplotlib.dates import MonthLocatortoday 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()month_formatter DateFormatter(%b %Y)fig plt.figure()fig.subplots_adjust(hspace.3)ax fig.add_subplot(211)ax.xaxis.set_minor_locator(alldays)ax.xaxis.set_major_locator(months)ax.xaxis.set_major_formatter(month_formatter)# 调大字号ax.tick_params(axisboth, whichmajor, labelsizex-large)amps np.abs(fftpack.fftshift(fftpack.rfft(y)))amps[amps 0.1 * amps.max()] 0plt.plot(dates, y, o, labeldetrended)plt.plot(dates, -fftpack.irfft(fftpack.ifftshift(amps)), labelfiltered)fig.autofmt_xdate()plt.legend(prop{size:x-large})ax2 fig.add_subplot(212)ax2.tick_params(axisboth, whichmajor, labelsizex-large)N len(qqq)plt.plot(np.linspace(-N/2, N/2, N), amps, labeltransformed)plt.legend(prop{size:x-large})plt.show()6、数学优化优化算法会尝试去寻求某一个问题的最优解, 比如说找到函数的最大值或者是找到函数的最小值, 这个函数可以是线性函数, 也可以是非线性函数, 解可能会有一些特定的约束条件, 比如说不允许出现负数, 在scipy.模块里面提供了一些优化算法, 最小二乘法函数就是其中的一种, 当我们在调用这个函数的时候, 我们需要给出一个残差函数, 也就是误差项函数, 这样的话呢, 就能够将会把残差的平方和给最小化。最终的解是和我们所使用的数学模型直接关联着的, 除此之外, 我们还得为这个算法去设置一个初始的起始点。这个起始点应该被当成是最佳的猜测, 也就是要尽可能地去接近那个真实的解, 如果不这样做的话, 就极有可能导致在执行了大概八百轮迭代之后程序直接就会停止运行。使用scipy模块对滤波后的信号拟合了一个正弦波函数。from matplotlib.finance import quotes_historical_yahooimport numpy as npimport matplotlib.pyplot as pltfrom scipy import fftpackfrom scipy import signalfrom matplotlib.dates import DateFormatterfrom matplotlib.dates import DayLocatorfrom matplotlib.dates import MonthLocatorfrom scipy import optimizestart (2010, 7, 25)end (2011, 7, 25)quotes quotes_historical_yahoo(QQQ, start, end)quotes np.array(quotes)dates quotes.T[0]qqq quotes.T[4]y signal.detrend(qqq)alldays DayLocator()months MonthLocator()month_formatter DateFormatter(%b %Y)fig plt.figure()fig.subplots_adjust(hspace.3)ax fig.add_subplot(211)ax.xaxis.set_minor_locator(alldays)ax.xaxis.set_major_locator(months)ax.xaxis.set_major_formatter(month_formatter)ax.tick_params(axisboth, whichmajor, labelsizex-large)amps np.abs(fftpack.fftshift(fftpack.rfft(y)))amps[amps amps.max()] 0def residuals(p, y, x):A,k,theta,b perr y-A * np.sin(2* np.pi* k * x theta) breturn errfiltered -fftpack.irfft(fftpack.ifftshift(amps))N len(qqq)f np.linspace(-N/2, N/2, N)p0 [filtered.max(), f[amps.argmax()]/(2*N), 0, 0]print P0, p0plsq optimize.leastsq(residuals, p0, args(filtered, dates))p plsq[0]print P, pplt.plot(dates, y, o, labeldetrended)plt.plot(dates, filtered, labelfiltered)plt.plot(dates, p[0] * np.sin(2 * np.pi * dates * p[1] p[2]) p[3], ^, labelfit)fig.autofmt_xdate()plt.legend(prop{size:x-large})ax2 fig.add_subplot(212)ax2.tick_params(axisboth, whichmajor, labelsizex-large)plt.plot(f, amps, labeltransformed)plt.legend(prop{size:x-large})plt.show()7、数值积分在SciPy这个库里, 有一个用来做数值积分的包叫作scipy.quad, 不过在NumPy这里, 是没有功能完全一样的包的。高斯积分是出现在误差函数定义里的内容, 那个在数学上一般是用erf来记的, 但是呢, 高斯积分本身它的积分区间是无穷大的, 它算出来的值正好就是pi的平方根, 所以我们接下来将会使用quad这个函数去计算它具体的数值结果。print Gaussian integral, np.sqrt(np.pi),integrate.quad(lambda x: np.exp(-x**2), -np.inf, np.inf)8、插值插值是这么一码事, 就是在那些数据集里头已经给出的那些数据点之间, 去把这些空白地方填补起来, scipy里的这个函数, 是可以依靠着手头的这些实验数据来进行插值的操作, 有个类可以专门来创建什么线性插值或者是三次插值的这样的一个函数, 通常情况下不特意声明的话它会默认就是创建一个线性的插值函数, 如果想要的是三次插值函数的话那你就得去设置一下那个叫kind的参数了才能实现, 至于那个另外的类嘛它们工作的原理和方式是完全一样的一个模子, 唯一的区别在于这个类啊是用在二维的那个插值情况里面的。我们利用函数构建了一个数据集, 并且在当中加入了噪音元素, 随后我们借助于模块当中所提供的类, 分别执行了线性插值操作以及三次插值操作。import numpy as npfrom seipy import interpolateimport matplotlib.pyplot as pltx np.linspaee(-18, 18, 36)noise 0.1 * np.random.random(len(x))signal np.sinc(x) noiseinterpreted interpolate.interpld(x, signal)x2 np.linspace(-18, 18, 180)y interpreted(x2)cubic interpolate.interpld(x, signal, kindcubic)y2 cubic(x2)plt.plot(x, signal, o, labeldata)plt.plot(x2, y, -, labellinear)plt.plot(x2, y2, -, lw2, labelcubic )plt.legend()plt.show()9、图像处理我们使用了scipy模块对Lena图像进行了一些处理。from scipy import miscimport numpy as npimport matplotlib.pyplot as pltfrom scipy import ndimageimage misc.lena().astype(np.float32)plt.subplot(221)plt.title(Original Image)img plt.imshow(image, cmapplt.cm.gray)plt.axis(off)plt.subplot(222)plt.title(Median Filter)filtered ndimage.median_filter(image, size(42,42))plt.imshow(filtered, cmapplt.cm.gray)plt.axis(off )plt.subplot(223)plt.title(Rotated)rotated ndimage.rotate(image, 90)plt.imshow(rotated, cmapplt.cm.gray)plt.axis(off)plt.subplot(224)plt.title(Prewitt Filter)filtered ndimage.prewitt(image)plt.imshow(filtered, cmapplt.cm.gray)plt.axis(off)plt.show()10、对音频进行加工和处理。我们读入了那个音频片段, 然后将它重复了四遍, 最后把那个新的数组写到了一个新的WAV文件里面去了。from scipy.io import wavfileimport matplotlib.pyplot as pltimport urllib2import numpy as npimport sysresponse urllib2.urlopen(http://www.thesoundarchive.com/austinpowers/smashingbaby.wav)print response.info()WAV_FILE smashingbaby.wavfilehandle open(WAV_FILE, w)filehandle.write(response.read())filehandle.close()sample_rate, data wavfile.read(WAV_FILE)print Data type, data.dtype, Shape, data.shapeplt.subplot(2, 1, 1)plt.title(Original )plt.plot(data)plt.subplot(2, 1, 2)plt.title(Repeated)plt.plot(repeated)wavfile.write(repeated_yababy.wav, sample_rate, repeated)plt.show ()