{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  非線形最小2乗法の原理
2  具体的な手順
3  Mapleによる解法の指針
4  python code
5  Gauss-Newton法に関するメモ
6  課題
6.1  一山ピークへのフィット
6.2  二山ピークのフィット
7  解答例
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "
\n", "非線形最小2乗法(NonLinearFit)\n", "
\n", "
\n", "
\n", "file:/Users/bob/Github/TeamNishitani/jupyter_num_calc/nonlinearfit\n", "
\n", "https://github.com/daddygongon/jupyter_num_calc/tree/master/notebooks_python\n", "
\n", "cc by Shigeto R. Nishitani 2017 \n", "
\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 非線形最小2乗法の原理\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "前章では，データに近似的にフィットする最小二乗法を紹介した．ここでは，フィット式が多項式のような線形関係にない関数の最小二乗法を紹介する．図のようなデータにフィットする場合を考えよう．\n", "\n", "![C9_NonLinearFitplot2d1.png](figs/C9_NonLinearFitplot2d1.png)\n", "\n", "このデータにあてはめるのはローレンツ関数，\n", "\n", "$$\n", "F \\left(x;\\mathbf{a} \\right)=a _{1}+ \\frac{a _{2}}{a _{3}+\\left(x -a _{4}\\right)^{2}}\n", "$$\n", "である．この関数の特徴は，今まで見てきた関数と違いパラメータが線形関係になっていない．誤差関数は，いままでと同様に\n", "\n", "$$\n", "\\chi ^{2}\\left(\\mathbf{a} \\right)={\\sum_i^N }d _{i }^{2}=\\sum_i^N \\left(F \\left(x _{i };a \\right)-y _{i }\\right)^{2}\n", "$$\n", "で，$a={a_0, a_1,..}$をパラメータとして変えた時に最小となる値を求める点もかわらない．しかし，線形の最小二乗法のように微分しても一元の方程式にならず，連立方程式を単に解くだけでは求まらない．\n", "\n", "そこで図のような2次関数の最小値を求める場合を考える．最小値の点$a_0$のまわりで，Taylor展開すると，$\\mathbf{d,D}$をそれぞれの係数とすると，\n", "\n", "$$\n", "\\chi^2 \\left( \\mathbf{a} \\right)= \\chi^2 \\left( \\mathbf{a_0} \\right) - \\mathbf{d} \\left(\\mathbf{a}-\\mathbf{a_0} \\right) +\\frac{1}{2} \\mathbf{D} \\left(\\mathbf{a}-\\mathbf{a_0} \\right)^{2}\n", "$$\n", "である．最小の点$a_0$は，微分が$0$になるので，\n", "\n", "$$\n", "\\mathbf{a _{0}}=\\mathbf{a} + \\mathbf{D} ^{-1} \\times (-\\mathbf{d})\n", "$$\n", "と予測される．図を参照して上の式を導け．またその意味を考察せよ．\n", "\n", "![C9_NonLinearFitplot2d2.png](figs/C9_NonLinearFitplot2d2.png)\n", "![C9_NonLinearFitplot2d3.png](figs/C9_NonLinearFitplot2d3.png)\n", "\n", "\n", "現実には高次項の影響で計算通りにはいかず，単に最小値の近似値を求めるだけである．これは，$ \\chi \\left(\\mathbf{a} \\right) ^{2}$の微分関数の解をNewton法で求める操作に対応する．つまり，この操作を何度も繰り返せばいずれ解がある精度で求まるはず．\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 具体的な手順\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "パラメータの初期値を\n", "\n", "$$\n", "a_{{0}}+\\Delta\\,a,\\,b_{{0}}+\\Delta\\,b,\\,c_{{0}}+\\Delta\\,c,\\,d_{{0}}+\\Delta\\,d\n", "$$\n", "とする．このとき関数$f$を真値$a_0, b_0, c_0, d_0$のまわりでテイラー展開し，高次項を無視すると\n", "\n", "$$\n", "\\Delta\\,f=f \\left( a_{{0}}+\\Delta\\,a_{{1}},b_{{0}}+\\Delta\\,b_{{1}},c_{{0}}+\\Delta\\,c_{{1}},d_{{0}}+\\Delta\\,d_{{1}} \\right) -f \\left( a_{{0}},b_{{0}},c_{{0}},d_{{0}} \\right)\n", "$$\n", "\n", "\n", "$$\n", "=\\left(\\frac{\\partial }{\\partial a }f \\right)_{0}\\Delta a _{1}+\\left(\\frac{\\partial }{\\partial b }f \\right)_{0}\\Delta b _{1}+\\left(\\frac{\\partial }{\\partial c }f \\right)_{0}\\Delta c _{1}+\\left(\\frac{\\partial }{\\partial d }f \\right)_{0}\\Delta d _{1}\n", "$$\n", "となる．\n", "\n", "課題でつくったデータはt = 1からt = 256までの時刻に対応したデータ点$f_{1},\\,f_{2},\\,\\cdots f_{256}$とする．各測定値とモデル関数から予想される値との差$\\Delta f_1,\\Delta f_2,\\cdots,\\Delta f_{256}$は，\n", "$$\n", "\\left(\\begin{array}{c}\\Delta f _{1} \\\\\\Delta f _{2} \\\\ \\vdots \\\\\\Delta f _{256} \\\\\\end{array}\\right)=J \\left(\\begin{array}{c}\\Delta a _{1} \\\\\\Delta b _{1} \\\\\\Delta c _{1} \\\\\\Delta d _{1} \\\\\\end{array}\\right)\n", "$$\n", "となる．ここで$J$はヤコビ行列と呼ばれる行列で，4列256行\n", "$$\n", "J =\\left(\\begin{array}{cccc}\\left(\\frac{\\partial }{\\partial a }f \\right)_{1} & \\left(\\frac{\\partial }{\\partial b }f \\right)_{1} & \\left(\\frac{\\partial }{\\partial c }f \\right)_{1} & \\left(\\frac{\\partial }{\\partial d }f \\right)_{1} \\\\ \\vdots & \\vdots & \\vdots & \\vdots \\\\\\left(\\frac{\\partial }{\\partial a }f \\right)_{256} & \\left(\\frac{\\partial }{\\partial b }f \\right)_{256} & \\left(\\frac{\\partial }{\\partial c }f \\right)_{256} & \\left(\\frac{\\partial }{\\partial d }f \\right)_{256} \\\\\\end{array}\\right)\n", "$$\n", "である．このような矩形行列の逆行列は転置行列$J^T$を用いて，`\n", "$$\n", "J ^{-1}=\\left(J ^{T }J \\right)^{-1}J ^{T }\n", "$$\n", "と表わされる．したがって，真値からのずれは\n", "$$\n", "\\left(\\begin{array}{c}\\Delta a_2 \\\\\\Delta b_2 \\\\\\Delta c_2 \\\\\\Delta d_2 \\\\\\end{array}\\right)\n", "=\\left(J ^{T }J \\right)^{-1}J ^{T }\n", "\\left(\\begin{array}{c}\\Delta f _{1} \\\\\\Delta f _{2} \\\\ \\vdots \\\\\\Delta f _{256} \\\\\\end{array}\\right)\n", "$$\n", "で求められる．理想的には$(\\Delta a_2,\\,\\Delta b_2,\\,\\Delta c_2,\\,\\Delta d_2)$は$(\\Delta a,\\,\\Delta b,\\,\\Delta c,\\,\\Delta d)$に一致するはずだが，測定誤差と高次項のために一致しない．初期値に比べ，より真値に近づくだけ．そこで，新たに得られたパラメータの組を新たな初期値に用いて，より良いパラメータに近付けていくという操作を繰り返す．新たに得られたパラメータと前のパラメータとの差がある誤差以下になったところで計算を打ち切り，フィッティングの終了となる．\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Mapleによる解法の指針\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "線形代数計算のためにサブパッケージとしてLinearAlgebraを呼びだしておく．\n", "```maple\n", "> restart; \n", " with(plots): \n", " with(LinearAlgebra):\n", "```\n", "\n", "\n", "データを読み込む．\n", "```maple\n", "> ndata:=8: \n", " f1:=t->subs({a1=1,a2=10,a3=1,a4=4},a1+a2/(a3+(t-a4)^2) );\n", "```\n", "$$\n", "{\\it f1}\\, := \\,t\\mapsto 1+10\\, \\left( 1+ \\left( t-4 \\right) ^{2} \\right) ^{-1}\n", "$$\n", "データの表示\n", "```maple\n", "> T:=[seq(f1(i),i=1..ndata)]:\n", " listplot(T); \n", " l1:=listplot(T):\n", "```\n", "\n", "![C9_NonLinearFitplot2d4.png](figs/C9_NonLinearFitplot2d4.png)\n", "\n", "\n", "\n", "ローレンツ型の関数を仮定し，関数として定義．\n", "```maple\n", "> f:=t->a1+a2/(a3+(t-a4)^2); nparam:=4:\n", "```\n", "$$\n", "$$\n", "ヤコビアンの中の微分を新たな関数として定義．\n", "```maple\n", "> for i from 1 to nparam do \n", " dfda||i:=unapply(diff(f(x),a||i),x); \n", " end do;\n", "```\n", "$$\n", "{\\it dfda1}\\, := \\,x\\mapsto 1 \\notag \\\\\n", "{\\it dfda2}\\, := \\,x\\mapsto \\left( {\\it a3}+ \\left( x-{\\it a4} \\right) ^{2} \\right) ^{-1}\n", "\\notag \\\\\n", "{\\it dfda3}\\, := \\,x\\mapsto -{\\frac {{\\it a2}}{ \\left( {\\it a3}+ \\left( x-{\\it a4} \\right) ^{2} \\right) ^{2}}} \\notag \\\\\n", "{\\it dfda4}\\, := \\,x\\mapsto -{\\frac {{\\it a2}\\, \\left( -2\\,x+2\\,{\\it a4} \\right) }{ \\left( {\\it a3}+ \\left( x-{\\it a4} \\right) ^{2} \\right) ^{2}}} \\notag\n", "$$\n", "ここで，\"$||$\"は連結作用素とよばれるMapleのコマンドで，$dfda||1 \\mapsto dfda1$と連結する．\n", "初期値を仮定して，データとともに関数を表示．\n", "```maple\n", "> g1:=Vector([1,8,1,4.5]): \n", " guess1:={}: \n", " for i from 1 to nparam do\n", " guess1:={op(guess1),a||i=g1[i]}; \n", " end do: \n", " guess1;\n", "```\n", "$$\n", "\\left\\{ {\\it a1}=1,{\\it a2}=8,{\\it a3}=1,{\\it a4}= 4.5\\right\\}\n", "$$\n", "\n", "```maple\n", "> p1:=plot(subs(guess1,f(x)),x=1..ndata): \n", " display(l1,p1); \n", "```\n", "\n", "![C9_NonLinearFitplot2d5.png](figs/C9_NonLinearFitplot2d5.png)\n", "\n", "見やすいように，小数点以下を3桁表示に制限する．\n", "```maple\n", "> interface(displayprecision=3):\n", "> df:=Vector([seq(subs(guess1,T[i]-f(i)),i=1..ndata)]);\n", "```\n", "$$\n", "{\\it df}\\, := \\, \\left[ \\begin {array}{c} 0.396\\\\ 0.897\\\\ 2.538\\\\ 3.600\\\\ - 1.400\\\\ - 0.462\\\\ - 0.103\\\\ - 0.016\\end {array} \\right]\n", "$$\n", "\n", "```maple\n", "> Jac:=Matrix(ndata,nparam): \n", " for i from 1 to ndata do \n", " for j from 1 to nparam do\n", " Jac[i,j]:=evalf(subs(guess1,dfda||j(i))); \n", " end do: \n", " end do:\n", " Jac;\n", "```\n", "$$\n", "\\left[ \\begin {array}{cccc} 1.0& 0.075&- 0.046&- 0.319\\\\ 1.0& 0.138&- 0.152&- 0.761\\\\ 1.0& 0.308&- 0.757&- 2.272\\\\ 1.0& 0.800&- 5.120&- 5.120\\\\ 1.0& 0.800&- 5.120& 5.120\\\\ 1.0& 0.308&- 0.757& 2.272\\\\ 1.0& 0.138&- 0.152& 0.761\\\\ 1.0& 0.075&- 0.046& 0.319\\end {array} \\right]\n", "$$\n", "\n", "```maple\n", "> tJac:=(MatrixInverse(Transpose(Jac).Jac)).Transpose(Jac);\n", "```\n", "$$\n", "{\\it tJac}\\, := \\, \\left[ \\begin {array}{cccccccc} 0.565& 0.249&- 0.354& 0.040& 0.040&- 0.354& 0.249& 0.565\\\\ - 2.954&- 0.506& 4.012&- 0.552&- 0.552& 4.012&- 0.506&- 2.954\\\\ - 0.352&- 0.029& 0.557&- 0.176&- 0.176& 0.557&- 0.029&- 0.352\\\\ - 0.005&- 0.012&- 0.035&- 0.080& 0.080& 0.035& 0.012& 0.005\\end {array} \\right] \n", "$$\n", "\n", "```maple\n", "> g2:=tJac.df; \n", " g1:=g1+g2;\n", "```\n", "$$\n", "{\\it g2}\\, := \\, \\left[ \\begin {array}{c} - 0.235\\\\ 5.592\\\\ 0.613\\\\ - 0.520\\end {array} \\right] \\notag \\\\\n", "{\\it g1}\\, := \\, \\left[ \\begin {array}{c} 0.765\\\\ 13.592\\\\ 1.613\\\\ 3.980\\end {array} \\right] \\notag\n", "$$\n", "これをまたもとの近似値(guess)に入れ直して表示させると以下のようになる．カーブがデータに近づいているのが確認できよう．この操作をずれが十分小さくなるまで繰\n", "り返す．\n", "```maple\n", "> guess1:={seq(a||i=g1[i],i=1..nparam)};\n", " p1:=plot(subs(guess1,f(x)),x=1..ndata):\n", " display(l1,p1);\n", "```\n", "$$\n", "guess1:=\\{a1=0.765, a2=13.592, a3=1.613, a4=3.980\\}\n", "$$\n", "\n", "![C9_NonLinearFitplot2d6.png](figs/C9_NonLinearFitplot2d6.png)\n", "\n", "\n", "4回ほど繰り返すと以下の通り，いい値に収束している．\n", "$$\n", "guess1:=\\{a1 = 1.006, a2 = 9.926, a3 = .989, a4 = 4.000\\}\n", "$$\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# python code\n", "\n", "幾つもの関数が用意されている．\n", "* curve_fit\n", "* curve_fit with bounds\n", "* least square fit\n", "\n", "全部を理解することはできないが，manualを見ながら求めることができるといいね．\n", "boundsとかparamsの初期値が重要．" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lFXe//H3SSUJoQSQXhVBBJQmICA1CARBxQXsbe2/\nLa5lLavuPq7l2XUtuz67yto7LMKqFDFMCL0jgnQFQZRAIDQhIWXO74+T0KVlJveUz+u65koymcz9\nHUg+Ofne55zbWGsREZHwF+N1ASIiEhgKdBGRCKFAFxGJEAp0EZEIoUAXEYkQCnQRkQihQBcRiRAK\ndBGRCKFAFxGJEHEVebCaNWvaJk2aVOQhRUTC3uLFi7dba2ud7HEVGuhNmjRh0aJFFXlIEZGwZ4zZ\neCqPU8tFRCRCKNBFRCKEAl1EJEIo0EVEIoQCXUQkQijQRUQihAJdRCRCKNAl9G3cCGPGeF2FSMhT\noEtomzMHOnWCESNgxQqvqxEJaQp0CV0ffgh9+kBysvs4M9PbekRCnAJdQo+18OSTcM01cNFFsHgx\nNG+uQBc5CQW6hJYDB+CGG+Dxx93bzEyoUQPS02H6dCgs9LpCkZClQJfQsX079OsH770Hf/4zvPUW\nJCa6z6Wnw759MHeupyWKhLIK3W1R5GetWQMZGbB5s+udjxzJhg0wbpy7/fRjb5bGxmIyM6FnT6+r\nFQlJGqGL97KzoWtX2LOH796cxp+/GUn79tCsGdx/vxu4L/uuKtubXaQ+usgJKNDFU/aNN7H90tkW\nW4d+qfNpek1XHnsMKlWC556Db7+FVaugbl2YFpsOixbBzp1ely0SkhToUuH8fpg908/Uzo9gbr2F\nqSW9OC9vDrZJU15+GX74wU0/v+8+N0qPiYErr4RX16e7L542zeuXIBKSThroxpg3jDHbjDFfH3Zf\nmjEm0xizrvRt9eCWKeGuqAh8Prj7bji7Xj4/XjKCfgue4fOGt7H51Ums2VoNnw/uuQfq1Tv264cN\ngxmFnSmqVFltF5GfcSoj9LeAAUfd9xDgs9Y2B3ylH4sc19y5rmXSrx9MfmsrXxT24irzMQV/fo4B\nG1/l5tvjqVnzxM/RowdUqxnP8rReCnSRn3HSQLfWzgDyjrp7KPB26ftvA5cHuC6JIK++6jolU1/8\nmm9rdab5ga8x48ZR6dH7wJhTeo64OLj8cvhoR7prrG/YEOSqRcLPmfbQa1trt5S+nwPU/rkHGmNu\nN8YsMsYsys3NPcPDSbiy1rVa7ms9hb6PXUxMUSHMmOHS+TQNGwafHUh3H2iULnKMcp8UtdZawJ7g\n86OstR2ttR1r1apV3sNJmFm3DlpuzuTh2RnuDOf8+dChwxk9V58+sKVKS/KS6yvQRY7jTAN9qzGm\nLkDp222BK0kiSVYWXMv72CrVYOZMaNjwjJ8rIQGGDDV8XpyO9fmgpCSAlYqEvzMN9E+BG0vfvxH4\nJDDlSKTx+aBP7HRi+vSE1NRyP9+wYTChMB2zcycsWRKACkUix6lMW/wQmAu0MMZsNsbcCjwLpBtj\n1gH9Sj8WOYLfD+umbqRRyXeYAC3X798f5ib3cx+o7SJyhFOZ5XK1tbautTbeWtvAWvu6tXaHtbav\ntba5tbaftfboWTAiLFsGbXdNdx/06hWQ50xKgosGn8XXcRdgM6cG5DlFIoVWikrQ+HzQi2z81dKg\ndeuAPe+wYTC5OB07ezbs3x+w5xUJdwp0CRqfD9Ljs4npdYlbvx8ggwbB9Pj0Q1MgRQRQoEuQFBXB\n+uxNNCzaELB2S5nKlSG5f3cOkID9Qn10kTIKdAmKBQugU35p/zwI+5dfNiKZWXRn/6cKdJEyCnQJ\nCtc/n46/WnVo2zbgz3/ZZZAVk07Kt8shJyfgzy8SjhToEhRZWdA/MZuYS3oEtH9eplo12NvFbQOg\n2S4ijgJdAm7/ftg0ZzMND3wb8P754S64qR3bqUHeGLVdRECBLkEwezZ0LQrs/PPjGXpFDFmmL/HZ\nmW4XMJEop0CXgPP5oI/JxlatGpT+eZmaNeH7FulU+WkLrFwZtOOIhAsFugSczwf9K03HXHIJxMYG\n9Vi1rnF99Jz31EcXUaBLQO3cCVsW/0jD/HVBma54tPRfNmYtzdkzTn10EQW6BNT06dDDBr9/XqZu\nXVhRpx8NvsmGwsKgH08klCnQJaB8PugXl42tUgUuvLBCjhk7IJ1k/z42j51XIccTCVUKdAmorCy4\nNCEb06NH0PvnZS68tzclxLDpdbVdJLop0CVgtmyBvJVbaLB/bYW0W8o0aluNFSkXkTpfgS7RTYEu\nAZOVBT0J3v4tJ7K3czqt9i3k++U7K/S4IqFEgS4Bk5UF/ROmY1NToV27Cj12o1vTicXP0hemVehx\nRUKJAl0CwtrS+eeJ2Zju3SEurkKP3/AXXdgXU5mSyWq7SPRSoEtArF8P+Ru30mDv6grtnx8UH8/m\nc3rROidTmy9K1FKgS0Ac0T/3ItCB1CvTOYdvmfrvDZ4cX8RrCnQJCJ8PMpKzsZUrQ/v2ntRQ93q3\nDcC2D7QNgEQnBbqUm9/vRuh946d70j8vY85ryZ7K9Wi4JpMdOzwpQcRTCnQptxUrgNxtNNi9ssKn\nKx7BGIp6pdPH+vh0fIl3dYh4RIEu5ebzwSXMcB941D8vkzYinRrksfTNLz2tQ8QLCnQpt6wsGFo1\nG1JSoEMHT2sx6f0ASJ2Xya5dnpYiUuEU6FIuxcVuh8U+sdOhWzeIj/e2oNq12XdOW/r4M5kwwdtS\nRCqaAl3KZdEiSNiTS/28rz1vt5RJHpJON2bz2ej9XpciUqEU6FIuWVmH9c+9PCF6GNM/nUQKyZ8y\ng59+8roakYqjQJdy8flgWM3pkJwMHTt6XY7Towf++AR6FmUyebLXxYhUHAW6nLH8fJg9G3qRDRdf\nDAkJXpfkJCdjundnQOxUPv7Y62JEKo4CXc7Y3LmQcmAH9bYvD5n+eRmT3o/zS5ax4LOtFBR4XY1I\nxShXoBtj7jXGrDDGfG2M+dAYUylQhUno8/mgd0xozD8/RrrbBqDr/ql88YXHtYhUkDMOdGNMfeDX\nQEdrbWsgFhgZqMIk9Pl8MLx2NiQlQadOXpdzpHbtsGlpZMRnqu0iUaO8LZc4IMkYEwckAz+WvyQJ\nB7t3w8KF0MM/PbT652ViYzF9+zIwLpNPP7EUFnpdkEjwnXGgW2t/AJ4DNgFbgN3WWv1xGyVmzICq\n/jzqbFsWMtMVj5GeTvX8H6m7exXTdCEjiQLlablUB4YCTYF6QIox5rrjPO52Y8wiY8yi3NzcM69U\nQorPB33jZ2KsDb3+eZnSPvrgBLVdJDqUp+XSD9hgrc211hYB44CLj36QtXaUtbajtbZjrVq1ynE4\nCSVZWTCyTjZUqgQXXeR1OcfXpAmccw4ja2YyfrzbpkAkkpUn0DcBXYwxycYYA/QFVgWmLAll27bB\n8uXQrWQ6dO0KiYlel/Tz0tNpsyOb3dsLWbLE62JEgqs8PfT5wFhgCbC89LlGBaguCWFZWVCNndTe\nsjR02y1l0tOJP7CPzszH5/O6GJHgKtcsF2vtE9baltba1tba6621BwJVmISurCwYkFzaPw/VE6Jl\neveGmBiuPSuTrCyvixEJLq0UldPm88HIutNdq6VzZ6/LObFq1aBTJwbEZDJrFhzQkEMimAJdTst3\n38H69dD1QDZ06eJOioa69HQab1tApYKdzJvndTEiwaNAl9OSlQVV2UWtH8Ogf15mwACM388A84Xa\nLhLRFOhyWnw+uKzaLIzfHz6B3qULpKVxfdpEnRiViKZAl1NmrRuhj6iT7Zb6h3r/vExsLAwYwCX7\nJ7NwXokueiERS4Eup2zVKsjJgc4F092oNynJ65JOXUYGlfO3065kIbNmeV2MSHAo0OWU+XxQhd3U\n3LQk9KcrHm3AAGxMDJfFTFIfXSKWAl1OWVYWDKs9O7z652XS0jBdu3JV8kQFukQsBbqckpISyM6G\n4Wdlu/55ly5el3T6MjJo8dMSfly8hbw8r4sRCTwFupySJUtg1y7otH+624wrOdnrkk5fRgYAA5nE\n9Oke1yISBAp0OSVZWZDKHtK+Wxx+7ZYybdpgGzRgSKzaLhKZFOhySnw+uLbxbExJSfidEC1jDCYj\ng3QymTlVewBI5FGgy0kdOACzZsGwmtMhPt5tmRuuMjJILvmJmqtnkpPjdTEigaVAl5OaNw/y86H9\n3mx3MeiUFK9LOnN9+uBPSCQDtV0k8ijQ5aR8Pkg1P1H920Xh2z8vk5KC6d2by2IU6BJ5FOhyUpmZ\ncEuL0v55uAc6YAZncI5/HeunrPO6FJGAUqDLCe3eDQsWwJVp2RAXBxcfc9nY8DNoEABtN09kwwaP\naxEJIAW6nNC0aeD3wwW7p4d//7xMs2YcaNqSQUxi2jSvixEJHAW6nFBmJtRK3keVNQvDd7ricSRc\nkUFPpjPrc229KJFDgS4nNHUq3Hb+HExxcUT0z8uYwRkkUojNnIq1XlcjEhgKdPlZmzbB2rVwZcoU\nN/88EvrnZbp3pzCpChfvmsjq1V4XIxIYCnT5WVOnAljarB4Dl14KqalelxQ48fEU9urPICaR5dMQ\nXSKDAl1+VmYmDK4xj4Sc72H4cK/LCbiUX2RQnx/ZMH6p16WIBIQCXY7L73cLiu6pNcZtlztkiNcl\nBZwZNBCA6nMmUlLicTEiAaBAl+Natgy25/rpsfU/MHAgVK3qdUmBV7s225t1ok/BRL76yutiRMpP\ngS7HNXUqXMwcUnb+EJHtljIJV2TQmfnM/TTX61JEyk2BLseVmQl3VR8NlSrBZZd5XU7QVBmZQQyW\n/PGfe12KSLkp0OUYBQUwe0YJgw+MdcvkI2l2y9Hat2dPUm2arJhIUZHXxYiUjwJdjjFnDnQomEXV\n/TkwYoTX5QRXTAx5XQbRt2QKC+YUe12NSLko0OUYmZlwtRmNTUo6eB3OSFbjhgyqs4t178z1uhSR\nclGgyzGmZRYzPO5jzODBkbEZ10mkXplOEXEkZE70uhSRclGgyxF27IDKi6eTVrQt8tstZapUYWPD\nHrTdPJH8fK+LETlz5Qp0Y0w1Y8xYY8xqY8wqY0wYX2xSwG2X+wvGUJKU4uafR4mi/hm0tl+zePwm\nr0sROWPlHaG/BHxurW0JXACsKn9J4iXflGKuMh9jhg6B5GSvy6kwje5y5wp2vKO2i4SvMw50Y0xV\n4BLgdQBrbaG1dlegChNv7J+QRQ27g5gRkbuY6HhS2rfgh8Rm1JivQJfwVZ4RelMgF3jTGPOlMeY1\nY0zkn0GLYOvXQ4+cMRRWSoUBA7wup2IZw6a2GbTflcXuHDXSJTyVJ9DjgPbAv6y17YB9wENHP8gY\nc7sxZpExZlFurpZXh7Kszwu5knHk9x/qVohGmeSrMkgmnzWv6Lp0Ep7KE+ibgc3W2vmlH4/FBfwR\nrLWjrLUdrbUda9WqVY7DSbBt/cBHGjupcmt0tVvKtLi9J/tIpvC/artIeDrjQLfW5gDfG2NalN7V\nF1gZkKqkwpWUQNOFY9ifUBVzaX+vy/FEpWqVWFqzH01XTkTXpZNwVN5ZLr8C3jfGLAMuBJ4uf0ni\nhaXzDzCocDw5nS+HxESvy/HM7m4Z1C/aSN4sjU0k/JQr0K21S0vbKW2ttZdba3cGqjCpWOtfzaQa\nu6l+Z5QsJvoZdW4ZBMD3oyZ5XInI6dNKUQGg2pTR7I6tTvWr+npdiqfaDmrAspgLSPKpjy7hR4Eu\n7M8roPPWT1h93hXucnNRLC4OVjXNoNmWWbBLyyokvCjQhTV/n0IV9hJ7dXS3W8r4B2YQRwk7PvzC\n61JETosCXWD0aLZTg1b39Pa6kpBw3k2d2U4Ndr2vtouEFwV6tMvPp8XaT5lbbxjJVeO9riYktG0X\ny7SEAdRaPBn8fq/LETllCvQot/vDSST797F3QHQuJjqemBj48cJBVCnIxS5Y6HU5IqdMgR7ldv97\nDFs5i+a/7Ol1KSEl9RcDKCGGne+p7SLhQ4Eezfbto/aiCUxIHEb7i+K8riakdLssjbl0peQzBbqE\nDwV6FLMTJpJYvJ9NnYcTG+t1NaHl3HNhZmoGtTYtgS1bvC5H5JQo0KPYT2+OYQt1qDeih9elhBxj\n4Kee7qIX/glaNSrhQYEerfbuJSlrImO5in6Xanh+PM2vbMP3NGDvaAW6hAcFerSaMIG4ogJm1h3B\n2Wd7XUxo6t3HMJEMKs3KhMJCr8sROSkFepTyfzSaH009qmdc7HUpIatxY1hSJ4PEA3th5kyvyxE5\nKQV6NNqzByZPZoz9Bf3661vgRBIG9KGARPya7SJhQD/N0ejTT4kpKmQMI+jTx+tiQluPASlMozcH\nxivQJfQp0KPR6NFsTWxIUfvO1KjhdTGhrVcv+IzLSNq0FhZq1aiENgV6tNm1CztlCu8XDadvuv77\nT6Z2bVh83vXsiU+DP/7R63JETkg/0dHmv//FFBXxkX846eleFxMeuqSn8lf7AEyaBPPnn/wLRDyi\nQI82Y8awo0oTlid2ols3r4sJD1ddBS8U/z/2JdfUKF1CmgI9muTlQWYm/40fTo9LDJUqeV1QeOjR\nA4ZcXZmnDjwAn38Oc+d6XZLIcSnQo8n48VBczD93DKdfP6+LCS/PPw9vJt/Drvha2Cee8LockeNS\noEeTMWPYc9bZLKG9+uenqU4dePTpFJ4s+j0mMxNmzfK6JJFjKNCjRW4u+HxMP2s4NWsaLrjA64LC\nz113wbwL72JbTG2KH9UoXUKPAj1ajB8PJSW8lDOCvn3dVXnk9MTGwoujknnG/xBxM7JgxgyvSxI5\ngn6so8Xo0Rxoci6+7W3VbimHTp3Af9sd/Ehd9t6nUbqEFgV6NNi6FbKzWXruCMDohGg5/ekvSfxf\n5YdIXZSN3zfN63JEDlKgR4MPPgC/n3cKhtO8udtFUM5ctWpw/ku38wP12HrnE2Ct1yWJAAr0yLd9\nOzz5JP6evXhn8fkanQfI1TdXYmzzR6j7zUx2fpzldTkigAI98j38MOzZw9JbX+anfUb98wAxBi79\nzy/5ngZsv/txjdIlJCjQI9m8efDaa/Db3/LJN+cTEwO9e3tdVORoeUEiiy99lOa5c1j+fKbX5Yhg\nbAWOLDp27GgXLVpUYceLaiUlcNFFkJMDq1dz8aWp+P0u4yVw9u8qJK9mc7Yn1KPVzjkkJBqvS5II\nZIxZbK3teLLHaYQeqUaNgiVL4Pnn2e1PZcEC1G4JguRqCey441EuzJ/HJ3d97nU5EuXKHejGmFhj\nzJfGmAmBKEgCIDcXHnkE+vSB4cPJznYDdp0QDY4LXriJbUmNafr2E2z8Tr108U4gRui/AVYF4Hkk\nUB56CH76CV5+Gb81PPcc1KgBXbt6XViESkgg9o+P0dG/kHdGTvK6Goli5Qp0Y0wDIAN4LTDlSLnN\nnQtvvAG/+x2cdx5vvun2kfrrXyEhweviIleNe29gZ1ozBs5/gs8+1ShdvFHeEfqLwIOAPwC1SHmV\nlMDdd0P9+vDYY2zbBg88AJdcAjfd5HVxES4+ntT//QMdWcynt33G/v1eFyTR6IwD3RgzGNhmrV18\nksfdboxZZIxZlJube6aHk1PxyiuwdCm88AJUrsx997nOy6uvunnTElxxN11Pfv2zuWvbH/nzkxql\nS8Urzwi9GzDEGPMd8BHQxxjz3tEPstaOstZ2tNZ2rFWrVjkOJye0bRs8+qg783nVVUydCu+959rp\nLVt6XVyUiIsj6enHac+XrP3rJ6zSmSWpYAGZh26M6QXcb60dfKLHaR56EN10k9uzZdky8hu3pG1b\nd/fy5ehScxWpuJjiFq1YvTGJX3f/Et+0GP11JOWmeejRZNYsePttuO8+aNmSp5+Gb75xHRiFeQWL\niyPuf56gdckyqk8fz/vve12QRBOtFA13xcXQoQPs3AmrVrFyYwoXXggjR8I773hdXJQqKcG2bs23\nG+PoXvkrVq2JoXp1r4uScKYRerT45z9h2TJ48UX8SSnceSekpsLf/uZ1YVEsNhbz+OOck/81vbaP\n5dFHvS5IooUCPZzl5MBjj8Gll8IVV/DWWzBzpptzrvPPHhs+HFq14oXqf2LUv0pYuNDrgiQaKNDD\n2YMPQkEB/OMfbMs13H8/9OgBN9/sdWFCbCw88QR181ZyW/X/cOedrjsmEkwK9HA1Ywa8+65bOdS8\nOfff7+acv/KK5pyHjKuugtat+d/kP7F0SQnPPed1QRLpFOjhqKgI7rkHGjWCRx7B53PZ/vvfQ6tW\nXhcnB8XEwBNPUOWH1bxw0Uc8/ribRioSLJrlEo5efBHuvRfGj6dgwOW0aePuXrYMkpK8LU2O4vdD\nu3YU7yug6d7l1KyXwPz52ldHTo9muUSqLVvg8cdh4EAYOvTgnPN//UthHpJiYuDpp4n7di1zW97E\nV0v9PPmk10VJpFKgh5sHHoADB+Dvf2fVasOzz8J112mv85CWkQHPPkuDGR8y5bx7eeZpy4IFXhcl\nkSjO6wLkNEyfDu+/D489hj37HO7sBZUra855WHjwQcjJIf3FF3mqah1uvPFhlizRX1USWBqhh4uy\nE6FNmsBDD/HWW26iy1/+Amed5XVxclLGuN+8117L73c/wsWrX9eCIwk4BXq4+PvfYcUKeOklcvcl\nc//90L073HKL14XJKYuJcRcfufRS/m1u59sXPmX6dK+LkkiiQA8HP/wAf/wjDB4MQ4Zw//2wd6/b\n5zxG/4PhJSEBxo7FdujIaDOC/xs5k717vS5KIoXiINRt2+aCvLgYXnqJrCy36daDD2rOediqXJnY\nyRPxN2zMqzlDeP4WTU6XwFCgh7JNm9xa/jVrYNw4Cuo148474eyzUf813NWsSfKMKcRUTuaXYwcw\n7a2NXlckEUCBHqrWrHFN8pwc+OILGDiQZ56Bdes05zxiNG5MpewpVI7ZT6Nf9mfnuu1eVyRhToEe\nir780o3MCwogOxu6d2f1anjmGbj2WkhP97pACZTEDq3JGfUZ9Uo2kddlkNuQR+QMKdBDzcyZ0KuX\nu9TQrFnQrh3Wwp13QkqK5pxHoha3dueTq0fTOG8JOd2HQWGh1yVJmFKgh5LJk93e5nXrwuzZcO65\nALz8sltT9Je/QO3aHtcoQTHs7SE83WQUdb76gvyrb3Z7wIicJgV6qBg9GoYMgZYt3Yqhhg0Bd6nQ\n3/zGrR6/9VaPa5SgiY+HqybewmOxT5M07gPs7+6DCtw4TyKDAj0U/PvfcPXV0KULTJt2cOnnRx+5\nhUN9+8LYsZpzHulatYJqzz7Ei/wG89KL7k8ykdOgiPDaX/4Ct98OAwbAlClQtSoA48e7Tbe6d4f/\n/te11CXy/fZew/juz/OfuKvhoYfgzTe9LknCiALdK9bCww+7q1KMGOFSOzkZgIkT3V2dOsGECe5k\nqESH2Fh4460Ybk94i8Vp6djbboPPPvO6LAkTCnQv+P1w993w7LNudP7++weveJCZCcOGQdu27hxp\naqrHtUqFO/tseOZvCfTK+5jc+u3cBadnz/a6LAkDCvSKVlTkeimvvOJG56+84oZluJksQ4dCixau\n+1Ktmse1imfuuAMu7p9Kp9xJFNZt5M6Kf/SRTpTKCSnQK1J+PlxxBXz4oRudP/vswSs6z53rfmab\nNHGj9Bo1vC1VvGUMvP467E6oxdVpX2Cbn+tOnA8ZAps3e12ehCgFekXZs8ed+Jw06dDovNSiRe5T\ndevC1Kna31ycBg3gH/+AcYsb87er5rpVZT6fmw7zr39prrocQxeJrghbt7rh91dfwbvvwsiRBz/1\n1VfQu7eb3HLY9HMRwHVYhg1zs57i4uDcuPW8XHwHvYunsiChO3+o/W9+SG1JQoI7DZOYyHHfr1XL\njSG0MC08nepFohXowZSf7y5M8fTTbjn3xx/DoEEHP71ypVvln5jowrxpU+9KldC1axeMGuXeFhZC\n4QFLh6/fZtic35FYvI+x5z3OmCYPkl8c7z5f6C47W/Z+YaHr0qSlufPvffp4/YrkdJ1qoGOtrbBb\nhw4dbFQoKbH2vfesbdTIWrD2ssusXbXqiIesXWttnTrutmaNR3VKeNuyxdpf/MJ9j7Vta+3ChT/7\n0GXLrG3Z0lpjrH3iCWuLiyuuTCk/YJE9hYxVDz3QZsyAzp3dTJaaNSErCz791C3pL7VhgxslFRe7\nlmjpli0ip6dOHRgzxq1h2L7dfd/dfz/s23fMQ9u0gYUL4frr4U9/cjt2btniQc0SVAr0QFm71s1g\n6dnT7WH+zjvuJ6h37yMe9v33Lsz37XMnQHXVISm3oUNd/+6229yJ0zZt3EjhKJUru72B3nwT5s2D\nCy90M6okcijQy2v7dvjVr+D8811CP/WUC/frrz9m85UtW1yY5+W5a1ZccIFHNUvkqVrVzZ7KznZn\nT/v1cxsB7dx5zENvusnNrKpVy23u+Yc/uL8WJfydcaAbYxoaY6YZY1YaY1YYY34TyMJCXkGB24fl\n7LPdFLJf/hK++QYeeeS4lxPats1tsrVlC3z+OXQ8+ekNkdPXs6ebOvXww+6vxPPOczu7HTX5oVUr\nWLAAbr7ZjUH69HHXIpcwdyqN9uPdgLpA+9L3U4G1QKsTfU1EnBQtKbH2gw+sbdzYnYwaPNjaFSt+\n9uHFxdYuX+7OWSUlWZudXXGlSpT78ktr27d336dDhlg7Z461fv8xD3v3XWtTUqytWdPayZM9qFNO\nimCfFLXWbrHWLil9fy+wCqhf3l8wIW3mTLfF7TXXuDlgPp/bOKm0EV5cDMuXw1tvwa9/Dd26QZUq\nrqW5Zg188okbQIlUiAsvhPnz3V+SWVlw8cWuz/fyy24OZKnrrnMtmLp1YeBAt8ljUZGHdcsZC8g8\ndGNME2AG0Npau+eoz90O3A7QqFGjDhs3htnVza2FxYvdXPLx46F+fXj6aQqHX8fXK2NYsgSWLHEP\nWbbMdWLA7ZB44YXQvj106ACXXKJ55uKhvXvdlhOjRrlv1qQkt6Xn7be7QYox5OfDvffCq6+67P/o\nIy10CxUVtrDIGFMZmA48Za0dd6LHhtXCorVr4YMP3G3dOkqSK7O430O8nXYv85Yls3z5oVFMlSou\nuMtuHTq7qWRrAAAKgUlEQVRA8+YH99wSCS2LF7tg/+ADd1HqNm1csF93HVSrxkcfuQ/j492smMGD\nvS5YKiTQjTHxwARgirX2+ZM9PtQDfcfyH9n16mhSJ3zAWRsX4ccwL7EXbxy4ho8Zxi6qk5Z2KLTL\nArxZM11NSMLQz43a77iDdWmdGTHS8OWXcN997g/U0h2eDzpwwE2iycuDHTvc25+7FRa6RdLXXuv2\nqJHTE/RAN8YY4G0gz1r721P5mlAIdL8fNm6E1ath1SrYuHQndeaOo9t3H9C9eBoxWBbRgbHx17C0\nxQhqtK3Peee5WYnt20OjRgc3SBSJHMcZtRfdcgePrLiW516rxvnnu03jDg/p46xfOig21p1mKrsV\nFrpDGONm1NxwA1x5pZsbLydXEYHeHZgJLAfKtn17xFo76ee+pqID3e93JymnTXNTtFatcicnbX4+\ng5nANXzAICaRSCFbq5zDhq7XUnTV1TRKb0HDhhp1SxQqG7W/+qo7OZSUxHedR/Bkzm2sq9GF6jVi\njgjqw281ahx6PzX12IHPN9/Ae++5/enWr3cX6LrySrdko29ftShPJCo357IWVqxwAT5tmrtgRF6e\n+1zThsVcfZaPy/d/wAUbxpNQsJeS2nWJvWakm7XSoYOG3iKHO3rUnpbmVj736+cS+JxzzuhnxlqY\nM8cF++jRbsJNvXrux/CGG1xLX44UFYFurRtxlwV4djbk5rrPNWlsuab9agZXm0XbPbNImfm5W91T\ntarbj/Taa90cQg0LRE5s7163H9HUqe5WdoGNhg0PhXvfvm5vmdNUUOCuofvOO+5SAcXFbmblDTe4\ngD+Dp4xIERno1sK33x4Z4GUbDDWtX8hNbRYzMHUW5++cRfKXs92ZGnBrnHv3dvuQDxwIlSqV/8WI\nRCNrYd06twZj6lT3g1i2vcD557tg79fPDZaqVDmtp87NdSP2d991LdKYGOjf34X70KEHr6EelSIq\n0CdOdP/R06YdGhw0P2s3t5w3l0tTZtFy+ywqLZuPKZsE3rw5dO9+6Na8udopIsFQUgJLlx4K+Jkz\n3bA7NhY6dTo0gu/a1W38f4pWrz7Ub9+0yf1h/c477gp80SiiAv3BB2HK65u5ufks+laaRfOts0hc\nswxjrfvGadcOevRw4d2tmy7LIuKVggJ3gVyfz90WLnShn5gIrVu7fkrbtoduJ7l4rt/vfkc88IBr\n6b/8Mtx1VwW9lhASUYFedNtdxL/2ivsgJcX9ti8bfXfurLlPIqFq9243O2HGDLdp2FdfHTrRBW7l\nddu2RwZ9ixZux8jD7NvnrpH92Wdua4KnnoquWWgRFehMnuxWbnbv7v7jj/rPFpEwsnWrC/Zlyw69\nXbXq0NLrxES3P9LhQd+mDcXVa/GrXxteecWdMH3jjdPq4oS1yAp0EYlshYWucb5s2ZFBn5Nz6DGV\nK2ObNmVdSTMmrmxKXPOm3PynplRu09RtlJSS4l39QXaqga6hroh4LyHhUMvlcNu2uWD/+mtYvx6z\nYQPnbviGZgmZxK3bD9cc9tizznLBfrxbo0Zuc5oIpxG6iIQfa5k5LpcnbthAy4T1PHbdBuoWbHAX\n7N2wwU2NOfwyTDExbvpynTpun+DD3x59XxDOyblN6c+876+Wi4hEvOXL3dKSPXtg3Dg3SxJwYb55\n86GA/+47177ZsuXQ261bj3/tvZSU44d+jRqHbofvdZCcTEGBu17wpk3utnHjke9//727UlmvXmf2\nOhXoIhIVNm92OzmuWgWvv+4WIp0Sv98tPszJOTbsj367Z8/PPk0+lcgjjR3UOPh2BzUoTHGhH3dW\nGpXq16Db/V1p3v3MplSrhy4iUaFBAzdX/cor4cYb3aj40UdPYS1hWRumVq3jbiCzYQN8/LG7LZ2X\nTxp51GAHaeRRN34H56Tl0SR1B/WTdlA7Po8adgeNivNILlhN/N48zI4dsKkINgGLgDsnAwOC8C9w\niAJdRMJe1apudvOtt8Jjj7lQ/+c/T3+G89q1h0J88WJ3X/v28Mj/JHH++fVp3Lg+jRpBzZqn8AvD\nWjeBvmyz+GbNzui1nQ4FuohEhIQEtz1Ao0bughw//OC2DDnZOc6VK2HsWHdbvtzd17kz/PWvbh+/\nM750pDHu4JUrQ+PGZ/gkp0eBLiIRwxi3irRRI7j7bncScuLEI3cDsdbNhCwL8dWr3dd16wYvvuha\nN+F6LVUFuohEnDvucLsKjBjhdgqZNMlt6V4W4t9+61roPXvCr34FV1zhJrSEOwW6iESkwYPdFtsZ\nGW4nAWtdT71PH/j97+Hyy9350EiiQBeRiNWpk9v88e9/d5uyDhnipo5HKgW6iES0s8+Gl17yuoqK\nEUUbUIqIRDYFuohIhFCgi4hECAW6iEiEUKCLiEQIBbqISIRQoIuIRAgFuohIhKjQC1wYY3KBjWf4\n5TWB7QEsJxzoNUcHveboUJ7X3Nhae9KNCio00MvDGLPoVK7YEUn0mqODXnN0qIjXrJaLiEiEUKCL\niESIcAr0UV4X4AG95uig1xwdgv6aw6aHLiIiJxZOI3QRETmBsAh0Y8wAY8waY8w3xpiHvK4n2Iwx\nDY0x04wxK40xK4wxv/G6popgjIk1xnxpjJngdS0VwRhTzRgz1hiz2hizyhjT1euags0Yc2/p9/TX\nxpgPjTGVvK4p0Iwxbxhjthljvj7svjRjTKYxZl3p2+rBOHbIB7oxJhb4P2Ag0Aq42hjTytuqgq4Y\nuM9a2wroAtwTBa8Z4DfAKq+LqEAvAZ9ba1sCFxDhr90YUx/4NdDRWtsaiAVGeltVULwFDDjqvocA\nn7W2OeAr/TjgQj7QgYuAb6y16621hcBHwFCPawoqa+0Wa+2S0vf34n7Q63tbVXAZYxoAGcBrXtdS\nEYwxVYFLgNcBrLWF1tpd3lZVIeKAJGNMHJAM/OhxPQFnrZ0B5B1191Dg7dL33wYuD8axwyHQ6wPf\nH/bxZiI83A5njGkCtAPme1tJ0L0IPAj4vS6kgjQFcoE3S9tMrxljUrwuKpistT8AzwGbgC3Abmvt\nF95WVWFqW2u3lL6fA9QOxkHCIdCjljGmMvAx8Ftr7R6v6wkWY8xgYJu1drHXtVSgOKA98C9rbTtg\nH0H6MzxUlPaNh+J+mdUDUowx13lbVcWzbmphUKYXhkOg/wA0POzjBqX3RTRjTDwuzN+31o7zup4g\n6wYMMcZ8h2up9THGvOdtSUG3GdhsrS37y2ssLuAjWT9gg7U211pbBIwDLva4poqy1RhTF6D07bZg\nHCQcAn0h0NwY09QYk4A7ifKpxzUFlTHG4Hqrq6y1z3tdT7BZax+21jaw1jbB/f9mWWsjeuRmrc0B\nvjfGtCi9qy+w0sOSKsImoIsxJrn0e7wvEX4i+DCfAjeWvn8j8EkwDhIXjCcNJGttsTHm/wFTcGfF\n37DWrvC4rGDrBlwPLDfGLC297xFr7SQPa5LA+xXwfulAZT1ws8f1BJW1dr4xZiywBDeT60sicMWo\nMeZDoBdQ0xizGXgCeBYYY4y5Fbfj7PCgHFsrRUVEIkM4tFxEROQUKNBFRCKEAl1EJEIo0EVEIoQC\nXUQkQijQRUQihAJdRCRCKNBFRCLE/wdBVxNgVaQZugAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit\n", "\n", "def func(t, a1, a2, a3, a4):\n", " return a1+a2/(a3+(t-a4)**2)\n", "\n", "xdata = np.linspace(0, 10, 20)\n", "y = func(xdata, 1, 10, 1, 4)\n", "y_noise = 0.2 * np.random.normal(size=xdata.size)\n", "ydata = y + y_noise\n", "plt.plot(xdata, ydata, 'b-', label='data')\n", "\n", "popt, pcov = curve_fit(func, xdata, ydata)\n", "plt.plot(xdata, func(xdata, *popt), 'r-', label='fit')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 8.70589402 42.10001551 400. 90.02564658 126.95109166]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd4VUX6xz+T3hMgIfQiVUBAmqDSlGJHV1HsuLruqr91\n17rqrrrqqqhrW1fdta1YUayIBRFpKqIgHaSHDgkkpPfM74/3Ts7JzU0hyU2dz/PkOff0uTfnfM/3\nvDPzjtJaY7FYLJbmS0BDF8BisVgs/sUKvcVisTRzrNBbLBZLM8cKvcVisTRzrNBbLBZLM8cKvcVi\nsTRzrNBbLBZLM8cKvcVisTRzqhR6pdRrSqlkpdR6r+V/VEr9qpTaoJR63LX8bqXUNqXUZqXUZH8U\n2mKxWCzVJ6ga27wO/Bt4wyxQSo0HpgCDtNb5Sqm2nuX9gGlAf6AD8I1SqrfWuriyE8THx+tu3brV\n6AtYLBZLS2XlypWHtdYJVW1XpdBrrZcopbp5Lb4BmKG1zvdsk+xZPgWY5Vm+Uym1DRgBLKvsHN26\ndWPFihVVFcVisVgsLpRSu6qzXU1j9L2B0Uqp5UqpxUqp4Z7lHYE9ru32epb5KuD1SqkVSqkVKSkp\nNSyGxWKxWKqipkIfBLQGRgJ3AO8rpdSxHEBr/ZLWepjWelhCQpVvHhaLxWKpITUV+r3AR1r4CSgB\n4oF9QGfXdp08yywWi8XSQNRU6D8BxgMopXoDIcBhYA4wTSkVqpTqDvQCfqqLglosFoulZlRZGauU\nehcYB8QrpfYC9wOvAa95mlwWAFdrSWy/QSn1PrARKAJuqqrFjcVisVj8i2oMA48MGzZM21Y3FovF\ncmwopVZqrYdVtZ3tGWuxWCzNHCv0lkbB5s3w7bcNXQqLpXlSnZ6xFovf+cc/YOlSSEpq6JJYLM0P\n6+gtjYLDhyEzs6FLYbE0T6zQWxoFqamQk9PQpbBYmidW6C2NgtRUyMuDYtsY12Kpc6zQWxoFqaky\nzc1t2HJYLM0RK/SWBqekBNLS5LMN31gsdY8Veku9sHYtpKf7Xnf0KJh+e9nZ9Vcmi6WlYIXe4nfy\n82HkSHj6ad/rTdgG6t/Rr1kDmzbV7zktlvrGCr3F72zaJLH3/ft9r3cL/ZEjcPfd9efsBw+Gfv3K\nL7/1Vnjqqfopg8Xib6zQW/zOes9ow25Bd+Ne/s03MGMGLFok8wsWwMyZfi0eAMnJZec//hi+/tr/\n57VY6gMr9Ba/s26dTKsj9Mb1m2X//jc89JD/ymaYP7/svO3AZWlOWKG3+J1jcfQHDpRdlp7u3yaX\nYWEy/eorZ1leHmRlWaG3NB+s0Fv8jnH0R474Xu9L6M226en+q6AtKhJRBydUBOLmwQq9pflghd7i\nV9LTYc8eCApqfI4+K0umISHgHp/eCr2luVGl0CulXlNKJXtGk/Jed5tSSiul4l3L7lZKbVNKbVZK\nTa7rAluaFmvWyHTECHHmeXmwcyfccYd0lAIR9TZt5LOpFHULfX6+s21dkpEh086d5RzG3fsS+oKC\nuj+/xVJfVMfRvw6c4b1QKdUZmATsdi3rB0wD+nv2eUEpFVgnJbU0SX7+WaaTJsk0LQ3GjYN//lME\nH0TUO3aUz27xB0eM/eHqzbG7dJHpnj1w2WVOqKmgQP62bIHISKeuwWJpalQp9FrrJYCvl+6ngTsB\n91iEU4BZWut8rfVOYBswoi4KammarFghjtm0VU9Nhd0ea2DayqemQtu2EkIxHDkiDts46WMV+uRk\nuOWWyp24cexG6BcuhHffhTffLLvNjh0Sz9++/djKYLE0FmoUo1dKTQH2aa3XeK3qCOxxze/1LPN1\njOuVUiuUUitS3AFSS7Pi559h+HBo3VrmV64ERQkRZJcK7a5d0L49REQ4+6Wmlk2ZcKxC/9VX8Mwz\ncr6K8Hb027dDODmsWe14l8xMJ5ZvphZLU+OYhV4pFQHcA9xXmxNrrV/SWg/TWg9LSEiozaEsjZTU\nVBFPI/SxHKX1PX8gnViyiWLA9GGkz/qS/fth4EAJj7j3dQv9sba8MaGfPXsq3sYdox/Pt1zzyink\nEEmWjuAVriWBZDIyrNBbmj41cfQ9gO7AGqVUEtAJ+EUp1Q7YB3R2bdvJs8zSAjFuetgwiFdH+JGR\nnLnvZWYzlQe4D5WTTeylZ/EQf2PQQF3G0R89WrY1zrE6erPv7t0Vb2OEfuwPj/ItpxORdYiH+Btv\ncBVX8iY/cDJ5Ow+UCrxthWNpqhzzmLFa63VAWzPvEfthWuvDSqk5wDtKqaeADkAv4Kc6KquliWHG\nf+3Ts5j2l19AEUlM4BsCTx/PggXQ+W/3MPyNP/K3Hx8me34QERF/L7P/rl3O55oKfWWOPjMTbuVJ\ner9+D+9wKX8KeoXDBfK0eTPwGr4qnkDCbWez6JrlQLB19JYmS3WaV74LLAP6KKX2KqWurWhbrfUG\n4H1gI/AVcJPW2o4Z1EIxoZeE2S8Q9MNSbgx8iUWM5+yzPevzQnmq93+YFT6dyH8+wFl5HwFSMQtS\nCWrwR+gmftV8nuAOCs+/iCt5k8M5zivFga4juYK3iN2+ihO/eQKwoRtL06U6rW4u1Vq311oHa607\naa1f9VrfTWt92DX/sNa6h9a6j9b6S38U2tI0SE+HRJVM6D/+BpMm8UWbKwE46yxZn5EBa9YF8NYp\n/4ERI7hn+2/pwD66d5f1pvkl+CF0c+gQ58++jF9VP4LefJ2gEGkFrJSs7t4dPuV8dg6byrilD9KJ\nPTZ0Y2my2J6xFr+Rng53hjyDysyEZ5+ldRtF+/bQu7e0sDl6FDZsgOMHh8LbbxOsC3iRG+jeTVq9\n1EbozYhVPh291nDDDYQWZHJDm/dRUZHExckq0wzUPGwWnvUEqqSYO3ncOnpLk8UKvaXWFBbCt9/6\nWJ5ylOsKnoepU6FvX846C667TlxzTIzE8AsKoGtXoGdPZp/wEOfxGVNyZwFlhd6Ebl58EQ4erLpM\nxtEnJ0uv1zJ89BF8/DGzBzzIgVai7EboBwyA//s/uFJePtgf3JUl3a7md7xM8OED1fo9LJbGhhV6\nS62ZOxdOP91Jd2AYvu41YnQG3HUXIL1hH3xQ1kVHOx2QTEx+8Yl/ZjkjmPLtzcSTUqaDUm6uiPaN\nN8KsWVWXKTUVoqLk8969rhWZmfCnP8GgQbzb/lZiYmSxEfr4eHjuORgzRjpwZWbCu13+Qhj5jN36\nCkVFUGxrnSxNDCv0llpj+rtt2OBaqDWn73yZ9dGj4MQTy+0TE+MIvelGER4VyG95jdDco/wr7E7A\nEevcXKd5Y1WxcjPY+MCBMl8mTv/gg7BvH7z4IunZQT6F3hAdLefaonsxnwmcue9lJp5WzO23V35+\ni6WxYYXeUmuM8G7e7Fr4ww90zfmVb477nc99oqOdcIxx9BERsJH+HLz8Ni7Ne51TWUq7drIuJ6f6\nHZfS0yUMf9JJMv/kk56EZevWycC1v/sdjBpFRoaUAyA2VqYmuZopo+kZ+19+T7uCPSSu+orVq6v8\nSSyWRoUVekutMR2Pygj966+TpaJY12eqz32MkwZH6E3P2Kw/30tm6668wI1kpRUSECCO3uTGqWo8\nWROfHzQIXngBPv8cHnqgBG64Qaz7o4+Wlrs6jj4rCz5lCkdUPFOy3qpw7FuLpbFihd5Sa8o5+oIC\n+PBDvgg+n7D4KJ/7GCcdEODkwTE9Y2PaR5L+j39xAuu54sgzhIeL0FfX0ZsWN61bi7ZPngz69Znw\n/fek//Vx/v5cGwoLfQu929HHxDhCX0Qw7+uLOI85HN1XTyOXWyx1hBV6S60xjn7LFk+a4QULIC2N\nt4suKQ2JeGMENj4eAj2JrC+6SAYGT0yEjn84jzmcywPq7/QM3eMzdFNS4tvdG0dvHiAXjT7EbQdv\nJ3fIybxWMp0HHoDPPpMQjxF4X0JvHL15kM1iGpHkMD77s9LvbLE0BazQW2qNEcKcHKnn5L330HFx\nfFkyqUKhN46+bVtnWdeu8Je/SPNLpWDAgn8RHqp5LPeP5OboUlE3Qv/887KPt8M3Qt+qlUynff9/\nRJLN+5NeZeOvcsnfequkHp4wQbYxIZvExLJlNEnNgoLgO05lHx2YymwbvrE0KazQW2pNRobTo3TL\nxiL47DNyJ5xLISFVOvrKEpced1o31EMPMjn3U07a+L9yjn7JEslb7x7YG7wc/YcfEvXlB7yYcD9v\nr+xb2jLIpEYeM0bmL71Umtd36uQcp00baZqpNbRrByUEMpdzmMTXHEjybpxvsTRerNBbak1mpjRl\nDAiApHd+gNRU0kZPATgmR++TW2/l56jxTF91M8G7tgGO0JuRoD78sOwupY6+KAVuugmGDGHPJbez\ndKk0ATXu/eKLnbBRdDRccEHZ4/To4YSGTOufuZxDNFkULVxaRcEtlsaDFXpLrcnIkJQBY8ZAwNw5\nEBLC/gEydmCthT4ggEf6zKRIBXPGW5cTRCHZ2RIm2rpVQipz5zrjvYJUxkZHFBN6zWWSZ+F//2Pc\nxGDy8qSst9wCU6ZI56vK6NnT+WyEfpE6jTxCifv+8yoKbrE0HqzQW2pNZqYI99SpcHLqZ2QOH09q\noSh5VaGbKoUeyGrVmceOe4lO+37iP/yBrEzNxo1SGTttmjh8d9v21COapwNug2++gX//GwYOlIeQ\n52ofNQo++URy7lSGW+jbt5dpbIdIlgaOp+sGK/SWpoMVekutMc0Upw7bSR+28EPsWaUpit3t5d1U\n29ED4eHwecRUPh96L9fyGn9NvZW1q2UUcZPyODnZs7HWnLHsfq7NelZSHVwrWbXj4mDoUNmkf//q\nfa/jjnM+G0ffti38GH82bY9ulVcKi6UJYIXeUiu0dhx9wur5AMwrmVgq9BU5etP00TjlyjDt6N87\n/gGe4U/8X+EzDH30QrqGJ5f2fk1JQYLzl1/OJVse4st218BTTzm1xMAll8DgwdV7uIC06+/QoWw5\n27aFX4/zPF0+t67e0jSwQm+pFfn5kr0yJgaYP5+U0I4sSe5bpdAPHQrvvQdnnln1OSIiPB2mshW3\n8DS38BTH7/yCTfnd6XzPldzLgwx/4RqpKJg9m+cSHuSVUa86sRoPt90Gq1Yd2/cz4Ru3o6d7d7YG\nH2+F3tJkqM4IU68ppZKVUutdy55QSv2qlFqrlPpYKRXnWne3UmqbUmqzUmqyvwpuaRyYNvQxkcWw\nYAHbuk9k23ZFerqYaROi8UYpafUSVI3BLMPDpfJVWsAonuEWJsavZmnHaQQtmMeD3E+PjXNY0moK\nn9y7kkcC76V1G1XVYauFEXq3o2/VCr4KPBsWL7bDTlmaBNVx9K8DZ3gtmw8M0FoPBLYAdwMopfoB\n04D+nn1eUEoF1llpLY0O00O0W+ovkJZG2lAJ22zYIOGZgDp4Z/ROgQCwJOV4Zk9+FZKT6dMtn9+e\nd4Txe97g3Q0DSUtzQkO1xVTYdukimTR79pR4/9z8SfIq8913dXMii8WPVGcowSVAqteyr7XWRZ7Z\nHwHTzWQKMEtrna+13glsA0bUYXktjQwj9N23SXxeTZSupp9/DiPq6D9fGrrxMs8mnBLXNoRVq6QV\nztatEk6qK6G//nr44AOJ1W/eLHW7rVrBUn0KOjjY94grFksjoy5i9L8FzNiwHQH34G17PcvKoZS6\nXim1Qim1IsUkNLc0OUzopv36+TB4MF2GSU1nQYHT67S2RESIiB85UjYXjTtubhrAbNwo07oS+lat\n4MIL5XOHDhAcLI4+lwjyh4wqJ/RLlsD06ZJewWJpLNRK6JVSfwWKgLePdV+t9Uta62Fa62EJlfWD\ntzRqMjIggmxiN3wPEyeWaZI4enTdnMPEx/ftc8QdnLw07svHDBtYV0LvC5MA7eiQ06R216TLRJKy\nzZwJb73lv/NbLMdKjYVeKTUdOAe4XGutPYv3AZ1dm3XyLLM0UzIzYQxLCCgqhIkTCQ8X5xsaCsOG\n1c053A8Pd9KxMi1hvPCn0JtkaQePP01eNZYsKV1nBhV/6CEJ4VssjYEaCb1S6gzgTuA8rXWOa9Uc\nYJpSKlQp1R3oBfxU+2JaGisZGTCWxRKvPuUUQJpOnn66iH1d0KOH89mX0Pt6IawPR7+7/UlSU+wK\n35h6hB07yo+ha7E0FFU2blNKvQuMA+KVUnuB+5FWNqHAfCUdUn7UWv9Ba71BKfU+sBEJ6dyktbZD\nKTdjMjNF6EuGDCfQM3LIrFnSkaquaN9eHhr5+ZU7erMNOK7bHxihT8sOkfiUS+jd49lWNbatxVJf\nVKfVzaVa6/Za62CtdSet9ata655a685a68Gevz+4tn9Ya91Da91Ha/1lZce2NH1yD2czjBUEjHNq\nXiMinGEB64KAACckYsQ9MtIZONwI/fDhzj714eg3bIBXdpwG69fDoUOAiLvpG1DVkIcWS31he8Za\nakX81mUEU4QaN9av5zFxeiPq7kpZE7oZNUqmwcF1+6DxxvT2/eQT+O+202Rm4UJAhN5UHtu+VJbG\nghV6S63ouGMJxQTAySf79TwmTh8dLSLuFvouXSAkROoGEhLEzau66Rjrk8BASfmwdSv8whAKwqJh\n0SJA6iyM0FtHb2ksWKG31Ioe+5awMWxIxWkq6wjj6E3Ixi308fGwc6ekSe7Y0b9hG4MJ35QQyM5O\noyUdAmUdvRV6S2PBCr2l5uTl0evIj6yNq6OeUZXQq5dM4+Lgj3+EK68su75DB4nljxwJJ57o9+KU\nqexdEzcWfv0VDh0iM9N5CNnQjaWxUI2UUhZLBfz8MyEl+WxO9G98HuCMMyTb5UknObF4X7z4ot+L\nAjiOHuC7gLFcDOjFS8jMnEqbNlIhax29pbFgHb2l5ixeTAmKpE6n+v1UgYGS7bIukqTVBW6hX5I1\nBCIjKfp2MSUlTj2CdfSWxkIjuW0sTZIlS9gcfAIB8fUQFG9kmNBNTAzs2u/pLOaJ08fESD2CdfSW\nxoIVekvNKCyEH35gacCYCnPON2eMoz/lFBl/vODksQT/up42HLaO3tLosEJvqRm//ALZ2XxbOMbf\nDW4aJUboT/VErQ71kXqK0SwlOto6ekvjwlbGWmpE/vwlhAILS8YwuAU6+vPPl6SVJuf+9tbD6RAa\nztj8xURHX2AdvaVRYR29pUYse2wxB+P6kExii3T0gwbBM89AJ8+QOzv3hZDWZxTjWFQaurGO3tJY\nsEJvOSa0BoqLGZz1HUsDJFzREmP0hu7dRexvvRVWRI5lIGuJLUmzoRtLo8IKvaXazJ0rzRv3fbmW\nONL5IlM6SrVER28IDYWlS2X6yLJxBKCJ37TUhm4sjQor9JZqc+utMk39RAba+KbQCj1At26Sf/8n\nRpBHKFG/LLaO3tKosEJvqRZaO+Oyxq5ZzA66s9czmFhLDt0Yhg6FfML4kZGELFtsHb2lUWGF3lIt\nnNGSNAkbl7AYJ+1BS3f0IEIPsCx4LGrVKloHppOfbwcJtzQOqhR6pdRrSqlkpdR617LWSqn5Sqmt\nnmkr17q7lVLblFKblVKT/VVwS/0yf75Mj2cT4TlHWIKTyMw6ekmkphSsjBoLJSX0TvkesOEbS+Og\nOo7+deAMr2V3AQu01r2ABZ55lFL9gGlAf88+LyilAuustJYG4/BhmY5FuvlbR1+WmBjo3Ru2tB4J\nwcF03y2/kxV6S2OgOkMJLgFSvRZPAWZ6Ps8Eznctn6W1ztda7wS2ASPqqKyWBiQrC8LCYAxLOBzW\nkZ3I2H4BATJ0oAUuvxxGT46AESPotH0RYIXe0jioaYw+UWt9wPP5IGCGbO4I7HFtt9ezrBxKqeuV\nUiuUUitSUlJqWAxLfZGdDQnxmrEsZkXEWECGcIqO9u9oTk2Je++F558Hxo6l9c6VRJFpK2QtjYJa\nV8ZqrTWga7DfS1rrYVrrYQlm0E9LoyUrC/qHbacDB1isbXy+UsaOJaCkmJP5wTp6S6OgpkJ/SCnV\nHsAzTfYs3weeNndCJ88ySxMnOxtOLVoEwJfZIvQRETY+75OTT0YHBjKWxdbRWxoFNRX6OcDVns9X\nA5+6lk9TSoUqpboDvYCfaldES2MgOxtOyl1IcmA71hT0BaSlicn1YnERFUVO/+GMZbF19JZGQXWa\nV74LLAP6KKX2KqWuBWYAE5VSW4EJnnm01huA94GNwFfATVrrYn8V3lJ/ZGVqTjy6iJVR4zDx+Xff\nhXfeadBiNVoKR41lOD+TeySnoYtisVSdplhrfWkFq06vYPuHgYdrUyhL4yPh6Fba5O9nbcfxkC7L\nOnSQIf4s5Qk8bSwh/32M0F+W8fvfn05qKsye3dClsrRUbD56S7UYnLYQgE3txsMOaWppRb5ioiaf\nQjEBxKxazA85p5OcXPU+Fou/sEJvqRYjchZyNLID6Qk9Acm3bqkYFRvDxrAhdN6+iKR8abWUng6x\nsQ1dMktLxOa6sVSJLtGcUriIpG7jiYyS+LztJFU1mxPH0vPIcoqycgEnKZzFUt9Yobf4pLAQzj4b\nfvwRCtduoh2H2Nd7PFFRst46+qo50HssIRRwEsuB6gn9hAnw3nt+LpilxWGF3uKTffvgiy/gq6+g\nYP4iAFIGjC8VeCv0VZMzdDQlqNL8QFu2VL59cTEsWAArV9ZD4SwtChujt/gk1ZPdaM8eCDi4kF10\nobhLd6I8CS5s6KZq2vaOYw2DGMtioqOrdvSmzX1Bgf/LZmlZWEdv8cmRIzLds6uEkGWLWMh4oqKV\ndfTHQOfOkuVzFMs4ZVh+lY7eCr3FX1iht/jEOPrQbRsISjvMIsYRGYmN0R8DXbrAIsYRTh4TYn+2\njt7SYFiht/jECH3vfdJ+fiFSEWsE3oZuqqZTJ1jKaACGZS/m6FHIqaSjrBV6i7+wQm8BID8frrsO\nfvc7mTehm1OLFpIZ353ddLWO/hiJiIDorm3YH38CPfcuAqCyjNxW6C3+wgq9BYDp0+HVV+GVV2Q+\nNRUCKGYsi9nUfjwg4m5j9MfGDz9AwiWn0277d4SRW2kPWSv0Fn9hhd4CwLJlMjWO/cgRGMYKWpPG\nwuBJpevMehu6qR4dOkDw2ZMILMhjNEut0FsaBCv0FgBypfMmWVnSWSo1FabGzKMExezUCYB19DVm\nzBh0SAiT+Lqc0BcXS9gMnPi9FXpLXWOF3gJAXp6TpCw9XYR+svqaX9QwVia1AayjrzGRkZSMOtWn\n0F9yiSSIA8fRG+G3WOoKK/QWQBx9+/by+ehRyE9Op1/Gj6xuK2GbwEAICZFtunSBE05owMI2QQLP\nmsxA1pGz/UCZ5R9+6Hy2oRuLv6iV0CulblFKbVBKrVdKvauUClNKtVZKzVdKbfVMW9VVYS3+obhY\nwjVG6NPSYEDytwTqYoonTAagpEQGAY+Kgl27YNy4hitvk2SSPDDbr5/vc7XWVugt/qPGQq+U6gjc\nDAzTWg8AAoFpwF3AAq11L2CBZ97SiDHxebfQj8qcR15INP1+OxIQIbLUgoEDSQ1qS6+dX5cuKnaN\nvVZYaIXe4j9qG7oJAsKVUkFABLAfmALM9KyfCZxfy3NY/Iy30O/epZmk57G312mMHB3ccAVrTgQE\nsDZxIoMPz5fXI2D3bmd1fr4Veov/qLHQa633Af8EdgMHgHSt9ddAotbaBCIPAom1LqXFr+TlybRD\nB5keXbGN7iSRMmQSwcFw6qkwZEjDla+5sL3HJNoUJcOqVUDZJGd5eVboLf6jNqGbVoh77w50ACKV\nUle4t9Faa8DnS79S6nql1Aql1IqUyroLWvyOt6NP+OlzALJPPQOAJUts6ty64ODgMyhBUTJXfl+3\n0FtHb/EntQndTAB2aq1TtNaFwEfAycAhpVR7AM/UZxcRrfVLWuthWuthCQkJtSiGpbYYoW/TBoKC\noO+2z9hAP8L7HwdIJayl9kQd15YfGUnJJ58BVugt9UdthH43MFIpFaGUUsDpwCZgDnC1Z5urgU9r\nV0SLvzFCHx4OnWPSGZK1hM84lx49GrZczY3ERJjLOQStXgEHDrBtm7POCr3Fn9QmRr8c+AD4BVjn\nOdZLwAxgolJqK+L6Z9RBOS1+xC3054bMI5giFoSfS6KtXalTRo0SoQfgiy/Yv99ZZ2P0Fn9SqxGm\ntNb3A/d7Lc5H3L2liWAqY8PDYXLhZxymDRn9RtqQTR3TtSuU9DuBQzu6kPjZZxw+fC3t2sHBg+Lo\nbQoEi7+wPWMtjqMPLuLUjC/4grPo2SewYQvVTDn7HMVHBeeg588nIzmPTp1kuTt0U1JSto29xVJb\nrNBbSoW+1a/LiClM5TPOpXfvhi1Tc+Wss2BOyTmonBxG5S8sFXp36Aasq7fULVboLaVCH7vkM4oC\ngpnHZCv0fmLQIBmtKz84kil8Ws7Rm3CZFXpLXWKFvhnz668wd2755f/4B0yc6MyL0Gsi533Eji7j\nySTGCr2fiI2FoMgwvo89iwv4mE7tJUZjhD4uTrazQm+pS6zQN2Oefhouu6x8npoNG2D1amc+NxcG\ns5rApO3sHTWVsDDo1at+y9pSUErGkp2ZdRGJJHNi9neAZAzVGlp5UgBaobfUJVbomzEZGZCZKUnK\n3OTlyTr3/FRmowMDOeWJ81m/HmJi6resLYmOHeHDvLPIJYx+mz4AnMHYrdBb/IEV+maMqdxLSiq7\nPC9PhMQMcJGbo5nKbNRppxHaMd52lPIznTpBNlF8yZm0+/5DFCXlhN4OPmKpS6zQN2OM0O/aVXa5\naTefmSnTVnvW0ottcNFF9Ve4FoypgP0oYCpByQcYxTKOHJFl1tFb/IEV+mZMZY4enPDNgE2zKSIQ\nLrig3srWkjFCvzz+bHRoKBfxQanQt24tUyv0lrrECn0zpiqhz8wEtGbw9g/4MWwc2ORy9ULHjjIN\nT4yByZNF6FMkR318vKyzQm+pS6zQN2OysmTqHbox7eYzMoBVq+iQsZn5sVPrtWwtGePo4+NBTZtG\nZ/Zy3J7FgGQQBSv0lrrFCn0zplqO/o03KAwIYVHbi+uzaC0aI/QJCcCUKWQQzaTkNwEburH4Byv0\nzZgqhT61EN55hx/iz6Mo2o7hXl/Ex0NoKLRtC0REMDfsIs7Ln004OdbRW/yCFfpmSkmJZEOMjIT0\ndPkzGKElUs69AAAgAElEQVSPWz4PUlL4vPVVhIc3TDlbIgEBMHs23HKLzH8SfRXRZDGFT63QW/yC\nFfpmionDd+4s06NHnXVG6LsteQPi41kUdoYV+nrm3HPhOBnAi9UxY9hFF67ijUpTICQnl/0/WizV\nxQp9M8WEbYxDNMKvtXTGiSONHhs+hcsuIys/mLCwhimnBULCAniby5nE10RkHATgqack06WbCy+E\nP/+5AQpoafLUSuiVUnFKqQ+UUr8qpTYppUYppVorpeYrpbZ6pjb42wCYFjemuZ5x8abH5TRmEVRc\nAFddRW4u1tE3IKGh8AZXEUgJcZ/OBODnn2Hx4rLb7dsHhw41QAEtTZ7aOvpnga+01n2BQciYsXcB\nC7TWvYAFnnlLPWMcvRF64+hF8DV/4D/sbj0YhgyxQt/AhIXBZvryc/gYot55CYW0qc/JgcJCZ7us\nLJsawVIzaiz0SqlYYAzwKoDWukBrfRSYAsz0bDYTOL+2hbQcO5UJ/SiWMYi1fNX9BlDKCn0DExoq\n04/b3UDgrh1MZH7pOnfyuaws583MYjkWauPouwMpwP+UUquUUq8opSKBRK31Ac82BwE7xHQD4C30\nRiDy8uAGXiSDaL5sdRkgDwEbo284jNCv6HwBOj6BP/Cf0nWm8rWoSP5PVugtNaE2Qh8EDAFe1Fqf\nCGTjFabRWmtA+9gXpdT1SqkVSqkVKSkptSiGxRcVVcYWHjjMxbzPG1xFck4UxcUSHrCOvuEwQh8W\nG4q+5recy2d0ZC/gNIs1/08r9JaaUBuh3wvs1Vov98x/gAj/IaVUewDPNNnXzlrrl7TWw7TWwxJs\njpU6p6LK2IjZrxNKAf9VN5CZCYcPy/KoqPovo0Uwb1NRUaB+fz0BlHAdrwCO0JtMozZGb6kJNRZ6\nrfVBYI9Sqo9n0enARmAOcLVn2dXAp7UqoaVG+IzRFxXR5r0XWMJoDsX3JyMDPvlE1p9+eoMU04Lj\n6KOjQfU4jnkBZ/IH/kMoeaVCbx7c1tFbakJtW938EXhbKbUWGAw8AswAJiqltgITPPOWeqSoqLyj\nz80FPvqIsP07eZpbaNtWXOKsWdC3rwxabWkY3EIP8HzobbTjEFfyphV6S50QVJudtdargWE+Vll/\n2ID07CliD67QTa6Gxx8nu2Nv5uw7jzEJMnbs4sVw//0ylqmlYTBCb8JnP4aP55fcE7mNJ5mXei0Q\nYEM3llphe8Y2M/LzJS3xvn0QFOS4xMRNi2DlSracdzslBJaOCRseDtdc02DFtVDe0YeEKp7gDvqy\nmcSf5wLW0VtqhxX6Zoa73XVkpIh9UBCctORxSExk68grARg/XtLlLlkCXbo0UGEtgFMZWyr0ITCb\nqSTRlRFLngAcoS8udt7WLJbqYoW+meEt9AAjQlbTZ8dXcPPN5JSIqpx/PuzZA0OHNkAhLWXwDt2E\nhkJEdBCvxtzCcfu+gyVLSkM3YF295dixQt+EmTcPdu8uu8wt9EY47iu6j+yQOLjxxlKRsB2kGg/l\nQjchMijJ3Ha/Iy2sHdx3H1mZTneUuojTf/QRvPVW7Y9jaRpYoW8ilJSUzXtSVATnnQfPPlt2u3KO\n/qefmFzwGV/1vx3i4ko7TlmhbzxUJPRhrSN4u8s9sHgxCeu+Ld2+Lhz9iy/CM8/U/jiWpoEV+ibC\ntddKL9d77pH53bslZ7l7QBGgzCt+ZCRw772kBsYzp9vNANbRN0LcHaYATjwRTj4ZYmNhVvTvoFMn\nRs+/D9PJvC6EPi/PhoBaElbo/cgHH8CiRXVzrM2bRcRnzJAKuW3bZLmppDMYRx8fD+MCl8LXX/N6\n27+QViR20dzcxkVaGp6ICJnGxsr0lVckH31sLBxIC2PpmL/Sbf8PnMFXQN2EbqzQtyys0PuRP/5R\nbti6wDh3rUXMt2+X+YqEft5Xmvuy7oR27ZjTSWLzzz0H69aJyNt2842HKVPg5Zel/4Ob2FjYsQNO\nf+e37FDH8Th3EkiRdfSWY8YKvZ9IToaDB8sL8bFQUgJ//zvs3y9CH+D5bx09WrWj77/mHYJX/giP\nPkpAVARHjsDNN8PHH9uwTWMjJgauu678w9cMK1hICLfrJziB9VzHK1boLceMFXo/sW6dTGsj9Dt2\nwAMPSAuJ9HTo2lWWewv9//4HZ54p8xkZEKWyCbnvLzBsGFx1FWFhsHevc1wr9E2DkBDn88dcwCLG\n8hD3UnS49gPH5uU5GU0tzR8r9H5i7VqZmuRiNcFkljx8WATdCH1aWlmh//57+PpreQPIyIB7Qx5D\n7dsnTXICAggPlzcMgxX6pkFamntOcW/UM7ThCJ1ff6jCfZKS5Pm+f3/lxzaOXvtMIm5pblih9xNG\n6Gvj6I3QJyXJ1Ah9amrZGH16uoh8ZiaE793KzfmPw6WXStMNygu7rYhtGvztb/D443DDDTJ/sN1g\nXuE6un76LKxZ43OfxYth5UpYtqzyY+flyTVje9m2DKzQ+4m6FPqdO2XarZtMN26UlhdhYfLGYCpq\nUw+XcOXS31EQEAZPPll6HO9BRY4cqXmZLPVH+/Zwxx1w/PEyHx8PdzGDgqg2EtQvLi63z9atMjXX\nTEW4RxyzNH+s0PuBkhLJDAl1E7oxN61x9Cb+36ePPEjMcHMB/3uVAYcX80L3f4pKeLBC37Tp3Vum\n8fGQRmuWXfYcrFhB0VP/YuBA+PRT0fz8fCekt2NHxccrKnKcvI3Ttwys0PuBjAy56eLiZOru0Xos\nGKHft0+mnTpJy5uNG2W+d2+5YZOToT376fDMHayIHs/inteWOY6NyTdt3EIPsHngVDjnHNR9fyN7\n3XaWLIFbboExY6rn6N3t8K2jbxlYofcDJpTSsaNMa+rqjfMuKZFpXJz8bdki80YADuwrYSZXo4oK\n+WvCS0THlG2nZxy96ZhjaVp06QIjRsApp8h8Xr6CF16gSAXzFldwYHchK1fCTz/B+vWyTWVC7xZ3\nK/Qtg1oLvVIqUCm1Sik11zPfWik1Xym11TNtVfti1p59+2DmzPo5V10JvXH0hthYEfqiIggMhO7d\nZfnNRU8ykW/47qJnWZ/XszTXvMEIvXGEdnzYpkVgICxfDpdfLvN5eUDnzrw+8r+M4kcmLHuotMK+\noEBSXyQlOQbBGyv0LY+6cPR/Aja55u8CFmitewELPPMNzsyZMH167WLm1cV0WurQQaY1rZD1JfSt\nPI/Ndu2ko81QVvAI9/ABF/J932vJyKCc0JvQTZs2kiOnsvitpfFiWkuZ0MtL6ZfwOldz9b6H6bF/\nael248fLNgcP+j6OFfqWR62EXinVCTgbPEPWC1MA451nAufX5hx1RX2O0ONvRw/yEGmlU3mPSzhI\nO67nJQ4fUWRllRd64+jbtIHOnSUzoqXpERAgnajy8kTI162DWwKfYyfdmcUldAk+AMDEibJ9RQ90\n9z1gK2NbBrV19M8AdwLul8RErfUBz+eDQGItz1EnGLGtT6GvC0dvXFxwsHw2jr5TuyKGPjGNTuxl\nKrNJo3VpbvqKHH3r1jUrh6XxEBoq1/D69VLJP3RcNL/hI2JJZ1HbqYwdVcBpp8m2N98Mn39e/hjW\n0bc8aiz0SqlzgGSt9cqKttFaa0xu1fL7X6+UWqGUWpGSklLTYlQbI/T1MbhyVUK/YkXV4ZOSEukY\n1aOHzMfGSi4U4+hv2HMPrVbM5yaeZzkjAadjlclrbnA7ekvTJixMxPmXX2T+/PNhHQO5llfpvu97\nFp14C336wLnnwq+/+q6XskLf8qiNoz8FOE8plQTMAk5TSr0FHFJKtQfwTJN97ay1fklrPUxrPSyh\nHmIJ/nT0338vvRENVYVurrgC7r238mMePSpib1rWmBS2cXHwW15l4uonOHrpH3iV6wBx/ObhYbY1\nWKFvPoSFiVlZuVKuhTFjZPnsgGkU33I7vPACwS88y5w5MHy47zh9fQl9aqoMV2lpeGos9Frru7XW\nnbTW3YBpwLda6yuAOcDVns2uBj6tdSnrAH8K/a23Snd1Q3q6DMhtWrl4O/ojR6CqlxgTn/cW+mGH\nPue//J69/SeR+6gzvFSXLk7HqRNOKHssd2WspWljQjcrV8KQIc7A7h06QOATM+CCC6RR/ezZJCbC\noUPlj+G+BzZtkjTJxxJeTEmR3rqmmW9F3HmnvHFYGh5/tKOfAUxUSm0FJnjmGxx/hm5ycuTPkJ4u\nwmxCKL5SCR+tIgGhT6FfvpzfvHcxaxjExgc+IKq1k97QNLVs3Vp6zLoxjt7G6Js+YWGS02jtWhnY\nPTZWmst27Yq0w3z7bclxdMUVnFq4sEqhnzcP5sxxOuFVh23bJCxk2uxXxIED5RsUWBqGOhF6rfUi\nrfU5ns9HtNana617aa0naK1T6+IctcUIsT8cvXdubyP0kZEy7w7d5OdLW+eqhN688vbtK9MhJStg\n8mTyW7XjbD6nY9/oMh2gjNCffLKTt97QoYMsMw8NS9MlLEzi8wUF4uiVEsEfPtyzQXi45ETo2ZM/\nfH4OJ6QvLXfNu+cPeJpNpPq4Sw8elIfI8uW+96/qXsrKsnUAjYUW0zPWn6Gb/HzfQh8WJgLrdvSm\njX3ZFLTlmTdPYrAnnQQjg1fywLKJ0Lo1od8v5M2v29G/vxg449ZNwjPTe9JNjx7yCj9yZI2/oqWR\nEBbmpCAeOlSmCxaUyWEnMbpvvyWnTRe+5EyOfv59mWP4Enpf1+O2bXLfeIdojkXovd+gV660cfuG\noMUJvT9CNxU5eqXE1fsS+qNHK84FXlICX3wBkydD0I/fMZ8JlETHwcKFBPfoUtpOGsRxRUc7beNP\nPdX3MU19gaVpY5rbtmvnDD0YGFj+LY7ERH6a8S376EjClZNh/vzSVe5rtaBApr4cvUnB4Q5LgnMP\nVSX02dnlt7n4Ynio4nT6Fj/R4oTe+8I7elRyhNSGihw9iBC7QzdG6IuKnBuopATuu8/JT7JqlTjw\n37efAxMnEtW9LdErFznpK11ERcm5zj8fHnkERo2q3XexNG6Cg2U6aVLV4/7GHd+ecSwiO7EHnH02\nvP8+4FugfQm9WeYt9Mfq6N2GJiPDaZVmqT9avNA//7w0USsqkvsgM/PYj+3t6DMyygq9L0cPTpx+\n1y5xOR9/DK+9BmecATepFxj33G+kCc133/kUeXP82Fhx9HffLe7O0nzZ5Ek2Mnly1dsmJsJB2vPp\nrYslbjdtGjz6KHm5orzuOp7KHL1379nqOnpz3Zu3BrOPjdvXPy1C6LWuOHRz6JAs27IFLrkEZs8+\n9mMXFJS9GdyOPjLSt6MHJy5qYq7Z2fDhrEKezLmBf+ubUGecIQHYSvoZREc7nagszZ9du2Q6YULV\n2yZ6+qTvyYyTSp9LLoF77uG8dy8lSmWX6VhX145ea99pR3Jzm6bQr10LvXo13VZELULo3a+PvkI3\n4IjtsTp6b3ej9bE7enPugH17eGT5aVyV8x/4y1+k9YR3N1cvHn5Y/iwtg5kz4coroW3bqrcNC5N0\nGAcPIrX277wDjz1G/43v8z2n0DMoqXRbX5WxtRH6/HxnACxzjxQVybLq5tdZurTxPBTWrpXKaTOg\nUFOjRQi921F7XzgmXmgGz67sIiwsLPsaCs5FXFgoF3FWlsTc3Y6+OkI/hU/48+uD6Jm1mieHvgMz\nZlQrDjNuHIwdW+VmlmbCVVfBG29Uf/synaaUgjvv5D/nfkFXkvj8wIlczHtA9Spj8/PlfqlO6MZ9\nzXs/GKoj3klJElIdNMjniIn1jtEF00qpqdHihN47dGPEtjpCf/nl4qbcuC9acyNA1ZWxpefOyWH4\n6zfxCReQHHkc53Zcxaq+l1b5nSyW6uCrd+yqxDM4M2Elu8L68h7TeIMrKUgpX0Pq7egffFD6aVQl\n2CtWlD2n94OhOo7ejKq2ZQs891zV2/sbK/RNgMocvRF6c2FW5ja2bpUWMW68h2XzJfTucJBb6KOW\nzYcBAzh59Qv8k9v4y+gfWJ/Xs6pojcVSbRITHRNjyMuDQ1E9uPGEpdzHA1zKu3y0baDE8V14V8Zu\n3izjGVTm6DMypOXXU0+VPZ+vaWW4U4SsWVP19v7Gl9Dn5ZX/bRsrLV7ojyV0k54unT3czcW8E0R5\nC310dPnQTbeow7zO1Zz/wiQIDuaWwQu5g3+Slh1CZmaVYXmLpdq0aVN+MPi8PInfh0QE8RD3MSni\ne7JLwqW519SppXba29GnpMhnc4/4Euw9eyQWv8k1FFFthL5t2+qLqdbw5z+X78lbF5jv7E4SN2OG\nhJYq6g/jzZIl8OqrdV+26tDihf5YQjcZGbK/u+bd29GbSi2TNz46Why91kBBASN/fIZVOb25jHeY\nf9LfYM0avswdV1qWvDw71J+l7oiPF6H3NidhYU6yu7TeJzGQNeTf9w+YO1fybvzzn2QelovbLfQl\nJc5bqbmX5s1zeubu3StT95i1NQndGKE//vjqC31mJjz7rO8c/LXFl6Pftk2E34SZquKpp8omP6xP\nWoTQu1sNuIXZfdGa0I37IiwultdGreXPuHXj6ufNK1/p5Evoi4o0he9+AP36cfmKW9gYMZwz2q5i\nVv+HICys9EIxbsE6ektd0aaNXMfp6WIkevYULXcL/XHHQQGhJF/3V8luNnYs3HEHq3N7cxUzycuW\n2lAjvuYaN8L98stw//1yTxihdzvfmjr6yEjpPpKcLD3Fb7658n3cZi05WUJNdYXRELfQm/O5314q\nY+vWmg9CVFtahNBX5OhLnTa+Hf0nn8DgweJWcnPllRRE6D/5RN50vWOR5g0hLg4oKWHQjo9ZwTBC\nLp8KERHcOfAr7jpxHsltB5CWJmUw/3xzc1hHb6krTGrqI0ekS8b27TIfGlpW6METquneHebO5fC7\n80khgZlM5421gyh+/wPSDhc72+HcSwcPyj2WnOwIvRtjrioL+XiTkiLdR0zoZtYsePHFysMk5kGU\nlwd//zucdlr1wypV4cvRm/P9+mvV+5eUyG+fnV3xoO3+pEUJfXh42YvMnUHSV2XsunUyveMO+PJL\nZ/mePfDtt/LZ3Z3bOPogCmk1bxYMGsRZr/yGaDJJfvx1WLWKb4MnExMjD4IdO+Cjj2TfmBjnQWId\nvaWuMDmODh8WoTeEhTkJ8cwoZjfe6FzXBwdMYAQ/cVngewSVFBB4yVQ20I/f8iqZR6SNsRFwI37b\nt/sWem8nX1wsoZ333qu43G6hz8uT9utFRZU/JNxm7cgRabZs+qjUhFNOgccec44J8pAz37s6jr6k\nRI7z5JNOfx7vfgn1QYsS+jZtyoZuvEUayjr6rVudbuJLljjL9+xxHgLuZFLFhw4z5OtH2amOI+jK\nS6G4mJ/+/Db92Mj+iVdDYCAZGSLqISESFpo+XfY16YjBCr2l7nA7eiPiINett6P/4QcZt2TrVk9c\nnwB+7HIxo9tsYu+T75FNJK9yHXM3dOd+/k5M5j60dt5Eqyv0AP/+t2RkqEi4Dx92hB6kwxKUbbXm\njdvRm/vYDLlYE9audc7r1oWDB0Wwq+PoDx6U3/WRR5xlDRG+abZCX1QE48fDN984Qt+6tVwoEyfK\n8H++csJ7C/2IEfLZvPKCCL25AFIPlzCeb5nJVZx8SSfO/u4edoT0hc8+g/XrST/rMooJKm1iaYT+\nmmvg0kslBPTqq3Dmmc7xbejGUlcYR792rcSsTzpJ5rdscYR++HA4/XT417+kj97ttzvhmU6dICs3\nkO1DL2YoK5nEPNZxAvfxIF9u6krRlN9wSs7XKErYsaPy0I1b1M0bdEXjMrgdPThvu5UlRDMO2y30\n3s2hq0tJieiGqY/wFvr0dOkkCZU7elMp7f6eWVnlW0L5m9oMDt5ZKbVQKbVRKbVBKfUnz/LWSqn5\nSqmtnmmruitu9dm/HxYtkuysRuhbtZJwyTffwIcf+r5ozD9Ua7kZ+vWT/YzQBwdLjrHWadt4kHv5\neHV3vuV0zmMOu8Zfwx9P28D/9Z4P55wDAQGl7txb6K+4QnqkT5kCv/1t2XFeraO31BXG0X/2mUzv\nuEOmW7fCsGFS79qqldwTf/yj9EbdscNpWda5s4QaxL0q5jOJSSVf0ZNtvBRzO+q7pXzNZLbTgxM/\nuIeYXevKGRVfjt67YteNccsJCeXTPFXH0efm1l7oc3OlHEagc3MdA3bggPNQOeEEEf6KHlhJSeWX\n/fSTfK+alq0m1MbRFwG3aa37ASOBm5RS/YC7gAVa617AAs+8X/n44/I91sz87t1yoYaFSRjGvGau\nXu37n2MuxsOH5UHQq5f8U3bsgF5s4YnWj/LhnuFsoxd/5WE2cTzTeJf2HGDZVS+yvqRfmSRjbqEv\nKpILJiam/HndN4d19Ja6IjZWwjQrVsi8eXOcNk3eKBctKrt9QoJc+0Y0u3UT5+o9yPhOjuOBsBks\ne38vl/IOW1Qfzlr/ON9nDmRNyQDu4WEGhMiIJb6E3gilL6E3eezdjt5wrI7eV+hm9Wp5qFVWKWrC\nK25Hbwb3OXTI+X3M4C8V9Zg1Qu9OKb12rTxETIK6+qA2g4Mf0Fr/4vmcCWwCOgJTgJmezWYCfh0e\nOCUFfvMbeOmlssvNhblnj1w4kZEi9iZvhlvo3WOpmgtk61YIpIgRxcu4K+c+VhYMYAt9+NOhe4iK\nDuDnqY9xUrvdTCz+iveYRh7hpZWxrVzvMG6h37FDPnfoUP57uMXdOnpLXREQIK6+oECaKkZEyP3w\n1lu+t4+PF6FPTpZr0twbvkQpLw8OpIYyi0u5f8RXdGA/N/I8OjaOh/kb6wr68Ct9OOXj22HxYvKz\ni0r3rczRm3XH6ujdlbHmPt69u3yY5OOPpY7ALP/wQ/j557Lb+BL6zp2d8plzmSE8K6pgTUqS3snj\nx4v7Bye85W4N6G/qJEavlOoGnAgsBxK11ub5dhBIrItzVIR5/fH+Z7odvVvoDWlpToWqW3hbZ+2G\nV16hw5+mkkICJ99+MlftfZjDxHMzz3J07W6Oz1jO8PfvJDO2U5lzViX0xlUNG1b+e1iht/gLE74x\nFf4RERXny0tIkDfPbdvks2mZU6HQe+6zU0+FFNryWuiNbPvfd3RhF492eI4kujH0h3/BuHFce09b\nZnMRN/ACbZI3AdrnW7Vb6MPDy94PlTl678rY9u1lfvVq39uZeoiLLpK6OHfCQiP07tBNXJyUJSXF\nOYYZJqIyoe/WTfrcmIerGUoxO1t+t9//vuLvVFfUWuiVUlHAh8CftdZlnrdaaw34bMmqlLpeKbVC\nKbUixZ3Y4hgxr2be2ffMBbh/v7h7M4arm8WLNEPCNnJV7n95kytIoiurUrvC735Hq19/4BP1G4re\nfo9br0hhPIt4jpuJ7Nu5dH9fscijR6kwdLNihZShX7/y38N9LPeAEBZLbTEVsu6WXVVtu3GjhE3M\ntbhrV/lxDwoK5P4KDpbY/zPPyANi1CjYQxc+7/5/TAmdx4M3H4EPPmBznykM52de4CbWFPZjPx0Y\n/vRl0uNq48bSWIpxy8bNt23rDKF4LI7+5JNl3jt8Y7ZLTXUqeaFs8jT3+BW5uSLk4eFSJrejr0ro\nd+4UoQ8KckK2bqHfs8c/w5t6E1SbnZVSwYjIv6219rQI55BSqr3W+oBSqj3gswOz1vol4CWAYcOG\n1bhbg3H03q+AJnRTXCxxyIsvhrZF+zmHlQz1/I3auox4jsB2OEgiSxnN0+o2nll7Glfe3Z8dOxXX\nXAYRHucfEeEM5QbyluAmO1suRLejDwoScc/MlNfDE0+UZd4YoY+K8jH+p8VSC4yj79On6m2NuCYl\nQf/+jtDv3i2i5u3ATWgiMRH+9CdZprU0H46KEoHO0NFw4YW8v/JCHl2tOY4djGchp/Et525eCNe/\nKzvGxIi1TjmJC0NOom/rk4C2tGsnD5mVKyt29O7mjsbRd+oEXbqUr/R0O3p3OpPbb5dm1J98UrYJ\nZFqaHM8IvanDiIlx7nVfYZjiYvndpk6VeXOPm57w2dlynvp4g6+x0CulFPAqsElr7eofyhzgamCG\nZ/pprUpYBT4dfXExbNvJeWxgCL8wNH8l4z9dSWSGqH8Jij0RfZmTcx4rw0fD6NG88HUPQIGGf/aF\ntescR2Aufu9KVPOPCw4WZ2CajHk7n+houUF++QWuu8739zDHsmEbS11TE0evtVz3RugPHpT7wWSS\nVEq2SUpyQiQGpcSFm3Bpbq606pG4uWIHPdhBD17lOm65TvPU7zfDjz/C8uUUfr+cM9bN4ByKoTfQ\nsSOfdx1IzoBB3LV2IGE7BkFRbwgK4ppr5Bz//rfcX8adG6EPDxdjVZmjN/fs669Lk+uXXxYhdwv9\n0aNyvIgI+U327pWHp/mO4NvRHzggFdmmEtfc4yZElJ0tBrA+Gl/UxtGfAlwJrFNKmSjYPYjAv6+U\nuhbYBVxcuyJWTPqRItS2HZzPBsZs2wiXb5QudL/+yn8970PFBLCJ48k+ZRLfFgzhsQVDyeg+mJ83\nRZH1XxiT6CRBCgyUZ8ShQ/KqamJn5uJ3N4EE5x8UFiYXmQkXtfJqUBodLWGbnBzf8XmzjfuYFktd\n4R2jrwxzrUNZoQdpgRYWJkIaEyPuOinJaZvv5v77xU2vWAFffSWNJQYPLr/d0XQlBevbF6ZP56nH\n4MG7stn45i90Pbgc1qwhdu1aYt/5hpmFhfA68G4o9O/PeZuPZ7Pugz61NxlhvYmkF2FtosjKEoE1\nQj9njgi3ubfcjt68+ffoIfftyy/Ld3ILfWqqfOfwcPl9Vq2Sbd2/jy+hNyEaU4kbGupoDMibQmFh\nI3f0WuvvgIrGoT+9psc9Fna+s4wtjJGZw8D3XSUAPmECt73Wn6CB/fj34gEUBEWS8SGsfAy+XwAj\nEuRHN0mSTBOz+HgR+Z9+kvmBA2ValaM3rXkqE3rTwcoc0xvr6C3+4oILxD22a1f1tu5WLm3bOpWx\nID1oIyNF9OLiROiTk30f17y5hoZKCzZwWp258Q65fvcddOsXSdcrRgOjnRUFBZzX+1cmt1vDTaPX\not9+dJYAAA7XSURBVFevYUj2Ui7gbbgUugJZwKGMjmwo7M0WetPnp15E9e/KZ7obm5Z0ZfiZ8RQV\nq9KGG25Hn5jo3Hu7dpUVevMwcMfoAwNh9GhH6H2Fbsyxze+jlPN2717f2B19g3PCFYPYn/0676zq\nx8Mf9iV1ZzRKSZ3Ov56F20dB0CrofZz8k0xlrHE4BvP6VZHQV+XoQ0PlnOaC8BW6MW12e/Xy/V1M\nGayjt9Q1I0fKX3WIjJTrOT+/vKPv0UPmjxyRa9y0xPEO3bhxN4DIyHDSdhvcQq+1vAFMmuTjQCEh\nHEgYyNxWA/n9o/KA6dYRwsnh3Qe30SVvC7Mf2cL4+M1EHdjCxbxP68/T4HNYAXA2UviOXfmcbiTR\njY7zuxK0pSPj6UCH9PYUJnQAYkhKUmXa2BsDZ4S+oEBi7717Vx66MSEid1+AqChH6OszW22TFvrA\nVjF0uOtq9BNw9H15qkZFyYVYVCQX4JQpUqkE1RN6kIELYmOlMgcqdvRmv7AwEXJTyeLL0YO8wlXU\noiY4WG4w6+gtDYlSThzaW+iPO86Zd5ueyt4UvFu6tWpVXugPHhRHvW+ffB4+3PexYmMlMhsdLdkp\nAXKJYO7ugZx++kAeBfZNdMbUfePZNCb12cXvz0jizot30S4viYC9u2hHEiexnDbLUmGZPAPwnDOL\nCLL/0YGc2A50pgP76UDPjxKZTgJ9tiZQ1CqB44gnhQR69YwmOFgRFFTW0X/4ISxc6Pwu7rckt5Gz\njv4YMcKamio/mnlStm9fdiBl00SrIqE3/5DVqyUCZHqzVcfRa+00k6pI6Hv3rvx7REVZobc0PPHx\n5YU+KKisUXG/tVYm9OaeM8TFiRsGeQjs2CGteZ5/3rnPKqrHiolx4t5vvinTzp2lEtXs4y5LYHwr\nYsa04lMGc9JgqSc4eBBMB90LJmTSI+IAe5bvZ9ZT++HAAT6csZ9uIfvpWLyf4WoF7fV+IhflcCbA\n87LfOZ79S6aHwB3xrCxJIPStBEhKgPh4Ape2IXhtKzqPi2NaZCtCf4oTUYiLIz68FRABKOvojxXT\ney8tTSqAXnpJRNo4eUN1HX1amlOBYtbHxJTvpecWekNERPlXWfOPrKp524QJTksfi6WhcLdfNzH6\nrl1F7H0JfXVDNyD3kWmx06WL5JMCMWSjR8s5Bg3yfSy30dqwQaYTJki6Y9PsMtHVPTM83Pk7csQj\n8h6V79QJ9mVEk0k0md16w2Wy/IPFUhk7cqTkB8rM0Iwfkc36RSm8OuMwHUNSeOTWFBJI4eGbDhOa\nkcLBt1PokZ8CPydBSgrnp6dLOoBvpdmhqUYEWAoUEMxR4kjLbcVR4mg3+yw44/6Kf8Q6oFkIvXHQ\nr7wC69dL5eqf/lS+Y1J1hR7K9pZVSlK8ms4RBndlrKFv3/Lt4Kvr6GfNqny9xVIfmPsgIcF5qzU5\n62sbujF1ZSalgBH6774Tpz9wYNkKYDfeodOICBkxKydHRDwoqGw6E3OcNm2k45I7106fPlLHkJvr\npGkGuccXL4YBA+T+DgpS7EiOIoko8k7oTnBfeONW+c7/fFr2uXGxtDx6+22ZP2tSEcvnp9M+NI2R\nfY/yyhNpEphPS+P1Z45yYFMacRylFTLtEub/jjPNSuiff14uwMsug4cfLr9dTYUenORFbtyO3twQ\nvnq9VtfRWyyNgS5d5F4IDxfnrZQjhrUN3ZghDN25Y/r3F4e+Zw+89lrFxzIPF/NG0Lmzcy/v3Cnr\n3Q8J8zk+XjreGpSSRhGrVsnDYdQoZ13XrlJpvHev3N9hYU7dm6mMhbKmLSKibGVsenYQqbQhNb8N\nfXsCE511Xy+Cdzc53wFg260Vf+e6oln0wTRPca1F5N9+u3yvVag4Rj9qlGQVdscGfSUe88bt6E1u\n6uOPL7+dcSJVOXqLpTFw990S9wYRpCefdPqUeDv6uLjyrt2N9zr3WLVdush0+nT4859lXIbTTqv4\nWOY+Gj9epm6h37GjYqFv00ZSMxjatBHBTk2VljHucI/p3LR+vdzfrVo5LYPCw2VZeHjZezkysqzQ\nuyubvbNvGs1wn9PG6KuJu/KzooockLjbVVc5g4kYOnWSeJxp7wvVE3rzMAkNdWKGvhz9RRfJg8Bk\nurNYGjOxsWVDM7fc4nz2Fvqq2uYbUTcO1j2E4eDBEuacNKni/iXe5QK4+moJpXbp4pi8HTtEfN0P\nFrfQm16zZ54p6RncIR536MY0f05Lc4TefTyl4N13y9b/mYygBncb/IqEvmPH+h0julk4eomlyeeK\nmmaB/HNnzvSdDx7KuoFjdfSmJYEvR9+tm7gkVVH3MoulieAduqmsIhact2hTv+XuzzJpkjQxrI7I\ngzj5adPEOJ16qswbR5+VVT5xoa9Q7SuvSC4bt9BfdJHzuXdv5z6NipKetQajD1OmSN2AwTt043b0\niV65e417N023AwIqrpOoS5qFo1dKnrzZ2dXr5l0RNRV6dxzSVFpZLM2Rmjr644+X1iwmdBMQIMfy\nFWKtiJ49xU0DLF0qU/fQhZWFbgwmxh4SItMpU8o66vBwMWY7d8ryIUOcdRX1galJ6KZjR2e+Pgxg\ns3D0IE/oijJDVhdzYVS3Pbvb0S9eDE8/XbvzWyyNHW9HfyxCb+ZNrLsuBM4t4t6O3lvoW7d2ss+e\ncQbcdJPvyl9jFity9N64QzeFhWXTDnsL/ckny5uIeSOor34zzUbon3gCZsyo3THMRVIdNw9lHf2Y\nMVKhZLE0Z9yO/sILyw5q7wvztmvqroyjryuBM+3kQR4+lTl6t+jGxUnWS3cIx2AeSu4RttzH88Yd\nujHxeaMh3qGbU06R+gXzoKyvlCfNxn+ee27tjxEQIK901RV6d2WsxdISMEIfFgYffFD19omJIpDu\nNCR1KfQgQr53b9Uxem/RrQi3o3dTUesid+jGhG2mTxfRd1f0eu8D9efom43Q1xXh4dUX+qAguPxy\np7mXxdLcueAC6fvj7jleGdOnS6WrMUPh4XD99U7Cr7rAl9CHhjodF03/GO8wSkV4C/0jj8C//lXx\ngEAREdIZq7jYEfqBA+GSSyo+R30nMbRC78X110uNfnWpaJBli6U50r493HNP9bcPDZVmxVrDs89K\nuMdURNYVJrziDt24wyy+QjeVMXCgtBIaMEDm775b/irCvOXk5jqhm6qcuhX6Bubxxxu6BBZL80Mp\nZ/yHusYIeWxs2TcHg3H01cnHb46TlFT987tTFRtHX12hb/KVsUqpM5RSm5VS25RSd/nrPBaLpWXj\nFvqgIPlzC31srOSRqmgYz9riHnzECH1VTr1ZOHqlVCCS1HMisBf4WSk1R2u9sfI9/7+9+wmxqozD\nOP59EHVRIpgh/qtGEMGFmIor0V2lLqZ27kQCNxW1aGG4EXcFtg2MBInITUYuhEgJclVa+Dcxx2lC\nB9OiRa4aq1+Lc64ep3uvjvece+557/OB4Z577jDzPvNzfp77njPvMTObmWKjhwevxGnpNl/eq+Lt\nBGc6ddP0I/qNwFhEjEfEFHAEGK3oe5nZEGs1+tYli8VlFvqhl6mbfh3RV9XolwLXC89v5PvMzEq1\nZUu2GFrrSqB2R/RVajd187BGP39+dlK63ZIpVajtZKyk3cBugGday9iZmc3Q+vVw8uT95/0+oi9O\n3dy582jr18ye/eDyDVWr6oh+Eiheabss33dPRByMiA0RseHp6bduMjN7TEuWPPrfwpShdW5gYiKb\no+/X+jUzUdUR/WlgpaQRsga/g3s36zIzq87Ro/1dc2rVqmxNnAMHsiUOBvG+z5X8OCLib0mvA18C\ns4BDEXGpiu9lZlZUXEO+HyTYvz9bhmV8vLcVdKtS2f97EXEcOF7V1zczGxTbt2c3Nz91Cqam6h7N\n/yWzeqWZWV2kbDVMyI7qB42XQDAzK8GaNdn69u2WPq6bG72ZWUl27ap7BO156sbMLHFu9GZmiXOj\nNzNLnBu9mVni3OjNzBLnRm9mljg3ejOzxLnRm5klThFR9xiQ9BvwSw9fYiHwe0nDGXTOmq5hyjtM\nWaG6vM9GxEPXeR+IRt8rSWciYkPd4+gHZ03XMOUdpqxQf15P3ZiZJc6N3swscak0+oN1D6CPnDVd\nw5R3mLJCzXmTmKM3M7POUjmiNzOzDhrd6CW9JOmKpDFJe+oeT9kkTUi6IOmspDP5vgWSvpJ0NX/s\n8x0yyyPpkKTbki4W9nXMJ+mdvNZXJL1Yz6gfT4es+yRN5vU9K2lb4bXGZgWQtFzS15J+lHRJ0pv5\n/uTq2yXr4NQ3Ihr5QXbT8WvACmAOcA5YXfe4Ss44ASyctu89YE++vQd4t+5x9pBvM7AOuPiwfMDq\nvMZzgZG89rPqztBj1n3A220+t9FZ8wyLgXX59jzgpzxXcvXtknVg6tvkI/qNwFhEjEfEFHAEGK15\nTP0wChzOtw8DL9c4lp5ExDfAH9N2d8o3ChyJiL8i4mdgjOzfQCN0yNpJo7MCRMTNiPgh374DXAaW\nkmB9u2TtpO9Zm9zolwLXC89v0P2H20QBnJD0vaTd+b5FEXEz3/4VWFTP0CrTKV+q9X5D0vl8aqc1\njZFUVknPAc8D35J4fadlhQGpb5Mb/TDYFBFrga3Aa5I2F1+M7H1gspdNpZ4P+IBs6nEtcBM4UO9w\nyifpSeAz4K2I+LP4Wmr1bZN1YOrb5EY/CSwvPF+W70tGREzmj7eBz8ne3t2StBggf7xd3wgr0Slf\ncvWOiFsR8U9E/At8yP2370lklTSbrPF9EhFH891J1rdd1kGqb5Mb/WlgpaQRSXOAHcCxmsdUGklP\nSJrX2gZeAC6SZdyZf9pO4It6RliZTvmOATskzZU0AqwEvqthfKVpNbzcK2T1hQSyShLwEXA5It4v\nvJRcfTtlHaj61n3Gusez3dvIznBfA/bWPZ6Ss60gOzN/DrjUygc8BZwErgIngAV1j7WHjJ+SvaW9\nSzZP+Wq3fMDevNZXgK11j7+ErB8DF4DzZL/8i1PImo9/E9m0zHngbP6xLcX6dsk6MPX1X8aamSWu\nyVM3Zmb2CNzozcwS50ZvZpY4N3ozs8S50ZuZJc6N3swscW70ZmaJc6M3M0vcf9gvnpj/a7SSAAAA\nAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit\n", "\n", "# ndata:=256; \n", "# f1:=t->subs({a=10,b=40000,c=380,d=128},a+b/(c+(t-d)^2) );\n", "# f2:=t->subs({a=10,b=40000,c=380,e=90},a+b/(c+(t-e)^2) );\n", "\n", "def func(t, a1, a2, a3, a4, a5):\n", " return a1+a2*1000/(a3+(t-a4)**2)+a2*1000/(a3+(t-a5)**2)\n", "\n", "xdata = np.linspace(0, 256, 256)\n", "y = func(xdata, 10, 40, 380, 90, 128)\n", "y_noise = 10 * np.random.normal(size=xdata.size)\n", "ydata = y + y_noise\n", "plt.plot(xdata, ydata, 'b-', label='data')\n", "\n", "popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [15,50,400,100,150]))\n", "plt.plot(xdata, func(xdata, *popt), 'r-', label='fit')\n", "\n", "print(popt)\n", "plt.show()\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 8.06534811 44.0829106 423.83792773 90.01430296 126.94177586]\n" ] } ], "source": [ "import scipy.optimize\n", "from numpy import *\n", "\n", "params0=[15,50,400,100,150]\n", "\n", "def fit_func(params,t,y):\n", " a1,a2,a3,a4,a5=params\n", " residual=y-(a1+a2*1000/(a3+(t-a4)**2)+a2*1000/(a3+(t-a5)**2))\n", " return residual\n", "\n", "params, cov=scipy.optimize.leastsq(fit_func,params0,args=(xdata, ydata))\n", "print(params)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "?curve_fit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Gauss-Newton法に関するメモ\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "このGauss-Newton法と呼ばれる非線形最小二乗法は線形問題から拡張した方法として論理的に簡明であり，広く使われている．しかし，収束性は高くなく，むしろ発散しやすいので注意が必要．2次の項を無視するのでなく，うまく見積もる方法を用いたのがLevenberg-Marquardt法である．明快な解説がNumerical Recipes in C(Ｃ 言語による数値計算のレシピ）WilliamH.Press 他著，技術評論社1993にある．\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 課題\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## 一山ピークへのフィット\n", "\n", "以下の256個のデータ\n", "```maple\n", "> ndata:=256; f1:=t->subs({a1=10,a2=40000,a3=380,a4=128},a1+a2/(a3+(t-a4)^2) );\n", "> T:=[seq(f1(i)*(0.6+0.8*evalf(rand()/10^12)),i=1..ndata)]:\n", "> f:=t->a1+a2/(a3+(t-a4)^2);\n", "```\n", "で近似したときのパラメータa1,a2,a3,a4を求めよ．ただし，パラメータの初期値は，ある程度近い値にしないと収束しない．\n", "\n", "## 二山ピークのフィット\n", "以下のように作成したデータ\n", "```maple\n", "> ndata:=256; f1:=t->subs({a=10,b=40000,c=380,d=128},a+b/(c+(t-d)^2) );\n", "> f2:=t->subs({a=10,b=40000,c=380,e=90},a+b/(c+(t-e)^2) );\n", "> T:=[seq((f1(i)+f2(i))*(0.6+0.2*evalf(rand()/10^12)),i=1..ndata)]:\n", "```\n", "を\n", "```maple\n", "> f:=t->a1+a2/(a3+(t-a4)^2)+a2/(a3+(t-a5)^2);\n", "```\n", "$$\n", "f\\, := \\,t\\mapsto {\\it a1}+{\\frac {{\\it a2}}{{\\it a3}+ \\left( t-{\\it a4} \\right) ^{2}}}+{\\frac {{\\it a2}}{{\\it a3}+ \\left( t-{\\it a5} \\right) ^{2}}}\n", "$$\n", "で近似したときのパラメータを求めよ．\n", "```maple\n", "> l1:=listplot(T): display(l1);\n", "```\n", "\n", "![C9_NonLinearFitplot2d8.png](figs/C9_NonLinearFitplot2d8.png)\n", "\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 解答例\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "2. ふた山ピークへのフィット．\n", "```maple\n", "> restart; with(plots): with(LinearAlgebra):\n", "> f1:=t->subs({a=10,b=40000,c=380,d=128},a+b/(c+(t-d)^2) );\n", "> f2:=t->subs({a=10,b=40000,c=380,e=90},a+b/(c+(t-e)^2) );\n", "> T:=[seq((f1(i)+f2(i))*(0.6+0.2*evalf(rand()/10^12)),i=1..256)]:\n", "```\n", "$$\n", "{\\it f1}\\, := \\,t\\mapsto 10+40000\\, \\left( 380+ \\left( t-128 \\right) ^{2} \\right) ^{-1} \\notag \\\\\n", "{\\it f2}\\, := \\,t\\mapsto 10+40000\\, \\left( 380+ \\left( t-90 \\right) ^{2} \\right) ^{-1} \\notag\n", "$$\n", "```maple\n", "> l1:=listplot(T):\n", "> f:=t->a1+a2/(a3+(t-a4)^2)+a2/(a3+(t-a5)^2); \n", " nparam:=5:\n", "```\n", "$$\n", "f\\, := \\,t\\mapsto {\\it a1}+{\\frac {{\\it a2}}{{\\it a3}+ \\left( t-{\\it a4} \\right) ^{2}}}+{\\frac {{\\it a2}}{{\\it a3}+ \\left( t-{\\it a5} \\right) ^{2}}}\n", "$$\n", "```maple\n", "> for i from 1 to nparam do \n", " dfda||i:=unapply(diff(f(x),a||i),x); \n", " end do;\n", "```\n", "$$\n", "{\\it dfda1}\\, := \\,x\\mapsto 1 \\notag \\\\\n", "{\\it dfda2}\\, := \\,x\\mapsto \\left( {\\it a3}+ \\left( x-{\\it a4} \\right) ^{2} \\right) ^{-1}+ \\left( {\\it a3}+ \\left( x-{\\it a5} \\right) ^{2} \\right) ^{-1} \\notag \\\\\n", "{\\it dfda3}\\, := \\,x\\mapsto -{\\frac {{\\it a2}}{ \\left( {\\it a3}+ \\left( x-{\\it a4} \\right) ^{2} \\right) ^{2}}}-{\\frac {{\\it a2}}{ \\left( {\\it a3}+ \\left( x-{\\it a5} \\right) ^{2} \\right) ^{2}}} \\notag \\\\\n", "{\\it dfda4}\\, := \\,x\\mapsto -{\\frac {{\\it a2}\\, \\left( -2\\,x+2\\,{\\it a4} \\right) }{ \\left( {\\it a3}+ \\left( x-{\\it a4} \\right) ^{2} \\right) ^{2}}} \\notag \\\\\n", "{\\it dfda5}\\, := \\,x\\mapsto -{\\frac {{\\it a2}\\, \\left( -2\\,x+2\\,{\\it a5} \\right) }{ \\left( {\\it a3}+ \\left( x-{\\it a5} \\right) ^{2} \\right) ^{2}}} \\notag\n", "$$\n", "\n", "```maple\n", "> g1:=Vector([10,1200,10,125,90]);\n", "```\n", "$$\n", "{\\it g1}\\, := \\, \\left[ \\begin {array}{c} 10\\\\ 1200\\\\ 10\\\\ 125\\\\ 90\\end {array} \\right] \n", "$$\n", "```maple\n", "> guess1:={seq(a||i=g1[i],i=1..nparam)};\n", "```\n", "$$\n", "guess1 := \\{a1 = 10, a2 = 1200, a3 = 10, a4 = 125, a5 = 90\\}\n", "$$\n", "\n", "```maple\n", "> p1:=plot(subs(guess1,f(x)),x=1..256): \n", " display(l1);\n", "```\n", "\n", "![C9_NonLinearFitplot2d9.png](figs/C9_NonLinearFitplot2d9.png)\n", "\n", "\n", "```maple\n", "> df:=Vector([seq(subs(guess1,T[i]-f(i)),i=1..256)]):\n", " Jac:=Matrix(1..256,1..nparam,sparse):\n", " for i from 1 to 256 do \n", " for j from 1 to nparam do\n", " Jac[i,j]:=evalf(subs(guess1,dfda||j(i))); \n", " end do:\n", " end do:\n", " tJac:=(MatrixInverse(Transpose(Jac).Jac)).Transpose(Jac):\n", " g2:=tJac.df; g1:=g1+g2;\n", "```\n", "$$\n", "{\\it g2}\\, := \\, \\left[ \\begin {array}{c} - 0.390553882992161205\\\\ 1584.55290636967129\\\\ 24.9577909601538366\\\\ - 0.0472041829705451138\\\\ - 0.00719532042503852940\\end {array} \\right] \\notag \\\\\n", "{\\it g1}\\, := \\, \\left[ \\begin {array}{c} 13.6348019182603064\\\\ 29567.3667677707381\\\\ 410.545681677467769\\\\ 128.512734548828887\\\\ 90.9223109918718678\\end {array} \\right] \\notag\n", "$$\n", "\n", "```maple\n", "> guess1:={seq(a||i=g1[i],i=1..nparam)};\n", " p1:=plot(subs(guess1,f(x)),x=1..256):\n", " display(l1,p1);\n", "```\n", "$$\n", "guess1 := \\{a1 = 30.251, a2 = 3854.136, a3 = 39.571, a4 = 124.800, a5 = 89.960\\}\n", "$$\n", "\n", "![C9_NonLinearFitplot2d10.png](figs/C9_NonLinearFitplot2d10.png)\n", "\n", "何回か繰り返せば，データ点に近づいてくるはず．\n" ] } ], "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.1" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "colors": { "hover_highlight": "#DAA520", "navigate_num": "#000000", "navigate_text": "#333333", "running_highlight": "#FF0000", "selected_highlight": "#FFD700", "sidebar_border": "#EEEEEE", "wrapper_background": "#FFFFFF" }, "moveMenuLeft": true, "nav_menu": { "height": "12px", "width": "252px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true, "widenNotebook": false } }, "nbformat": 4, "nbformat_minor": 2 }