カーブフィッティング手法 scipy.optimize.curve_fit の使い方を理解する.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pythonを使ってカーブフィッティング（曲線近似）する方法として、 scipy.optimize.curve_fit を使う方法がありますが、使い方が少し理解しにくいと思ったので整理してみました。\n",
    "\n",
    "# 用いる実験値\n",
    "\n",
    "<a href=\"https://qiita.com/maskot1977/items/f9e84fb6e62f7015d6c4\">Numpy.polyfit を使ったカーブフィッティング</a> で用いたのと同じデータを用いて、比較してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_observed = [9, 28, 38, 58, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198, 208, 218, 228, 238, 278, 288, 298]\n",
    "y_observed = [51, 80, 112, 294, 286, 110, 59, 70, 56, 70, 104, 59, 59, 72, 87, 99, 64, 60, 74, 151, 157, 57, 83]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python の List 型から、Numpy の Array 型に変換しておきましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "x_observed = np.array(x_observed)\n",
    "y_observed = np.array(y_observed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# １次式近似\n",
    "\n",
    "まず、１次式に近似してみましょう。関数 func1 を定義します。a と b が、求めたい値になります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func1(X, a, b): # １次式近似\n",
    "    Y = a + b * X\n",
    "    return Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "関数 func1 の使用例はこちらです。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 20,  58,  78, 118, 178, 198, 218, 238, 258, 278, 298, 318, 338,\n",
       "       358, 378, 398, 418, 438, 458, 478, 558, 578, 598])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func1(x_observed, 2, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さて、scipy.optimize.curve_fit を使って近似してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([125.77023172,  -0.1605313 ])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(func1,x_observed,y_observed) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで得られた popt が最適推定値を格納しています。<a href=\"https://qiita.com/maskot1977/items/f9e84fb6e62f7015d6c4\">Numpy.polyfit を使ったカーブフィッティング</a> で得られた推定値と比較してみましょう。\n",
    "\n",
    "# 2次式近似\n",
    "\n",
    "次に、２次式に近似してみましょう。関数 func2 を定義します。a と b と c が、求めたい値になります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func2(X, a, b, c): # 2次式近似\n",
    "    Y = a + b * X + c * X ** 2\n",
    "    return Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "関数 func2 の使用例はこちらです。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   182,   1626,   2966,   6846,  15666,  19406,  23546,  28086,\n",
       "        33026,  38366,  44106,  50246,  56786,  63726,  71066,  78806,\n",
       "        86946,  95486, 104426, 113766, 155126, 166466, 178206])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func2(x_observed, 2, 2, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さて、scipy.optimize.curve_fit を使って近似してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.39787684e+02, -4.07809516e-01,  7.97183226e-04])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(func2,x_observed,y_observed) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで得られた popt が最適推定値を格納しています。<a href=\"https://qiita.com/maskot1977/items/f9e84fb6e62f7015d6c4\">Numpy.polyfit を使ったカーブフィッティング</a> で得られた推定値と比較してみましょう。\n",
    "\n",
    "次のようにすれば、最適推定値を用いた近似曲線が描けます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([136.1819702 , 128.9940092 , 125.44205497, 118.81645644,\n",
       "       110.07383349, 107.47849913, 105.04260142, 102.76614035,\n",
       "       100.64911593,  98.69152815,  96.89337701,  95.25466253,\n",
       "        93.77538468,  92.45554348,  91.29513893,  90.29417102,\n",
       "        89.45263976,  88.77054514,  88.24788717,  87.88466585,\n",
       "        88.02614699,  88.46010889,  89.05350743])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func2(x_observed, 1.39787684e+02, -4.07809516e-01,  7.97183226e-04)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3次式近似\n",
    "\n",
    "同様に、３次式に近似してみましょう。関数 func3 を定義します。a と b と c と d が、求めたい値になります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func3(X, a, b, c, d): # 3次式近似\n",
    "    Y = a + b * X + c * X ** 2 + d * X ** 3\n",
    "    return Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "さて、scipy.optimize.curve_fit を使って近似してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(func3,x_observed,y_observed) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで得られた popt が最適推定値を格納しています。<a href=\"https://qiita.com/maskot1977/items/f9e84fb6e62f7015d6c4\">Numpy.polyfit を使ったカーブフィッティング</a> で得られた推定値と比較してみましょう。\n",
    "\n",
    "次のようにすれば、最適推定値を用いた近似曲線が描けます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 93.97825188, 118.32181643, 126.93087413, 136.59795028,\n",
       "       135.72770915, 132.27386543, 127.62777811, 122.02323006,\n",
       "       115.69400413, 108.87388316, 101.79665001,  94.69608754,\n",
       "        87.80597859,  81.36010602,  75.59225267,  70.73620141,\n",
       "        67.02573508,  64.69463654,  63.97668864,  65.10567422,\n",
       "        92.76660851, 106.63700433, 123.75703075])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func3(x_observed, 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 何次式になっても使える汎用的な関数を作りたい\n",
    "\n",
    "今まで１次式、２次式、３次式と別々の関数を作りましたが、これだと面倒すぎるので、何次式になっても使える汎用的な関数を作りましょう。何次式に近似するかはパラメータの数で自動的に判定するようにします。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "def func(X, *params):\n",
    "    Y = np.zeros_like(X)\n",
    "    for i, param in enumerate(params):\n",
    "        Y = Y + np.array(param * X ** i)\n",
    "    return Y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次の２つのコードを実行して、上記の func2 による計算結果と、func3 による計算結果と同じであることを確認しましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([136.1819702 , 128.9940092 , 125.44205497, 118.81645644,\n",
       "       110.07383349, 107.47849913, 105.04260142, 102.76614035,\n",
       "       100.64911593,  98.69152815,  96.89337701,  95.25466253,\n",
       "        93.77538468,  92.45554348,  91.29513893,  90.29417102,\n",
       "        89.45263976,  88.77054514,  88.24788717,  87.88466585,\n",
       "        88.02614699,  88.46010889,  89.05350743])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func(x_observed, 1.39787684e+02, -4.07809516e-01,  7.97183226e-04)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 93.97825188, 118.32181643, 126.93087413, 136.59795028,\n",
       "       135.72770915, 132.27386543, 127.62777811, 122.02323006,\n",
       "       115.69400413, 108.87388316, 101.79665001,  94.69608754,\n",
       "        87.80597859,  81.36010602,  75.59225267,  70.73620141,\n",
       "        67.02573508,  64.69463654,  63.97668864,  65.10567422,\n",
       "        92.76660851, 106.63700433, 123.75703075])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "func(x_observed, 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 何次式になっても使える汎用的な関数による多項式近似\n",
    "\n",
    "ここまでくれば、何次式にでも多項式近似できます。ただし、パラメータの数を指定するために p0= で初期値を必要数だけ入力する必要が生じます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([125.77023172,  -0.1605313 ])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1]) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.39787684e+02, -4.07809516e-01,  7.97183226e-04])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1]) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1, 1]) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-7.54495613e+01,  1.13179580e+01, -1.50591743e-01,  7.02970546e-04,\n",
       "       -1.07313869e-06])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit  \n",
    "popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1, 1, 1]) # poptは最適推定値、pcovは共分散\n",
    "popt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上記の計算結果を、<a href=\"https://qiita.com/maskot1977/items/f9e84fb6e62f7015d6c4\">Numpy.polyfit を使ったカーブフィッティング</a> で得られた推定値と比較してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}