Scipy是一個用于數學、科學、工程領域的常用軟件包，可以處理插值、積分、優化、圖像處理、常微分方程數值解的求解、信號處理等問題。它用于有效計算Numpy矩陣，使Numpy和Scipy協同工作，高效解決問題。

Scipy 由不同科學計算領域的子模塊組成：

子模塊描述

cluster	聚類算法
constants	物理數學常數
fftpack	快速傅里葉變換
integrate	積分和常微分方程求解
interpolate	插值
io	輸入輸出
linalg	線性代數
odr	正交距離回歸
optimize	優化和求根
signal	信號處理
sparse	稀疏矩陣
spatial	空間數據結構和算法
special	特殊方程
stats	統計分布和函數
weave	C/C++ 積分

在使用 Scipy 之前，為了方便，假定這些基礎的模塊已經被導入：

import numpy as np import scipy as sp import matplotlib as mpl import matplotlib.pyplot as plt

由于 Scipy 以 Numpy 為基礎，因此很多基礎的 Numpy 函數可以在scipy 命名空間中直接調用

1.插值

樣條插值法是一種以可變樣條來作出一條經過一系列點的光滑曲線的數學方法

from scipy.interpolate import interp1d np.set_printoptions(precision=2, suppress=True) #設置 Numpy 浮點數顯示格式 #從文本中讀入數據 data = np.genfromtxt("JANAF_CH4.txt", delimiter="\t", # TAB 分隔skiprows=1, # 忽略首行names=True, # 讀入屬性missing_values="INFINITE", # 缺失值filling_values=np.inf) # 填充缺失值 ch4_cp = interp1d(data['TK'], data['Cp']) #默認情況下，輸入值要在插值允許的范圍內，否則插值會報錯

我們可以通過 kind 參數來調節使用的插值方法，來得到不同的結果:

• nearest 最近鄰插值
• zero 0階插值
• linear 線性插值
• quadratic 二次插值
• cubic 三次插值
• 4,5,6,7 更高階插值

對于徑向基函數，其插值的公式為：

通過數據點 xj來計算出 nj 的值，來計算 x處的插值結果

from scipy.interpolate.rbf import Rbf cp_rbf = Rbf(data['TK'], data['Cp'], function = "multiquadric") plt.plot(data['TK'], data['Cp'], 'k+') p = plt.plot(data['TK'], cp_rbf(data['TK']), 'r-')

高維 RBF 插值:

from mpl_toolkits.mplot3d import Axes3D x, y = np.mgrid[-np.pi/2:np.pi/2:5j, -np.pi/2:np.pi/2:5j] z = np.cos(np.sqrt(x**2 + y**2))zz = Rbf(x, y, z) xx, yy = np.mgrid[-np.pi/2:np.pi/2:50j, -np.pi/2:np.pi/2:50j] fig = plt.figure(figsize=(12,6)) ax = fig.gca(projection="3d") ax.plot_surface(xx,yy,zz(xx,yy),rstride=1, cstride=1, cmap=plt.cm.jet)

2.概率統計方法

Scipy 中的子庫 scipy.stats 中包含很多統計上的方法

Numpy 自帶簡單的統計方法:mean(),min(),max(),std()

import scipy.stats.stats as st print 'median, ', st.nanmedian(heights) # 忽略nan值之后的中位數 print 'mode, ', st.mode(heights) # 眾數及其出現次數 print 'skewness, ', st.skew(heights) # 偏度 print 'kurtosis, ', st.kurtosis(heights) # 峰度

（1）連續分布，以正態分布為例：

from scipy.stats import norm

它包含四類常用的函數：

• norm.cdf 返回對應的累計分布函數值
• norm.pdf 返回對應的概率密度函數值
• norm.rvs 產生指定參數的隨機變量
• norm.fit 返回給定數據下，各參數的最大似然估計（MLE）值

#產生正態分布，通過 loc 和 scale 來調整參數 p = plot(x, norm.pdf(x, loc=0, scale=1)) p = plot(x, norm.pdf(x, loc=0.5, scale=2)) p = plot(x, norm.pdf(x, loc=-0.5, scale=.5))

