1. 最小二乘擬合#
假設有一組實驗數據 (x [i], y [i])，我們知道它們之間的函數關係 =
f (x)，通過這些已知信息，需要確定函數中的一些參數項。例如，如果 f 是一個線型函數 f (x) =
k*x+b，那麼參數 k 和 b 就是我們需要確定的值。如果將這些參數用 p 表示的話，那麼我們就是要找到一組 p 值使得如下公式中的 S 函數最小：
這種算法被稱之為最小二乘擬合 (Least-square fitting)。
用已知數據集（x，y）去確定一組參數（k，b）使得實際的 y 值減去預測的 y 值的平方和最小。
scipy 中的子函數庫 optimize 已經提供了實現最小二乘擬合算法的函數 leastsq。下面是用 leastsq 進行數據擬合的一個例子：
# -*- coding: utf-8 -*-
# 注意如果代碼裡有中文，文件開頭必須添加上面這句話，或者添加“#encoding:utf-8”。
import numpy as np
from scipy.optimize import leastsq
import pylab as pl
def func(x, p):
"""
數據擬合所用的函數: A*sin(2*pi*k*x + theta)
"""
A, k, theta = p
return A*np.sin(2*np.pi*k*x+theta)
def residuals(p, y, x):
"""
實驗數據x, y和擬合函數之間的差，p為擬合需要找到的係數
"""
return y - func(x, p)
x = np.linspace(0, -2*np.pi, 100)
A, k, theta = 10, 0.34, np.pi/6 # 真實數據的函數參數
y0 = func(x, [A, k, theta]) # 真實數據
y1 = y0 + 2 * np.random.randn(len(x)) # 加入噪聲之後的實驗數據
p0 = [7, 0.2, 0] # 第一次猜測的函數擬合參數
# 調用leastsq進行數據擬合
# residuals為計算誤差的函數
# p0為擬合參數的初始值
# args為需要擬合的實驗數據
plsq = leastsq(residuals, p0, args=(y1, x))
print u"真實參數:", [A, k, theta]
print u"擬合參數", plsq[0] # 實驗數據擬合後的參數
pl.plot(x, y0, label=u"真實數據")
pl.plot(x, y1, label=u"帶噪聲的實驗數據")
pl.plot(x, func(x, plsq[0]), label=u"擬合數據")
pl.legend()
pl.show()
這個例子中我們要擬合的函數是一個正弦波函數，它有三個參數 A, k, theta ，分別對應振幅、頻率、相角。假設我們的實驗數據是一組包含噪聲的數據 x,
y1，其中 y1 是在真實數據 y0 的基礎上加入噪聲的到了。
通過 leastsq 函數對帶噪聲的實驗數據 x, y1 進行數據擬合，可以找到 x 和真實數據 y0 之間的正弦關係的三個參數： A, k, theta。
最後我們可以發現擬合參數雖然和真實參數完全不同，但是由於正弦函數具有周期性，實際上擬合參數得到的函數和真實參數對應的函數是一致的。
2. 函數最小值#
optimize 庫提供了幾個求函數最小值的算法：fmin, fmin_powell, fmin_cg,
fmin_bfgs。下面的程序通過求解卷積的逆運算演示 fmin 的功能。
對於一個離散的線性時不變系統 h, 如果它的輸入是 x，那麼其輸出 y 可以用 x 和 h 的卷積表示：
y = x*h
現在的問題是如果已知系統的輸入 x 和輸出 y，如何計算系統的傳遞函數 h；或者如果已知系統的傳遞函數 h 和系統的輸出 y，如何計算系統的輸入 x。這種運算被稱為反卷積運算，是十分困難的，特別是在實際的運用中，測量系統的輸出總是存在誤差的。
下面用 fmin 計算反卷積，這種方法只能用在很小規模的數列之上，因此沒有很大的實用價值，不過用來評價 fmin 函數的性能還是不錯的。
# -*- coding: utf-8 -*-
# 本程序用各種fmin函數求卷積的逆運算
import scipy.optimize as opt
import numpy as np
def test_fmin_convolve(fminfunc, x, h, y, yn, x0):
"""
x (*) h = y, (*)表示卷積
yn為在y的基礎上添加一些干擾噪聲的結果
x0為求解x的初始值
"""
def convolve_func(h):
"""
計算 yn - x (*) h 的power
fmin將通過計算使得此power最小
"""
return np.sum((yn - np.convolve(x, h))**2)
# 調用fmin函數，以x0為初始值
h0 = fminfunc(convolve_func, x0)
print fminfunc.__name__
print "---------------------"
# 輸出 x (*) h0 和 y 之間的相對誤差
print "error of y:", np.sum((np.convolve(x, h0)-y)**2)/np.sum(y**2)
# 輸出 h0 和 h 之間的相對誤差
print "error of h:", np.sum((h0-h)**2)/np.sum(h**2)
print
def test_n(m, n, nscale):
"""
隨機產生x, h, y, yn, x0等數列，調用各種fmin函數求解b
m為x的長度, n為h的長度, nscale為干擾的強度
"""
x = np.random.rand(m)
h = np.random.rand(n)
y = np.convolve(x, h)
yn = y + np.random.rand(len(y)) * nscale
x0 = np.random.rand(n)
test_fmin_convolve(opt.fmin, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_powell, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_cg, x, h, y, yn, x0)
test_fmin_convolve(opt.fmin_bfgs, x, h, y, yn, x0)
if __name__ == "__main__":
test_n(200, 20, 0.1)
if __name__ == "__main__"
上面這句代碼的功能是在當前文件下編譯運行的時候，執行 if 下面的語句。如果在其他文件內通過 import 調用本文件，if 內語句不執行。
3. 非線性方程組求解#
optimize 庫中的 fsolve 函數可以用來對非線性方程組進行求解。它的基本調用形式如下：
fsolve(func, x0)
func (x) 是計算方程組誤差的函數，它的參數 x 是一個矢量，表示方程組的各個未知數的一組可能解，func 返回將 x 代入方程組之後得到的誤差；x0 為未知數矢量的初始值。如果要對如下方程組進行求解的話：
f1(u1,u2,u3) = 0
f2(u1,u2,u3) = 0
f3(u1,u2,u3) = 0
那麼 func 可以如下定義：
def func(x):
u1,u2,u3 = x
return [f1(u1,u2,u3), f2(u1,u2,u3), f3(u1,u2,u3)]
下面是一個實際的例子，求解如下方程組的解：
5*x1 + 3 = 0
4*x0*x0 - 2*sin(x1*x2) = 0
x1*x2 - 1.5 = 0
程序如下：
from scipy.optimize import fsolve
from math import sin,cos
def f(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
5*x1+3,
4*x0*x0 - 2*sin(x1*x2),
x1*x2 - 1.5
]
result = fsolve(f, [1,1,1])
print result
print f(result)
輸出為：
[-0.70622057 -0.6 -2.5 ]
[0.0, -9.1260332624187868e-14, 5.3290705182007514e-15]
由於 fsolve 函數在調用函數 f 時，傳遞的參數為數組，因此如果直接使用數組中的元素計算的話，計算速度將會有所降低，因此這裡先用 float 函數將數組中的元素轉換為 Python 中的標準浮點數，然後調用標準 math 庫中的函數進行運算。
在對方程組進行求解時，fsolve 會自動計算方程組的雅可比矩陣，如果方程組中的未知數很多，而與每個方程有關的未知數較少時，即雅可比矩陣比較稀疏時，傳遞一個計算雅可比矩陣的函數將能大幅度提高運算速度。
雅可比矩陣
雅可比矩陣是一階偏導數以一定方式排列的矩陣，它給出了可微分方程與給定點的最優線性逼近，因此類似於多元函數的導數。例如前面的函數 f1,f2,f3 和未知數 u1,u2,u3 的雅可比矩陣如下：
使用雅可比矩陣的 fsolve 實例如下，計算雅可比矩陣的函數 j 通過 fprime 參數傳遞給 fsolve，函數 j 和函數 f 一樣，有一個未知數的解矢量參數 x，函數 j 計算非線性方程組在矢量 x 點上的雅可比矩陣。由於這個例子中未知數很少，因此程序計算雅可比矩陣並不能帶來計算速度的提升。
# -*- coding: utf-8 -*-
from scipy.optimize import fsolve
from math import sin,cos
def f(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
5*x1+3,
4*x0*x0 - 2*sin(x1*x2),
x1*x2 - 1.5
]
def j(x):
x0 = float(x[0])
x1 = float(x[1])
x2 = float(x[2])
return [
[0, 5, 0],
[8*x0, -2*x2*cos(x1*x2), -2*x1*cos(x1*x2)],
[0, x2, x1]
]
result = fsolve(f, [1,1,1], fprime=j)
print result
print f(result)
4. 數值積分#
數值積分是對定積分的數值求解，例如可以利用數值積分計算某個形狀的面積。下面讓我們來考慮一下如何計算半徑為 1 的半圓的面積，根據圓的面積公式，其面積應該等於 PI/2。單位半圓（X 軸上面的半圓）曲線可以用下面的函數表示：
def half_circle(x):
return (1-x**2)**0.5
下面的程序使用經典的分小矩形計算面積總和的方式，計算出單位半圓的面積：
>>> N = 10000
>>> x = np.linspace(-1, 1, N)
>>> dx = 2.0/N # 這個值其實就是上面x每相鄰兩點之間的距離，也就是舉行的寬
>>> y = half_circle(x) # 每一點的縱坐標，每一點到下一點極端的看成高一致。這樣每一點的y值乘以dx就是當前點到下一點之間矩形的面積。
>>> dx * np.sum(y) * 2 # 求面積的兩倍方便看pi數值
3.1412751679988937
利用上述方式計算出的圓上一系列點的坐標，還可以用 numpy.trapz 進行數值積分：
>>> import numpy as np
>>> np.trapz(y, x) * 2 # 面積的兩倍
3.1415893269316042
此函數計算的是以 x,y 為頂點坐標的折線與 X 軸所夾的面積。同樣的分割點數，trapz 函數的結果更加接近精確值一些。
如果我們調用 scipy.integrate 庫中的 quad 函數的話，將會得到非常精確的結果：
>>> from scipy import integrate
>>> pi_half, err = integrate.quad(half_circle, -1, 1)
>>> pi_half*2
3.1415926535897984
這一節我搬過來的內容少很多，想學更多內容的讀者建議去原地址學習。
本文參考自： http://old.sebug.net/paper/books/scipydoc/scipy_intro.html#id1