1. 最小二乗フィッティング#
仮定として、一組の実験データ (x [i], y [i]) があり、それらの間の関数関係 =
f (x) が分かっているとします。これらの既知の情報を通じて、関数内のいくつかのパラメータを特定する必要があります。例えば、f が線形関数 f (x) =
k*x+b である場合、パラメータ k と b は私たちが特定する必要がある値です。これらのパラメータを p で表すと、次の式の S 関数を最小にする p の値の組を見つける必要があります：
このアルゴリズムは最小二乗フィッティング（Least-square fitting）と呼ばれます。
既知のデータセット（x，y）を使用して、実際の y 値から予測された y 値の平方和を最小にするパラメータ（k，b）を特定します。
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()
この例では、フィッティングする関数は正弦波関数で、3 つのパラメータ A, k, theta があり、それぞれ振幅、周波数、位相角に対応します。私たちの実験データはノイズを含むデータ x,
y1 であり、y1 は実データ y0 にノイズを加えたものです。
leastsq 関数を使用してノイズのある実験データ x, y1 をフィッティングすることで、x と実データ y0 の間の正弦関係の 3 つのパラメータ：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のパワーを計算します
fminはこのパワーを最小にするように計算します
"""
return np.sum((yn - np.convolve(x, h))**2)
# fmin関数を呼び出し、x0を初期値として使用します
h0 = fminfunc(convolve_func, x0)
print fminfunc.__name__
print "---------------------"
# x (*) h0とyの間の相対誤差を出力します
print "yの誤差:", np.sum((np.convolve(x, h0)-y)**2)/np.sum(y**2)
# h0とhの間の相対誤差を出力します
print "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 は fsolve に fprime パラメータを介して渡され、関数 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の隣接する2点間の距離、つまり矩形の幅です
>>> y = half_circle(x) # 各点の縦座標、各点から次の点までの高さを一様と見なします。したがって、各点のy値にdxを掛けることで、現在の点から次の点までの矩形の面積を得ます。
>>> dx * np.sum(y) * 2 # 面積の2倍を求めてpiの数値を確認します
3.1412751679988937
上記の方法で計算された円上の一連の点の座標を使用して、numpy.trapz を使用して数値積分を行うこともできます：
>>> import numpy as np
>>> np.trapz(y, x) * 2 # 面積の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