（2）離散分布，沒有概率密度函數，但是有概率質量函數?PMF

from scipy.stats import binom, poisson, randint

二項分布:

num_trials = 60 x = arange(num_trials)plot(x, binom(num_trials, 0.5).pmf(x), 'o-', label='p=0.5') plot(x, binom(num_trials, 0.2).pmf(x), 'o-', label='p=0.2')legend()

自定義離散分布:

from scipy.stats import rv_discrete xk = [1, 2, 3, 4, 5, 6] pk = [.3, .35, .25, .05, .025, .025] #定義離散分布,此時， loaded 可以當作一個離散分布的模塊來使用 loaded = rv_discrete(values=(xk, pk))#產生100個隨機變量，將直方圖與概率質量函數進行比較 samples = loaded.rvs(size=100) bins = linspace(.5,6.5,7)hist(samples, bins=bins, normed=True) stem(xk, loaded.pmf(xk), markerfmt='ro', linefmt='r-')

(3)舉例：

配對樣本t檢驗：指的是兩組樣本之間的元素一一對應

例如，假設我們有一組病人的數據

pop_size = 35pre_treat = norm(loc=0, scale=1) n0 = pre_treat.rvs(size=pop_size)

經過某種治療后，對這組病人得到一組新的數據：

effect = norm(loc=0.05, scale=0.2) eff = effect.rvs(size=pop_size)n1 = n0 + eff

新數據的最大似然估計：

loc, scale = norm.fit(n1) post_treat = norm(loc=loc, scale=scale) fig = figure(figsize=(10,4))ax1 = fig.add_subplot(1,2,1) h = ax1.hist([n0, n1], normed=True) p = ax1.plot(x, pre_treat.pdf(x), 'b-') p = ax1.plot(x, post_treat.pdf(x), 'g-')ax2 = fig.add_subplot(1,2,2) h = ax2.hist(eff, normed=True)

t_val, p = ttest_rel(n0, n1)print 't = {}'.format(t_val) print 'p-value = {}'.format(p) t = -1.89564459709 p-value = 0.0665336223673

配對 t 檢驗的結果說明，配對樣本之間存在顯著性差異，說明治療時有效的，符合我們的預期。

3.曲線擬合

（1）多項式：

#導入線多項式擬合工具： from numpy import polyfit, poly1dx = np.linspace(-5, 5, 100) y = 4 * x + 1.5 noise_y = y + np.random.randn(y.shape[-1]) * 2.5coeff = polyfit(x, noise_y, 1)p = plt.plot(x, noise_y, 'rx') p = plt.plot(x, coeff[0] * x + coeff[1], 'k-') p = plt.plot(x, y, 'b--')#還可以用 poly1d 生成一個以傳入的 coeff 為參數的多項式函數： #f = poly1d(coeff) #p = plt.plot(x, noise_y, 'rx') #p = plt.plot(x, f(x))

（2）最小二乘擬合

當我們使用一個 N-1 階的多項式擬合這 M 個點時，有這樣的關系存在：XC=Y

?

from scipy.linalg import lstsq from scipy.stats import linregressx = np.linspace(0,5,100) y = 0.5 * x + np.random.randn(x.shape[-1]) * 0.35 X = np.hstack((x[:,np.newaxis], np.ones((x.shape[-1],1)))) C, resid, rank, s = lstsq(X, y) p = plt.plot(x, y, 'rx') p = plt.plot(x, C[0] * x + C[1], 'k--')

slope, intercept, r_value, p_value, stderr = linregress(x, y) p = plt.plot(x, slope * x + intercept, 'k--')

還有許多高級擬合方式，用時自查

4.最小化函數

（1）minimize 函數

已知斜拋運動的水平飛行距離公式：

d=2v20gsin(θ)cos(θ)d=2v02gsin?(θ)cos?(θ)

• dd 水平飛行距離
• v0v0 初速度大小
• gg 重力加速度
• θθ 拋出角度

