=============================================================
フィッティング３(最小二乗法，optimizeの利用)
=============================================================


.. toctree::
   :maxdepth: 1
   
　前項では，最小二乗法がうまく機能して，測定データによく沿う傾き，切片を決めることができた．
ところでpythonでは様々な機能がモジュールとして提供されているはずであった．最小二乗法もありそうである．

.. index:: optimize

optimize
========

検索などをして調べてみると，scipyのoptimizeにleastsqという機能があるのが見つかる．
　`>> SciPy Optimize <https://docs.scipy.org/doc/scipy/reference/optimize.html>`_ 

optimize.leastsqを使うとコードは以下のようになる．

.. code-block:: python
    :linenos:
    
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy import optimize
    
    def f(x):
        return 5.*x + 10.
    
    #Pusedo measured data by random number
    xdata = np.linspace(-10, 10, num=101)
    np.random.seed(1234)
    ydata = f(xdata) + 5.*np.random.randn(xdata.size)

    #Least squares method with scipy.optimize 
    def fit_func(parameter,x,y):
        a = parameter[0]
        b = parameter[1]
        residual = y-(a*x+b)
        return residual
    
    parameter0 = [0.,0.]
    result = optimize.leastsq(fit_func,parameter0,args=(xdata,ydata))
    print(result)
    a_fit=result[0][0]
    b_fit=result[0][1]

    print(a_fit,b_fit)

    #PLot
    plt.figure(figsize=(8,5))
    plt.plot(xdata,ydata,'bo', label='Exp.') 
    plt.plot(xdata,a_fit*xdata+b_fit,'k-', label='fitted line', linewidth=10, alpha=0.3) 
    plt.plot(xdata,f(xdata),'r-', label='mother line', linewidth=3) 
    plt.xlabel('Time')
    plt.ylabel('Verocity')
    plt.legend(loc='best',fancybox=True, shadow=True)
    plt.grid(True)
    plt.show()

１４～１８行目：近似する関数形を記述する．a,bが決めたい係数である．

２１行目：主要部分．a,bの初期値parameter0から初めて，(xdata,ydata)に対してfit_funcが最小になるような，a,bを返してくれる．

２３，２４行目：得られた結果をa_fit,b_fitに格納．

得られたa,bの値は，

::

    (4.9868427708984067, 10.188239312198032)

やはり測定データによく一致した直線が得られている．また当然だが，測定データを生成するのにつかった　:math:`f(x) = 5x + 10`　にもぴたりと一致している．

.. figure:: ../img/212_fitted_line2.png
   :scale: 70 %
   :alt: fitted line by optimize

optimize.leastsqは，１次関数以外にも使える．「データの描画」項で読み込んだデータはいかにも２次関数で近似できそうである．
　:math:`f(x) = ax^2 + bx + c`　

.. code-block:: python
    :linenos:

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

    file_id = 'sample_data.csv'
    #file_path = '/Users/kentaro/work_space/python/'
    file_path = './'

    rfile = file_path + file_id
    data = np.loadtxt(rfile, comments='#' ,delimiter=',')

    x_csv = data[:,0]
    y_csv = data[:,1]

    #Least squares method with scipy.optimize 
    def fit_func(parameter,x,y):
        a = parameter[0]
        b = parameter[1]
        c = parameter[2]
        residual = y-(a*x**2+b*x+c)
        return residual
    
    parameter0 = [0.,0.,0.]
    result = optimize.leastsq(fit_func,parameter0,args=(x_csv,y_csv))
    a_fit=result[0][0]
    b_fit=result[0][1]
    c_fit=result[0][2]
    print(a_fit,b_fit,c_fit)
    
    #PLot
    plt.figure(figsize=(8,5))
    plt.plot(x_csv,y_csv,'bo', label='Exp.') 
    plt.plot(x_csv,a_fit*x_csv**2+b_fit*x_csv+c_fit,'k-', label='fitted parabora', linewidth=10, alpha=0.3) 
    plt.xlabel('X-axis')
    plt.ylabel('Y-axis')
    plt.legend(loc='best',fancybox=True, shadow=True)
    plt.grid(True)
    plt.show()

::

    (0.88720506883977268, 0.41856724361772535, 5.7365571795694734)

もともと測定データを作るのにつかった一次関数の傾きa=5,切片b=10であった．これにとても近い数値が得られている．
グラフを描けば以下の通り．最小二乗法により得られた１次関数がデータによく沿っている．

.. figure:: ../img/212_fitted_parabora.png
   :scale: 70 %
   :alt: fitted parabora