NumPy的扩展：SciPy

最新推荐文章于 2026-03-23 07:00:00 发布

原创最新推荐文章于 2026-03-23 07:00:00 发布 · 603 阅读

3 ·

本内容遵循CC 4.0 BY-SA版权协议

标签

#Scipy #python

Python 专栏收录该内容

15 篇文章

订阅专栏

本文介绍了SciPy作为Python科学计算库的重要扩展，涵盖了保存和加载.mat文件、统计分析、信号处理、趋势分析、傅里叶分析、数学优化、数值积分、插值、图像处理和音频处理等多个方面，提供了丰富的函数和工具，便于进行数据处理和分析。

SciPy是世界著名的Python开源科学计算库，建立在NumPy之上。它增的功能包括数值积分、最优化、统计和一些专用函数。

1、保存和加载.mat 文件

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

a = np.arange(7) 
io.savemat("a.mat", {"array": a})

2、统计

scipy.stats包中的统计函数分析生成的数据。

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

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

(4) 峰度（kurtosis）描述的是概率分布曲线的陡峭程度。我们来做一个峰度检验。该检验与
偏度检验类似，当然这里是针对峰度。
print "Kurtosistest", "pvalue", stats.kurtosistest(generated)

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

(6) 使用SciPy我们可以很方便地得到数据所在的区段中某一百分比处的数值print "95 percentile", stats.scoreatpercentile(generated, 95)

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

我们按正态分布生成了一个随机数据集，并使用scipy.stats模块分析了该数据集：

from scipy import stats 
import matplotlib.pyplot as plt 
generated = stats.norm.rvs(size=900) 
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.signal模块中包含滤波函数和B样条插值（B-spline interpolation）函数。

4、趋势分析

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 
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() 
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.fftpack模块中找到（NumPy也有自己的傅里叶工具包，即numpy.fft）。这个模块包含快速傅里叶变换、微分算子和拟微分算子以及一些辅助函数。MATLAB用户会很高兴，因为scipy.fftpack模块中的很多函数与MATLAB对应的函数同名，且功能也很相近。

去除了一个信号的趋势，并使用scipy.fftpack模块对其应用了一个简单的滤波器。

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 scipy import fftpack 
from matplotlib.dates import DateFormatter 
from matplotlib.dates import DayLocator 
from matplotlib.dates import MonthLocator 
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() 
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(axis='both', which='major', labelsize='x-large') 
amps = np.abs(fftpack.fftshift(fftpack.rfft(y))) 
amps[amps < 0.1 * amps.max()] = 0 
plt.plot(dates, y, 'o', label="detrended") 
plt.plot(dates, -fftpack.irfft(fftpack.ifftshift(amps)), label="filtered") 
fig.autofmt_xdate() 
plt.legend(prop={'size':'x-large'}) 
ax2 = fig.add_subplot(212) 
ax2.tick_params(axis='both', which='major', labelsize='x-large') 
N = len(qqq) 
plt.plot(np.linspace(-N/2, N/2, N), amps, label="transformed") 
plt.legend(prop={'size':'x-large'}) 
plt.show()

6、数学优化

优化算法（optimization algorithm）尝试寻求某一问题的最优解，例如找到函数的最大值或最小值，函数可以是线性或者非线性的。解可能有一些特定的约束，例如不允许有负数。在scipy.optimize模块中提供了一些优化算法，最小二乘法函数leastsq就是其中之一。当调用这个函数时，我们需要提供一个残差（误差项）函数。这样，leastsq将最小化残差的平方和。得到的解与我们使用的数学模型有关。我们还需要为算法提供一个起始点，这应该是一个最好的猜测——尽可能接近真实解。否则，程序执行800轮迭代后将停止。

使用scipy.optimize模块对滤波后的信号拟合了一个正弦波函数。

from matplotlib.finance import quotes_historical_yahoo 
import numpy as np 
import matplotlib.pyplot as plt 
from scipy import fftpack 
from scipy import signal 
from matplotlib.dates import DateFormatter 
from matplotlib.dates import DayLocator 
from matplotlib.dates import MonthLocator 
from scipy import optimize 
start = (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(axis='both', which='major', labelsize='x-large') 
amps = np.abs(fftpack.fftshift(fftpack.rfft(y))) 
amps[amps < amps.max()] = 0 
def residuals(p, y, x): 
  A,k,theta,b = p 
  err = y-A * np.sin(2* np.pi* k * x + theta) + b 
  return err 
filtered = -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", p0 
plsq = optimize.leastsq(residuals, p0, args=(filtered, dates)) 
p = plsq[0] 
print "P", p 
plt.plot(dates, y, 'o', label="detrended") 
plt.plot(dates, filtered, label="filtered") 
plt.plot(dates, p[0] * np.sin(2 * np.pi * dates * p[1] + p[2]) + p[3], '^', label="fit") 
fig.autofmt_xdate() 
plt.legend(prop={'size':'x-large'}) 
ax2 = fig.add_subplot(212) 
ax2.tick_params(axis='both', which='major', labelsize='x-large') 
plt.plot(f, amps, label="transformed") 
plt.legend(prop={'size':'x-large'}) 
plt.show()

7、数值积分

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

print "Gaussian integral", np.sqrt(np.pi), 
integrate.quad(lambda x: np.exp(-x**2), -np.inf, np.inf)

8、插值

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

我们用sinc函数创建了一个数据集并加入了噪音，然后使用scipy.interpolate模块中的interp1d类进行了线性插值和三次插值。

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 
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()

9、图像处理

我们使用scipy.ndimage模块对Lena图像进行了一些处理。

from scipy import misc 
import numpy as np 
import matplotlib.pyplot as plt 
from scipy import ndimage 
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") 
plt.subplot(224) 
plt.title("Prewitt Filter") 
filtered = ndimage.prewitt(image) 
plt.imshow(filtered, cmap=plt.cm.gray) 
plt.axis("off") 
plt.show()

10、音频处理

我们读入了一个音频片段将其重复四遍并将新数组写入了一个新的WAV文件。

from scipy.io import wavfile 
import matplotlib.pyplot as plt 
import urllib2 
import numpy as np 
import sys 
response = urllib2.urlopen('http://www.thesoundarchive.com/austinpowers/ 
smashingbaby.wav') 
print response.info() 
WAV_FILE = 'smashingbaby.wav' 
filehandle = open(WAV_FILE, 'w') 
filehandle.write(response.read()) 
filehandle.close() 
sample_rate, data = wavfile.read(WAV_FILE) 
print "Data type", data.dtype, "Shape", data.shape 
plt.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 ()