希望找到使 dd 最大的角度 θ

def dist(theta, v0):"""calculate the distance travelled by a projectile launchedat theta degrees with v0 (m/s) initial velocity."""g = 9.8theta_rad = pi * theta / 180return 2 * v0 ** 2 / g * sin(theta_rad) * cos(theta_rad) theta = linspace(0,90,90)#最大化距離就相當于最小化距離的負數： def neg_dist(theta, v0):return -1 * dist(theta, v0) from scipy.optimize import minimize result = minimize(neg_dist, 40, args=(1,)) print "optimal angle = {:.1f} degrees".format(result.x[0])

minimize 接受三個參數：第一個是要優化的函數，第二個是初始猜測值，第三個則是優化函數的附加參數，默認 minimize 將優化函數的第一個參數作為優化變量，所以第三個參數輸入的附加參數從優化函數的第二個參數開始。

（2）Rosenbrock 函數

Rosenbrock 函數是一個用來測試優化函數效果的一個非凸函數：

from scipy.optimize import rosen from mpl_toolkits.mplot3d import Axes3Dx, y = meshgrid(np.linspace(-2,2,25), np.linspace(-0.5,3.5,25)) z = rosen([x,y]) fig = figure(figsize=(12,5.5)) ax = fig.gca(projection="3d") ax.azim = 70; ax.elev = 48 ax.set_xlabel("X"); ax.set_ylabel("Y") ax.set_zlim((0,1000)) p = ax.plot_surface(x,y,z,rstride=1, cstride=1, cmap=cm.jet) rosen_min = ax.plot([1],[1],[0],"ro")

不同方法的計算開銷量是不同的

5.積分

（1）符號積分

符號運算可以用 sympy 模塊完成。

from sympy import init_printing #方便其顯示 from sympy import symbols, integrate import sympyx, y = symbols('x y') sympy.sqrt(x ** 2 + y ** 2)

z = sympy.sqrt(x ** 2 + y ** 2) z.subs(x, 3) #將其中的 x 利用 subs 方法替換為 3

from sympy.abc import theta y = sympy.sin(theta) ** 2 y #對 y 進行積分 Y = integrate(y)

查看具體數值:

integrate(y, (theta, 0, sympy.pi)).evalf() integrate(y, (theta, 0, np.pi))

(2)數值積分

導入貝塞爾函數

from scipy.special import jv def f(x):return jv(2.5, x)

積分到無窮

from numpy import inf interval = [0., inf]def g(x):return np.exp(-x ** 1/2) value, max_err = quad(g, *interval) x = np.linspace(0, 10, 50) fig = plt.figure(figsize=(10,3)) p = plt.plot(x, g(x), 'k-') p = plt.fill_between(x, g(x)) plt.annotate(r"$\int_0^{\infty}e^{-x^1/2}dx = $" + "{}".format(value), (4, 0.6),fontsize=16) print "upper bound on error: {:.1e}".format(max_err)

6.解微分方程

（1）簡單一階：

from scipy.integrate import odeint def dy_dt(y, t):return np.sin(t) t = np.linspace(0, 2*pi, 100)result = odeint(dy_dt, 0, t)fig = figure(figsize=(12,4)) p = plot(t, result, "rx", label=r"$\int_{0}^{x}sin(t) dt $") p = plot(t, -cos(t) + cos(0), label=r"$cos(0) - cos(t)$") p = plot(t, dy_dt(0, t), "g-", label=r"$\frac{dy}{dt}(t)$") l = legend(loc="upper right") xl = xlabel("t")

（2）高階微分方程：

def dy_dt(y, t):"""Governing equations for projectile motion with drag.y[0] = positiony[1] = velocityg = gravity (m/s2)D = drag (1/s) = force/velocitym = mass (kg)"""g = -9.8D = 0.1m = 0.15dy1 = g - (D/m) * y[1]dy0 = y[1] if y[0] >= 0 else 0.return [dy0, dy1]position_0 = 0. velocity_0 = 100 t = linspace(0, 12, 100) y = odeint(dy_dt, [position_0, velocity_0], t) p = plot(t, y[:,0]) yl = ylabel("Height (m)") xl = xlabel("Time (s)")

