最小二乗法による多項式近似
$n$ 個の観測点 $(x_i, y_i) (i = 1, 2, ..., n)$ を $m$ 次 $(m < n)$ の多項式
$$ y(x) = a_0 + a_1 x + a_2 x^2 + ... + a_m x^m$$
で近似する(curve fitting)ことを考えます。
近似式と観測点 $y_i$ との残差を $e_i$ とすると
e_i = y_i - y(x_i) = y_i - \sum_{j=0}^{m} a_j x_i^j
となります。ここでの残差の二乗和を $f$ とすると、
f = \sum_{i=1}^{n} e_i^2 = \sum_{i=1}^{n} \Bigl(y_i - \sum_{j=0}^{m} a_j x_i^j\Bigr)^2
となります。この $f$ を最小化するような $a_0, a_1, a_2, ..., a_m$ を決定するのが最小二乗法による多項式近似の目的になります。上式を $a_k$ について偏微分すると
\frac{\partial f}{\partial a_k} = -2 \sum_{i=1}^{n} \Biggl( \Bigl(y_i - \sum_{j=0}^{m} a_j x_i^j \Bigr) x_i^k \Biggr) = 0 \quad (k = 0, 1, ..., m)
すなわち
\sum_{i=1}^{n}y_i x_i^k = a_0 \sum_{i=1}^{n}x_i^k + a_1 \sum_{i=1}^{n}x_i^{k+1} + ... + a_m \sum_{i=1}^{n}x_i^{k+m}
となります。ここで
\sum_{i=1}^{n}y_i x_i^k = T_k \\
\sum_{i=1}^{n}x_i^k = S_k \\
とおけば、
\sum_{j=0}^{m}a_j S_{k+j} = T_k \quad (k=0, 1, ..., m)
という正規方程式になります。これを連立方程式の形で表せば
\left\{
\begin{array}{ll}
S_0 a_0 + S_1 a_1 + S_2 a_2 + ... + S_m a_m = T_0 \\
S_1 a_0 + S_2 a_1 + S_3 a_2 + ... + S_{m+1} a_m = T_1 \\
S_2 a_0 + S_3 a_1 + S_4 a_2 + ... + S_{m+2} a_m = T_2 \\
... \\
S_m a_0 + S_{m+1} a_1 + S_{m+2} a_2 + ... + S_{2m} a_m = T_m \\
\end{array}
\right.
となり、この方程式を解けば未定係数 $a_k (k = 0, 1, 2, ..., m)$ が求められます。
最小二乗法による直線近似
直線近似の場合、 $y_i$は誤差を含んでいるが$x_i$は誤差を含んでいないという前提で、直線 $y = a + b x$ に当てはめることを考えます。
\bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_i , \quad
\bar{y} = \frac{1}{n}\sum_{i=1}^{n}y_i \\
S_x = \sum_{i=1}^{n}(x_i - \bar{x})^2 , \quad
S_{xy} = \sum_{i=1}^{n}(x_i - \bar{x})(y_i - \bar{y})
とすると、
b = \frac{S_{xy}}{S_x} , \quad a = \bar{y} - b \bar{x}
のように求められます。
課題80:多項式近似
以下のように、それぞれ3つの関数 $y = f(x)$, $y = g(x)$, $y = h(x)$ に基づいた3つの現象が存在するものとします。ただし、$y = f(x)$, $y = g(x)$, $y = h(x)$ は未知であるとします。それを観測するための実験は、下記の x_observed
で示された11個の観測点でしか行えず、得られた観測値をそれぞれ fx_observed
, hx_observed
, gx_observed
とします。
import numpy as np
f = lambda x: 1/(1 + np.exp(-x))
g = lambda x: 1.0/(1.0+x**2)
h = lambda x: np.sin(x)
x_observed = np.linspace(-10, 10, 11) # 観測点
fx_observed = f(x_observed) # f(x) の観測値
gx_observed = g(x_observed) # g(x) の観測値
hx_observed = h(x_observed) # h(x) の観測値
$y = f(x)$, $y = g(x)$, $y = h(x)$ を多項式近似してください。ただし、 Scipy
および Numpy.polyfit
は用いないものとします(答え合わせに使うのは良い)。
課題提出方法
-
基本的にGoogle Colaboratoryを用いてプログラミングしてください。どうしても Google Colaboratory を用いることができない場合のみ、Jupyter Notebook または Jupyter Lab を用いてください。
-
課題1つごとに、ノートブックを新規作成してください。1つのノートブックで複数の課題を解かないでください。
-
ノートブックを新規作成すると「Untitled.ipynb」のような名前になりますが、それを「学籍番号・氏名・課題番号」のような名前に変更してください。
-
質問・感想・要望などございましたらぜひ書き込んでください。
-
もし課題を解くにあたって参考になったウェブサイトがあれば、それについても触れてください。
-
課題を計算し終わった ipynb ファイルを提出するときは、指定したメールアドレスに Google Drive で共有する形で授業担当者に提出してください。