数値計算 - Intro to Python

目次

重み付き非線形最小二乗法を行う

scipyパッケージのleast_squares関数を使う。第一引数に残差を返す関数を指定するが、重み無しの場合は「測定値-回帰モデルによる計算値」を返す関数を、重み付きの場合は「重み×（測定値-回帰モデルによる計算値）」を返す関数を指定すればよい。

以下の例では、与えられた8組の測定値を使用して、非線形最小二乗法により回帰モデルy = b1 + b2 * x * exp(b3 * x)（b1, b2, b3は回帰係数）を作成した例。重み無し、重み付きでそれぞれ計算している。図を見てのとおり、重みを考慮しないモデル（青実線）は全点の近くをまんべんなく通っているが、重みを考慮したモデル（x = 3, 6のみ重みが他の100倍、赤実線）は、x = 3, 6の測定値を無理矢理通るようなモデル（回帰係数）が求められていることがわかる。

>>> import numpy as np
>>> import math
>>> from scipy.optimize import least_squares
>>> import matplotlib.pyplot as plt
>>>
>>> # 計算の準備
>>> x = np.array([1, 2, 3, 5, 6, 7, 8, 9])
>>> y = np.array([1, 4, 10, 25, 32, 49, 64, 75])
>>> w = np.array([1] * len(x))
>>> w[2] = 100
>>> w[4] = 100
>>>
>>> # 回帰モデル
>>> def f(x, b):
...     return b[0] + b[1] * x * np.exp(b[2] * x)
...
>>> # 回帰係数の初期値
>>> b0 = [1, 1, 1]
>>>
>>> # 残差平方和を計算する関数（重み無し）
>>> def res1(b, x, y):
...     return y - f(x, b)
...
>>> # 残差平方和を計算する関数（重み付き）
>>> def res2(b, x, y, w):
...     return np.sqrt(w) * (y - f(x, b))
...
>>> # 重みを考慮しない計算（図の青実線）
>>> r1 = least_squares(res1, b0, args = (x, y))
>>> print(r1.x)  # 推定した回帰係数
[-3.89077136  3.39493809  0.1084658 ]
>>>
>>> # 重みを考慮した計算（図の赤実線）
>>> r2 = least_squares(res2, b0, args = (x, y, w))
>>> print(r2.x)  # 推定した回帰係数
[-0.0304514   2.06079139  0.15929171]
>>>
>>> plt.plot(x, f(x, r1.x), '-', color = 'blue')
>>> plt.plot(x, f(x, r2.x), '-', color = 'red')
>>> plt.scatter(x, y, marker = 'o', color = 'black')
>>> plt.show()