Zenn
Open14

測定データをフィットしたい

Python
pandas
numpy
scipy
shotakahashotakaha

曲線フィット

scipy.optimize.curve_fit

import pandas as pd
import numpy as np
from scipy.optimize import curve_fit

f = フィット関数
xdata = X軸
ydata = Y軸
p0 = 初期パラメーター
popt, pcov = curve_fit(f, xdata, ydata, p0)
perr = np.sqrt(np.diag(pcov))
shotakahashotakaha

曲線フィットのオプション

  1. f, xdata, ydata
  2. p0=None
  3. sigma=None
  4. absolute_sigma=False
  5. check_finite=None
  6. bounds=(-inf, inf)
  7. method=None
  8. jac=None
  9. *
  10. full_output=False
  11. nan_policy=None
  12. **kwargs
shotakahashotakaha

p0: 初期パラメーター

  • フィットに使う関数の数だけ、初期パラメーターを用意する
  • p0 = Noneの場合、各要素が1の初期値が適用される
  • p0に指定した値は、numpy.atleast_1dで1次元配列（ndarray）に変換される
shotakahashotakaha

bounds: パラメーターの上限と下限

  • デフォルトは上限／下限なし（(-inf, inf)）
  • 2要素のタプル or scipy.optimize.Boundsオブジェクト
  • p0 == Noneの場合は、適当な値を設定する
if p0 is None:
        p0 = _initialize_feasible(lb, ub)
shotakahashotakaha

method: フィットの方法

  • デフォルトはlm
  • 上限／下限がある場合はtrf 
bounded_problem = np.any((lb > -np.inf) | (ub < np.inf))
if method is None:
    if bounded_problem:
        method = 'trf'
    else:
        method = 'lm'
  • 詳細はscipy.optimize.least_squaresを参照
  • trf : Trust Region Reflective algorithm
    • particularly suitable for large sparse problems with bounds. Generally robust method.
  • dogbox : dogleg algorithm with rectangular trust regions
    • typical use case is small problems with bounds
    • Not recommended for problems with rank-deficient Jacobian.
  • lm : Levenberg-Marquardt algorithm as implemented in MINPACK.
    • Doesn’t handle bounds and sparse Jacobians.
    • Usually the most efficient method for small unconstrained problems.
shotakahashotakaha

lm

if method == 'lm':
        # if ydata.size == 1, this might be used for broadcast.
        if ydata.size != 1 and n > ydata.size:
            raise TypeError(f"The number of func parameters={n} must not"
                            f" exceed the number of data points={ydata.size}")
        res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
        popt, pcov, infodict, errmsg, ier = res
        ysize = len(infodict['fvec'])
        cost = np.sum(infodict['fvec'] ** 2)
        if ier not in [1, 2, 3, 4]:
            raise RuntimeError("Optimal parameters not found: " + errmsg)
shotakahashotakaha

trf

else:
        # Rename maxfev (leastsq) to max_nfev (least_squares), if specified.
        if 'max_nfev' not in kwargs:
            kwargs['max_nfev'] = kwargs.pop('maxfev', None)

        res = least_squares(func, p0, jac=jac, bounds=bounds, method=method,
                            **kwargs)

        if not res.success:
            raise RuntimeError("Optimal parameters not found: " + res.message)

        infodict = dict(nfev=res.nfev, fvec=res.fun)
        ier = res.status
        errmsg = res.message

        ysize = len(res.fun)
        cost = 2 * res.cost  # res.cost is half sum of squares!
        popt = res.x

        # Do Moore-Penrose inverse discarding zero singular values.
        _, s, VT = svd(res.jac, full_matrices=False)
        threshold = np.finfo(float).eps * max(res.jac.shape) * s[0]
        s = s[s > threshold]
        VT = VT[:s.size]
        pcov = np.dot(VT.T / s**2, VT)
  • scipy.optimize.least_squaresでフィット
    • Solve a nonlinear least-squares problem with bounds on the variables.
  • scipy.optimise.leastsqとは別物？
    • こちらのほうが新しいインターフェース？
shotakahashotakaha

check_finite: 入力値のチェック

  • check_finite=Trueにすると、入力した配列にinfやnanがないかを確認する
    • infやnanがある場合はValueErrorになる
  • nan_policyを設定していない場合ははcheck_finite=Trueになる
shotakahashotakaha

nan_polity: NaNの扱い方

  • raise : エラー
  • omit : NaN値を除外して計算する
  • None: なにもしない
  • nan_police=None以外にすると、check_finite=Falseになる
shotakahashotakaha

sigma: ydataの誤差

  • データの誤差を考慮してフィットする（ということ？）
 if sigma is not None:
        sigma = np.asarray(sigma)
  • None or scalar or M-length sequence or MxM array
  • 残差（residuals）: r = ydata - f(xdata, *popt)（実データと推定値の差）
  • 1次元の場合：固定値 or ydataの標準偏差
# 1-D scalar
transform = 1.0 / sigma
func = _lightweight_memoizer(_wrap_func(f, xdata, ydata, transform))
  • 2次元の場合：ydataの誤差の共分散行列
# 2-D scalar
# scipy.linalg.cholesky requires lower=True to return L L^T = A
transform = cholesky(sigma, lower=True)
func = _lightweight_memoizer(_wrap_func(f, xdata, ydata, transform))
shotakahashotakaha

absolute_sigma

  • pcovの計算で使われる
  • デフォルトはabsolute_sigma=False 
warn_cov = False
if pcov is None or np.isnan(pcov).any():
    # indeterminate covariance
    pcov = zeros((len(popt), len(popt)), dtype=float)
    pcov.fill(inf)
    warn_cov = True
elif not absolute_sigma:
    if ysize > p0.size:
        s_sq = cost / (ysize - p0.size)
        pcov = pcov * s_sq
    else:
        pcov.fill(inf)
        warn_cov = True
shotakahashotakaha

full_output

  • full_output=Trueにすると、返り値が増える
if full_output:
        return popt, pcov, infodict, errmsg, ier
    else:
        return popt, pcov
  • lmの場合
 if method == 'lm':
    popt, pcov, infodict, errmsg, ier = res
  • lmじゃない場合
else:
    infodict = dict(nfev=res.nfev, fvec=res.fun)
    ier = res.status
    errmsg = res.message
    popt = res.x
    pcov = np.dot(VT.T / s**2, VT)
shotakahashotakaha

最小二乗法

  • 残差（residual）: r = y_{\text{data}} - y_{\text{fit}}
  • 残差平方和 : \sum_{\text{lower}}^{\text{upper}} r^{2}
  • 最小二乗法: \min \left( \sum r^{2} \right)