Scipy
numpy 其实是scipy的一部分，后来才从scipy中分离出来。scipy函数库在numpy库的基础上征集了众多数学、科学以及工程计算中常用的函数库。例如线性代数、常微分方程值求解、信号处理、图形处理、稀疏矩阵等等。由于其中涉及的领域众多，只能想到哪里，学到哪里
1 插值
一维插值和二维插值，是最常用的scipy功能之一，也很容易上手
1.1 一维插值和样条插值
下面的例子清楚的展示了线性插值和样条插值之后的数据形态。
import numpy as np
from scipy import interpolate
x = np.arange(0, 10)
y = np.exp(-x/3.0) x array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]) y array([1. , 0.71653131, 0.51341712, 0.36787944, 0.26359714,
0.1888756 , 0.13533528, 0.09697197, 0.06948345, 0.04978707]) f_linear = interpolate.interp1d(x, y)
x_new = np.arange(0,9,0.2)
f_linear(x_new) # 线性插值 array([1. , 0.94330626, 0.88661252, 0.82991879, 0.77322505,
0.71653131, 0.67590847, 0.63528563, 0.5946628 , 0.55403996,
0.51341712, 0.48430958, 0.45520205, 0.42609451, 0.39698698,
0.36787944, 0.34702298, 0.32616652, 0.30531006, 0.2844536 ,
0.26359714, 0.24865283, 0.23370852, 0.21876422, 0.20381991,
0.1888756 , 0.17816754, 0.16745947, 0.15675141, 0.14604335,
0.13533528, 0.12766262, 0.11998996, 0.11231729, 0.10464463,
0.09697197, 0.09147426, 0.08597656, 0.08047886, 0.07498115,
0.06948345, 0.06554417, 0.0616049 , 0.05766562, 0.05372634]) bs = interpolate.splrep(x, y)
interpolate.splev(x_new, bs) # 样条插值 array([1. , 0.93571489, 0.8754193 , 0.8189194 , 0.76602135,
0.71653131, 0.67025545, 0.62699994, 0.58657094, 0.54877461,
0.51341712, 0.48031625, 0.44933621, 0.42035284, 0.39324198,
0.36787944, 0.34414628, 0.32194438, 0.30118082, 0.28176271,
0.26359714, 0.24659576, 0.23068846, 0.2158097 , 0.20189393,
0.1888756 , 0.17669225, 0.16529366, 0.15463272, 0.1446623 ,
0.13533528, 0.1266067 , 0.11844022, 0.11080172, 0.10365702,
0.09697197, 0.09071432, 0.0848594 , 0.07938446, 0.07426673,
0.06948345, 0.06501185, 0.06082918, 0.05691267, 0.05323955]) import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['FangSong'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图片是‘-’的显示问题
plt.rcParams['figure.dpi'] = 200 #分辨率
plt.plot(x, y, 'o', label=u'原始数据')
plt.plot(x_new, f_linear(x_new), label=u'线性插值')
plt.plot(x_new, interpolate.splev(x_new, bs), label=u"B-spline插值")
plt.legend()
plt.show() 
样条插值中包含很多的默认参数，下面简单说明
scipy.interpolate.splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None, full_output=0, per=0, quiet=1) - 参数s用来确定平滑点数，通常是m-SQRT(2m),m是曲线点数。如果在插值中不需要平滑应该设定s=0。splrep()输出的是一个3元素的元胞数组（t,c,k）,其中t是曲线点，c是计算出来的系数，k是样条阶数，通常是3阶，但可以通过k改变。
scipy.interpolate.splev(x, tck, der=0) - 其中的der是进行样条计算是需要实际计算到的阶数，必须满足条件der<=k。
1.2 二维插值
在一个房间的地板上按九宫格的位置放置9个传感器，测的温度如下：
import numpy as np
temp = np.random.randint(20, 30, (3, 3))
temp array([[26, 28, 26],
[25, 22, 23],
[28, 29, 25]]) 在不增加传感器的前提下是数据变得平滑
x = np.arange(3)
y = np.arange(3)
ip = interpolate.interp2d(x, y, temp)
x_new = np.linspace(0, 2, 5)
y_new = np.linspace(0, 2, 5)
temp_new = ip(x_new, y_new)
temp_new array([[26. , 27. , 28. , 27. , 26. ],
[25.5 , 25.25, 25. , 24.75, 24.5 ],
[25. , 23.5 , 22. , 22.5 , 23. ],
[26.5 , 26. , 25.5 , 24.75, 24. ],
[28. , 28.5 , 29. , 27. , 25. ]]) import matplotlib.pyplot as plt
plt.subplot(121)
plt.imshow(temp)
plt.subplot(122)
plt.imshow(temp_new)
plt.show() 
2 拟合
在工作中，我们常常需要在途中描绘某些实际数据观察的同时，使用一个曲线来拟合这些实际数据。所谓拟合，计时找出符合数据变化趋势的曲线方程，或者直接绘制出拟合曲线。
2.1 使用numpy.polyfit拟合
下面的代码，基于numpy模块，可以直接绘制处拟合曲线，但无法的到曲线方程，尽管输出了一堆曲线参数。
import numpy
import pylab
def plot_polynomail_fit(y, *deg):
'''
这个函数一次拟合一组数据，但可以对一组数据同时拟合多条曲线并显示
y是一个list，为代拟和的数据
*deg是一个元组，长度不定，里面存放拟合的次数，可以对一组数据拟合出多条直线进行比较
'''
x = range(len(y))
COLOR = ['c', 'm', 'y', 'k', 'r', 'p', 'o', 'g', 'b']
temp = []
numOfLineFit = len(deg) # 需要拟合的次数列表
for index, item in enumerate(deg):
param = numpy.polyfit(x, y, item) # 曲线参数
equation = numpy.poly1d(param) # 曲线方程
temp.extend(param[:]) # 后去曲线参数
print(param)
pylab.subplot(numOfLineFit, 1, index+1)
pylab.plot(x, equation(x), '%s--' % COLOR[index])
pylab.plot(x, y, 'b--')
pylab.show()
return temp
if __name__ == "__main__":
y = [17, 19, 21, 28, 33, 38, 37, 37, 31, 23, 19, 18]
plot_polynomail_fit(y, 2, 3, 4) [-0.68406593 7.70304695 13.22802198]
[-0.01981352 -0.35714286 6.32600733 14.20879121]
[ 0.03095862 -0.70090326 4.33530012 -4.0849359 17.71153846] 
2.2 使用scipyoptimize.curve_fit拟合
scipy提供的拟合，貌似需要先确定带参数的曲线方程，然后由scipy求解方程，返回曲线参数。我们试试 以上面的一组数据为例使用scipy拟合曲线。
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
x = np.arange(1, 13, 1)
y = np.array([17,19,21,28,33,38,37,37,31,23,19,18])
plt.plot(x, y)
plt.show() 
可以看出，曲线近似正弦函数。构建函数y=asin(xpi/6+b)+c, 使用scipy的optimize.cure_fit可以求出a,b,c的值：
def fmax(x, a, b, c):
return a*np.sin(x*np.pi/6+b)+c
fita, fitb = optimize.curve_fit(fmax, x, y, [1, 1, 1])
xn = np.arange(1, 13, 0.1)
plt.plot(x, y)
plt.plot(xn, fmax(xn, fita[0], fita[1], fita[2]))
plt.show() 
3 求解非线性方程组
在数学建模中，需要对一些稀奇古怪的方程（组）求解， matlab自然是首选，但Matlab不是免费的，scipy则为我们提供了免费的午餐！scipy.optimize库中的fsolve函数可以用来对非线性方程组（组）进行求解。它的基本形式如下：fsolve(func, x0)
func(x)是计算方程组误差的函数，它的参数x是一个矢量，表示方程组的各个位置数的一组可能的解，func返回将x代入方程组之后的得到的误差；x0为未知数矢量的初始值。
我们先来求解一个简单的方程：$\sin(x)-\cos(x)=0.2$
from scipy.optimize import fsolve
import numpy as np
def f(A):
x = float(A[0])
return [np.sin(x) - np.cos(x) - 0.2]
result = fsolve(f, [1])
print("结果：", result, "误差：", f(result)) 结果： [0.92729522] 误差： [-1.6653345369377348e-16] So easy. Now let’s try some difficult
$4x^2−2sin(yz)=0$
$5y+3=0 $
$yz−1.5=0$
from scipy.optimize import fsolve
import numpy as np
def f(A): # 误差函数
x = float(A[0])
y = float(A[1])
z = float(A[2])
return [4*x*x-2*np.sin(y*z), 5*y+3, y*z-1.5] # 返回各个式子对应的误差。
result = fsolve(f, [1, 1, 1]) # 对应的各个未知数的解
deviation = f(result)
print("结果: %s\n误差：%s" % (str(result), str(deviation))) 结果: [-0.70622057 -0.6 -2.5 ]
误差：[-9.126033262418787e-14, 0.0, 5.329070518200751e-15] 4 数值积分
数值积分是对定积分的数值求解，例如可以利用数值积分求解某个形状的面积。我们知道，圆的方程可写成：$x^2+y^2=1$， 下面让我们来考虑一下如何计算半径为1的半圆的面积，根据圆的面积公式，其面积应等于Π。单位半圆曲线可用以下函数表示$y=\sqrt{(1-x^2)}$
我们先定义一个根据x计算y的函数：
def half_circle(x):
return (1 - x**2) ** 0.5 4.1 经典微分法
下面的程序使用经典的分小矩形的方法计算面积总和的方式，计算出半圆的面积：
N = 10000
x = np.linspace(-1, 1, N)
dx = 2.0 / N
y = half_circle(x)
dx * np.sum(y[:-1] + y[1:]) # 面积的两倍 3.1412751679989044 4.2 使用定积分求解函数
如果我们调用scipy.integrate库中的quad函数的话将会的到非常精确的结果
from scipy import integrate
pi_half, err = integrate.quad(half_circle, -1, 1)
print("圆的面积：%f\n误差：%f" % (pi_half*2, err)) 圆的面积：3.141593
误差：0.000000 5 图形处理
在scipy.misc模块中，有一个函数可以载入Ascent幅图像是被用作图像处理的经典示范图。简单演示在该图像上的操作。
- 载入Ascent图像，并显示灰度图像
- 应用中值滤波扫描每一个数据点，并替换为相邻数据点的中值
- 旋转图像
- 应用Prwitt滤波器（基于图形强度的梯度计算）
from scipy import misc
from scipy import ndimage
img = misc.ascent().astype(np.float32)
plt.subplot(221)
plt.title("Original Image")
plt.imshow(img, cmap=plt.cm.gray) # 1
plt.subplot(222)
plt.title("Median Filter")
filtered = ndimage.median_filter(img, size=(42, 42))
plt.imshow(filtered, cmap=plt.cm.gray) # 2
plt.subplot(223)
plt.title("Rotated")
rotated = ndimage.rotate(img, 90)
plt.imshow(rotated, cmap=plt.cm.gray)# 3
plt.subplot(224)
plt.title("Prewitt Filter")
filtered = ndimage.prewitt(img)
plt.imshow(filtered, cmap=plt.cm.gray) # 4
plt.show() 
写个博客不容易，可怜可怜博主，点个广告再走呗(✿◕‿◕✿)。