7.稀疏矩陣

稀疏矩陣主要使用 位置 + 值 的方法來存儲矩陣的非零元素，根據存儲和使用方式的不同，有如下幾種類型的稀疏矩陣：

類型描述

bsr_matrix(arg1[, shape, dtype, copy, blocksize])	Block Sparse Row matrix
coo_matrix(arg1[, shape, dtype, copy])	A sparse matrix in COOrdinate format.
csc_matrix(arg1[, shape, dtype, copy])	Compressed Sparse Column matrix
csr_matrix(arg1[, shape, dtype, copy])	Compressed Sparse Row matrix
dia_matrix(arg1[, shape, dtype, copy])	Sparse matrix with DIAgonal storage
dok_matrix(arg1[, shape, dtype, copy])	Dictionary Of Keys based sparse matrix.
lil_matrix(arg1[, shape, dtype, copy])	Row-based linked list sparse matrix

在這些存儲格式中：

• COO 格式在構建矩陣時比較高效
• CSC 和 CSR 格式在乘法計算時比較高效

from scipy.sparse import * import numpy as npcoo_matrix((2,3)) #也可以使用一個已有的矩陣或數組或列表中創建新矩陣： A = coo_matrix([[1,2,0],[0,0,3],[4,0,5]])

(0, 0) 1 (0, 1) 2 (1, 2) 3 (2, 0) 4 (2, 2) 5

可以轉化為普通矩陣：

C = A.todense()

matrix([[1, 2, 0], [0, 0, 3], [4, 0, 5]])

還可以傳入一個 (data, (row, col)) 的元組來構建稀疏矩陣：

I = np.array([0,3,1,0]) J = np.array([0,3,1,2]) V = np.array([4,5,7,9]) A = coo_matrix((V,(I,J)),shape=(4,4)) print A

(0, 0) 4 (3, 3) 5 (1, 1) 7 (0, 2) 9

8.線性代數

import numpy as np import numpy.linalg import scipy as sp import scipy.linalg import matplotlib.pyplot as plt from scipy import linalg

scipy.linalg 包含 numpy.linalg 中的所有函數，同時還包含了很多 numpy.linalg 中沒有的函數，在使用時，我們一般使用 scipy.linalg 而不是 numpy.linalg

（1）求解方程

A = np.array([[1, 3, 5],[2, 5, 1],[2, 3, 8]]) b = np.array([10, 8, 3]) x = linalg.solve(A, b)

（2）計算行列式

A = np.array([[1, 3, 5],[2, 5, 1],[2, 3, 8]])linalg.det(A)

（3）矩陣分解

對于給定的 N×NN×N 矩陣 AA，特征值和特征向量問題相當與尋找標量 λλ 和對應的向量 vv 使得：

?

Av=λvAv=λv

矩陣的 NN 個特征值（可能相同）可以通過計算特征方程的根得到：

?

|A?λI|=0|A?λI|=0

然后利用這些特征值求（歸一化的）特征向量。

問題求解

• linalg.eig(A)
  • 返回矩陣的特征值與特征向量
• linalg.eigvals(A)
  • 返回矩陣的特征值
• linalg.eig(A, B)
  • 求解 Av=λBvAv=λBv 的問題

（4）對稀疏矩陣

• scipy.sparse.linalg.inv
  • 稀疏矩陣求逆
• scipy.sparse.linalg.expm
  • 求稀疏矩陣的指數函數
• scipy.sparse.linalg.norm
  • 稀疏矩陣求范數

對于特別大的矩陣，原來的方法可能需要太大的內存，考慮使用這兩個方法替代：

• scipy.sparse.linalg.eigs
  • 返回前 k 大的特征值和特征向量
• scipy.sparse.linalg.svds
  • 返回前 k 大的奇異值和奇異向量

?

?

?

?

?

?

?

?

?

?

?

?

?

總結

以上是生活随笔為你收集整理的python之scipy库详解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。