Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

This article is a Private article. Only a writer and users who know the URL can access it.
Please change open range to public in publish setting if you want to share this article with other users.

More than 3 years have passed since last update.

Pythonデータ解析お百度参り80:多項式近似

Last updated at Posted at 2020-07-02

最小二乗法による多項式近似

$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 で共有する形で授業担当者に提出してください。


0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?