相談事例: Python による最適化（最小二乗法を例に）

前回の吸着がらみで、具体的な最適化プログラムについて相談がありました。
自分で１からプログラムをするのは大変なので、簡単な例（toy example）を渡して、それを自分の問題に合わせて変更する、という形を取ることにしました. 以下、Python (scipy.optimize) を使ったプログラム例です.

# 最小二乗法については, 関数を与えて直接フィッティングする関数もあるかもしれない.
# ここに示す方法は, 最小二乗法に限らず, 最適化を行う一般的な方法である.


import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize

def func(x, p):   # Fitting する関数を定義する
    return p[0] + p[1] * x**2

    
def calc_ls(p, x, y):  # 最小化する関数を定義する. 最初の引数は　fitting する p とする.

    return np.sum((y - func(x, p))**2)
    
x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 3, 1, 6, 4])
p_0 = [1, 1]  # p の初期値を与える.
#
# 最小化する関数の引数を, fitting する p と その他の (x, y) に分け,  
# 最小化する関数名, fitting する p の初期値 p_0, その他の引数を args=() の形で与える.
# 問題によっては, オプションとして最適化手法を指定することを検討する (method='')
#
res = scipy.optimize.minimize(calc_ls, p_0, args=(x, y)) 
print(res)  # 結果のプリントアウト. 結果には様々な情報が含まれている.
p_est = res.x  # fitting された p　は, 結果のオブジェクトの　x という属性に入っている.
print(p_est)  # fitting された p　のプリントアウト.
#
#
# 以下はおまけとして, グラフにしてみた.
#
fig, ax = plt.subplots()
xs = np.linspace(x.min()-1, x.max()+1)
ys = func(xs, p_est)
ax.scatter(x, y)
ax.plot(xs, ys)