【初心者向け】Numpy polyfitで多項式近似

はい、承知いたしました。Numpyのpolyfit関数を使った多項式近似について、初心者向けの詳細な解説記事を約5000語で記述します。記事の内容を直接出力します。


Numpy polyfit で学ぶ多項式近似:初心者向け完全ガイド

はじめに

データ分析や科学技術計算において、「与えられたデータ点に最もよく合う曲線を見つける」という作業は非常に重要です。この作業をカーブフィッティングと呼びます。カーブフィッティングは、データの傾向を理解したり、将来の値を予測したり、複雑なデータを簡単なモデルで表現したりするために利用されます。

カーブフィッティングの方法はたくさんありますが、その中でも広く使われているのが多項式近似です。多項式近似は、データ点群を最もよく表すような多項式($y = ax^n + bx^{n-1} + … + cx + d$のような形)を見つけ出す手法です。

Pythonで多項式近似を行う際に、最も手軽で強力なツールの1つが、数値計算ライブラリNumpyに含まれる numpy.polyfit 関数です。この記事では、プログラミングやデータ分析の経験が少ない初心者の方でも理解できるよう、多項式近似の考え方から numpy.polyfit の使い方、そしてその応用までを徹底的に解説します。

この記事を読めば、あなたは以下のことができるようになります。

  • 多項式近似の基本的な考え方を理解する。
  • numpy.polyfit 関数を使ってデータ点にフィットする多項式を見つけ出す。
  • 求めた多項式を使って近似曲線を生成し、データを可視化する。
  • 多項式の「次数」の選び方について考える。
  • polyfit の応用的な機能(重み付き近似、詳細情報の取得)を利用する。
  • 多項式オブジェクト numpy.poly1d を活用する。
  • 多項式近似を使う上での注意点を理解する。

さあ、Numpy polyfit の世界に飛び込みましょう!

1. 多項式近似の基本

まず、多項式近似がどのようなものか、基本的な考え方から見ていきましょう。

1.1. 「近似」とは何か?

「近似」とは、複雑なものや未知のものを、より単純な、または既知のもので大まかに表すことです。データ分析における近似は、多数のデータ点が示す複雑な関係性を、比較的単純な数式(モデル)で表現することを指します。これにより、データの背後にあるパターンやトレンドを抽出しやすくなります。

例えば、時間とともに変化するある物の値を計測したデータがあるとします。これらのデータ点は完全に一直線やきれいな曲線を描くわけではなく、測定誤差や他の要因による「ノイズ」が含まれているのが普通です。多項式近似のような手法は、このノイズの影響をできるだけ排除しつつ、データ点全体が示すおおよその傾向を捉えるための「最適な」曲線を求めます。

1.2. 「多項式」とは何か?

多項式とは、変数(例えば $x$)とその非負の整数べき乗($x^0=1, x^1=x, x^2, x^3, …$)に係数を掛けた項の和で表される式のことです。

例えば、

  • $y = 2x + 1$ (1次多項式)
  • $y = -x^2 + 3x – 5$ (2次多項式)
  • $y = 0.5x^3 + 2x^2 – x + 4$ (3次多項式)

などがあります。一般的に、$n$次の多項式は以下の形で表されます。

$P(x) = a_n x^n + a_{n-1} x^{n-1} + … + a_1 x + a_0$

ここで、$a_n, a_{n-1}, …, a_1, a_0$ は係数と呼ばれる定数です。$n$は多項式の次数と呼ばれ、0以上の整数です。次数が高いほど、多項式はより複雑な形状を表現できます。

  • 次数 0: $y = a_0$ (定数関数、水平な直線)
  • 次数 1: $y = a_1 x + a_0$ (線形関数、直線)
  • 次数 2: $y = a_2 x^2 + a_1 x + a_0$ (2次関数、放物線)
  • 次数 3: $y = a_3 x^3 + a_2 x^2 + a_1 x + a_0$ (3次関数、S字のような曲線など)

1.3. なぜ多項式近似を使うのか?

多項式近似が広く使われる理由はいくつかあります。

  1. シンプルさ: 多項式は比較的単純な関数のクラスであり、微分や積分といった操作も簡単に行えます。これは、分析や計算の観点から非常に扱いやすい性質です。
  2. 柔軟性: 次数を適切に選ぶことで、直線的な関係から複雑な非線形な関係まで、幅広いデータパターンを表現できます。
  3. 普遍性: テーラー展開などの数学的な概念からも分かるように、多くの滑らかな関数は、ある点のごく近傍であれば多項式で非常に良く近似できる性質を持っています。データ全体を近似する場合も、多項式は強力なツールとなります。
  4. 最小二乗法: 多項式近似は通常、最小二乗法という強力な数学的手法を使って行われます。これは、「データ点と近似曲線との間の縦方向の距離(誤差)の二乗の合計」が最小になるような多項式の係数を求める方法です。最小二乗法は数学的に確立されており、信頼性の高い近似を提供します。

1.4. 線形回帰と多項式近似の関係

統計学やデータ分析でよく聞く線形回帰は、実は多項式近似の特殊なケースです。線形回帰は、データ点群を直線で近似する手法です。この直線は1次多項式 $y = a_1 x + a_0$ で表されます。つまり、線形回帰は1次多項式近似に他なりません。

Numpyの polyfit 関数で次数を 1 に指定すると、まさに線形回帰と同じ結果が得られます。

2. Numpy polyfit 関数の紹介

Numpyの polyfit 関数は、与えられたデータ点 $(x_i, y_i)$ の組に対して、次数 deg の多項式 $P(x) = a_{deg} x^{deg} + … + a_1 x + a_0$ を、最小二乗法の意味で最もよくフィットするように係数 $a_{deg}, …, a_0$ を計算してくれる関数です。

具体的には、以下の誤差関数 $E$ を最小にする係数 $a_k$ を求めます。

$E = \sum_{i=1}^{N} (y_i – P(x_i))^2 = \sum_{i=1}^{N} (y_i – (a_{deg} x_i^{deg} + … + a_1 x_i + a_0))^2$

ここで、$N$ はデータ点の数です。

2.1. polyfit 関数の基本的な使い方と引数

numpy.polyfit 関数の定義は以下のようになっています。

python
numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)

引数の意味は以下の通りです。

  • x: データのx座標を含む配列ライクなオブジェクト(例: Numpy配列、リスト)。
  • y: データのy座標を含む配列ライクなオブジェクト(例: Numpy配列、リスト)。x と同じ長さである必要があります。
  • deg: 近似する多項式の次数を指定する整数。
  • rcond: (オプション) 特異値分解 (SVD) における小さい特異値を判定するためのパラメータ。通常はデフォルト値 (None) で問題ありません。計算上の精度に関わる高度な設定です。
  • full: (オプション) ブール値 (True または False)。デフォルトは False
    • full=False (デフォルト): 多項式の係数のみを返します。
    • full=True: より詳細な情報(係数、残差、特異値、ランクなど)を含むタプルを返します。これは近似の質を評価するのに役立ちます。後述します。
  • w: (オプション) 重み。各データ点 $(x_i, y_i)$ の重み $w_i$ を含む配列ライクなオブジェクト。重み付き最小二乗法を行う際に使用します。デフォルトは None(すべての点の重みが1、つまり通常の最小二乗法)。後述します。
  • cov: (オプション) ブール値 (True または False) または文字列 ('unscaled')。デフォルトは False
    • cov=False (デフォルト): 共分散行列は計算されません。
    • cov=True: 係数の推定共分散行列を計算し、返り値に含めます。これは係数の不確実性(推定誤差)を評価するのに使われます。後述します。
    • cov='unscaled': スケールされていない共分散行列を計算します。

2.2. 返り値

polyfit 関数の最も基本的な返り値は、計算された多項式の係数を含むNumpy配列です。係数は次数の高い順に並んでいます。

例えば、deg=2(2次多項式)で近似した場合、返り値 coeffs[a_2, a_1, a_0] のようになります。これは多項式 $a_2 x^2 + a_1 x + a_0$ を意味します。

full=Truecov=True を指定した場合の返り値については、それぞれの引数の説明で詳しく見ていきます。

3. polyfit の使い方(基本編):具体的な例とコード

それでは、実際のデータを使って polyfit の基本的な使い方を見ていきましょう。まずはシンプルな例から始め、次数を変えながら近似を行ってみます。

必要なライブラリをインポートします。データの可視化のために matplotlib.pyplot も使用します。

python
import numpy as np
import matplotlib.pyplot as plt

3.1. 例1: 1次多項式近似(線形回帰)

まずは、直線的な関係を持つデータにノイズを加えたものを用意し、1次多項式(直線)で近似してみましょう。これは線形回帰と同じです。

“`python

データの生成

0から10までのx値

x = np.linspace(0, 10, 50)

真の直線関係 y = 2x + 5 にランダムなノイズを加える

ノイズは平均0、標準偏差2の正規分布に従う

true_slope = 2
true_intercept = 5
noise = np.random.normal(0, 2, size=x.shape)
y = true_slope * x + true_intercept + noise

データをプロット

plt.figure(figsize=(8, 6))
plt.scatter(x, y, label=’Generated Data’, s=20)
plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Generated Data with Noise’)
plt.grid(True)
plt.legend()
plt.show()
“`

上のコードで、直線的な傾向を持つデータ点群が生成され、散布図として表示されます。これらの点に最もよくフィットする直線(1次多項式)を polyfit で求めてみましょう。

polyfit を使います。次数 deg1 を指定します。

“`python

1次多項式で近似 (次数 deg=1)

polyfitは係数を返します

coeffs = np.polyfit(x, y, 1)

print(“Calculated coefficients (slope, intercept):”, coeffs)
“`

出力結果は、おそらく [1.9xxxx, 5.yyyyy] のような2つの数値を含む配列になるはずです。これは、近似によって得られた直線の傾き(a_1)と切片(a_0)です。私たちがデータを生成した際の真の値(傾き=2, 切片=5)に近い値が得られていることが分かります。

3.2. 近似曲線の生成と可視化

polyfit で係数が求まりましたが、この係数を使って近似曲線をプロットするにはどうすれば良いでしょうか?求めた多項式に、プロットしたいxの値を代入して対応するyの値(近似値)を計算すれば良いのです。

Numpyには、polyfit で得られた係数を簡単に扱えるようにするための便利な関数 numpy.poly1d があります。これは、係数配列を「多項式オブジェクト」に変換してくれます。

“`python

求めた係数を使って多項式オブジェクトを作成

coeffs[0] が傾き、coeffs[1] が切片

多項式は y = coeffs[0] * x + coeffs[1]

poly_function = np.poly1d(coeffs)

近似曲線を描画するためのxの値を用意(元のxの範囲で滑らかな曲線を描くため、点の数を増やす)

x_fit = np.linspace(0, 10, 100)

求めた多項式を使って、これらのx_fitに対応するy_fitを計算

y_fit = poly_function(x_fit)

元データと近似曲線を一緒にプロット

plt.figure(figsize=(8, 6))
plt.scatter(x, y, label=’Generated Data’, s=20)
plt.plot(x_fit, y_fit, color=’red’, label=f’Polyfit (deg=1): y = {coeffs[0]:.2f}x + {coeffs[1]:.2f}’) # 係数を小数点以下2桁で表示
plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Data and 1st Degree Polynomial Fit’)
plt.grid(True)
plt.legend()
plt.show()
“`

これで、元データの上に、polyfit によって計算された最適な直線(1次多項式)が描画されます。この直線が、データ全体の傾向をうまく捉えていることが視覚的に確認できます。

3.3. 例2: 2次多項式近似

次に、データが直線ではなく、放物線のような曲線的な傾向を持つ場合を考えます。このような場合は、2次以上の多項式で近似するのが適切です。今回は2次多項式で近似してみましょう。

“`python

データの生成 (2次関数 + ノイズ)

y = -x^2 + 10x + 5 + noise

x_quad = np.linspace(0, 10, 50)
true_coeffs_quad = [-1, 10, 5] # [a_2, a_1, a_0]
y_quad_true = true_coeffs_quad[0] * x_quad**2 + true_coeffs_quad[1] * x_quad + true_coeffs_quad[2]
noise_quad = np.random.normal(0, 5, size=x_quad.shape) # ノイズを少し大きくする
y_quad = y_quad_true + noise_quad

データをプロット

plt.figure(figsize=(8, 6))
plt.scatter(x_quad, y_quad, label=’Generated Data (Quadratic)’, s=20)
plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Generated Data with Quadratic Trend and Noise’)
plt.grid(True)
plt.legend()
plt.show()
“`

このデータに対して、polyfit を使って2次多項式で近似します。deg2 に指定します。

“`python

2次多項式で近似 (次数 deg=2)

coeffs_quad = np.polyfit(x_quad, y_quad, 2)

print(“Calculated coefficients (a_2, a_1, a_0):”, coeffs_quad)
“`

出力される coeffs_quad は3つの数値を含む配列 [a_2, a_1, a_0] になります。これは、$y = a_2 x^2 + a_1 x + a_0$ という多項式の係数です。生成に使った真の係数 [-1, 10, 5] に近い値が得られているはずです。

求めた係数を使って近似曲線を描画します。ここでも poly1d が便利です。

“`python

求めた係数を使って多項式オブジェクトを作成

poly_function_quad = np.poly1d(coeffs_quad)

近似曲線を描画するためのxの値

x_fit_quad = np.linspace(0_quad, 10, 100)

求めた多項式を使ってy_fitを計算

y_fit_quad = poly_function_quad(x_fit_quad)

元データと近似曲線を一緒にプロット

plt.figure(figsize=(8, 6))
plt.scatter(x_quad, y_quad, label=’Generated Data (Quadratic)’, s=20)

多項式の文字列表示を生成

poly_str = f’y = {coeffs_quad[0]:.2f}x^2 + {coeffs_quad[1]:.2f}x + {coeffs_quad[2]:.2f}’
plt.plot(x_fit_quad, y_fit_quad, color=’red’, label=f’Polyfit (deg=2): {poly_str}’)
plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Data and 2nd Degree Polynomial Fit’)
plt.grid(True)
plt.legend()
plt.show()
“`

結果のプロットを見ると、2次多項式がデータ点群の曲線的な傾向をよく捉えていることが分かります。

3.4. より高次の多項式近似

同様の方法で、3次、4次、あるいはそれ以上の高次の多項式で近似することも可能です。polyfitdeg 引数に適切な次数を指定するだけです。

例えば、3次多項式で近似する場合は deg=3 とします。返り値は4つの係数 [a_3, a_2, a_1, a_0] になります。

“`python

例: 3次多項式近似

データの生成 (3次関数 + ノイズ)

x_cube = np.linspace(-5, 5, 50)
true_coeffs_cube = [0.1, -0.5, -2, 10] # [a_3, a_2, a_1, a_0]
y_cube_true = true_coeffs_cube[0]x_cube3 + true_coeffs_cube[1]x_cube*2 + true_coeffs_cube[2]x_cube + true_coeffs_cube[3]
noise_cube = np.random.normal(0, 8, size=x_cube.shape)
y_cube = y_cube_true + noise_cube

3次多項式で近似

coeffs_cube = np.polyfit(x_cube, y_cube, 3)
print(“Calculated coefficients (a_3, a_2, a_1, a_0):”, coeffs_cube)

近似曲線の描画

poly_function_cube = np.poly1d(coeffs_cube)
x_fit_cube = np.linspace(-5, 5, 100)
y_fit_cube = poly_function_cube(x_fit_cube)

plt.figure(figsize=(8, 6))
plt.scatter(x_cube, y_cube, label=’Generated Data (Cubic)’, s=20)
plt.plot(x_fit_cube, y_fit_cube, color=’red’, label=’Polyfit (deg=3)’)
plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Data and 3rd Degree Polynomial Fit’)
plt.grid(True)
plt.legend()
plt.show()
“`

このように、polyfit を使えば、データの形に応じて様々な次数の多項式で近似を行うことができます。

4. 最適な多項式の次数の選び方

多項式近似において、最も重要な決定の1つが「多項式の次数 deg をいくつにするか」です。次数が高すぎても低すぎても問題が生じます。

4.1. 次数が低すぎる場合(アンダーフィッティング)

次数が低すぎる場合、多項式はデータの複雑な傾向を十分に捉えきれません。データ点と近似曲線の間の誤差が大きくなり、データのパターンを単純化しすぎてしまいます。これをアンダーフィッティング (Underfitting) と呼びます。

例えば、本当は2次関数のような形をしているデータに対して、1次関数(直線)で近似しようとする場合などです。直線では曲線的な傾向を表現できないため、多くのデータ点から離れてしまいます。

4.2. 次数が高すぎる場合(オーバーフィッティング)

逆に、次数が高すぎる場合、多項式はデータ点一つ一つに過剰に適合しようとします。これには、本来捉えたいデータの基本的な傾向だけでなく、測定誤差やランダムなノイズまで含まれてしまいます。結果として、得られた曲線は元データ点群を通り過ぎるように、不自然にギザギザしたり、急激な変動を示したりすることがあります。これをオーバーフィッティング (Overfitting) と呼びます。

オーバーフィッティングしたモデルは、元のデータ点に対しては誤差が小さく見えますが、未知の新しいデータに対しては予測性能が非常に悪くなる傾向があります。これは、モデルがノイズに適合してしまったため、真のパターンを捉えられていないからです。

例として、先ほどの2次関数のデータに、敢えて高すぎる次数(例えば10次)で近似してみましょう。

“`python

先ほどの2次関数のデータ (x_quad, y_quad) を使用

高すぎる次数 (deg=10) で近似

coeffs_high_deg = np.polyfit(x_quad, y_quad, 10)

近似曲線の描画

poly_function_high_deg = np.poly1d(coeffs_high_deg)
x_fit_high_deg = np.linspace(x_quad.min(), x_quad.max(), 200) # 滑らかに描くため点を多く
y_fit_high_deg = poly_function_high_deg(x_fit_high_deg)

plt.figure(figsize=(8, 6))
plt.scatter(x_quad, y_quad, label=’Generated Data (Quadratic)’, s=20)
plt.plot(x_fit_high_deg, y_fit_high_deg, color=’red’, label=’Polyfit (deg=10)’)
plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Data and 10th Degree Polynomial Fit (Overfitting Example)’)
plt.ylim(y_quad.min() – 10, y_quad.max() + 10) # Y軸の範囲を適切に調整
plt.grid(True)
plt.legend()
plt.show()
“`

このプロットを見ると、10次多項式は多くのデータ点を通り過ぎるように、複雑で不自然な曲線を描いていることが分かります。特にデータ点が存在しない区間では、急激な変化を見せる可能性があります。これがオーバーフィッティングの典型的な兆候です。

4.3. 最適な次数の選び方のアプローチ

最適な多項式の次数を選ぶための絶対的なルールはありませんが、いくつかの有用なアプローチがあります。

  1. データの形状を目視で確認する: まずはデータをプロットして、どのような形状をしているか観察します。直線的なら1次、放物線のようなら2次、S字のようなら3次など、おおよその見当をつけます。
  2. ドメイン知識を利用する: 解析対象のデータが物理現象などに基づいている場合、その現象がどのような種類の関数で記述されるべきか(例: 物理法則から2乗に比例することが分かっているなど)に関する知識があれば、次数選択の手がかりになります。
  3. 異なる次数で近似し、比較する: いくつかの異なる次数で近似を行い、それぞれの結果を比較します。
    • プロットによる比較: 近似曲線をプロットして、視覚的に最もデータの傾向を捉えつつ、不自然な振動が少ないものを選びます。
    • 残差の確認: polyfit(..., full=True) を使って残差の情報を取得し、残差のパターンや大きさを確認します。良い近似であれば、残差はランダムに分布し、特定のパターン(例えば、残差がx軸に対して曲線を描くなど)が見られないはずです。また、残差の合計平方和(最小二乗法で最小化される値)を比較することもできます。次数を上げると残差平方和は必ず減少しますが、減少率が小さくなったところが目安になることがあります。
    • 評価指標を使う: 近似の性能を定量的に評価するための指標(例: 平均平方根誤差 RMSE – Root Mean Squared Error, 決定係数 R^2 – Coefficient of Determination)を計算し、比較します。ただし、これらの指標は次数を上げると訓練データに対する値は改善する傾向があるため、単に指標が良いものを選ぶのではなく、オーバーフィッティングに注意が必要です。
  4. 検証データや交差検証を使う: データを「訓練データ」と「検証データ」に分割します。訓練データで様々な次数のモデルをフィットさせ、検証データでそれぞれのモデルの予測性能を評価します。検証データに対する誤差が最も小さくなる次数を選ぶことで、オーバーフィッティングを防ぎつつ、汎化性能の高いモデルを選択できます。これは機械学習でよく用いられる手法です。
  5. 情報量規準 (AIC, BIC) など: モデルの複雑さ(次数)とデータのフィットの良さのバランスを考慮した統計的な指標です。これらの指標が最小となるモデルを選択します。polyfit は直接これらの指標を計算しませんが、残差の情報などから自分で計算することは可能です。

初心者の方は、まずデータのプロットを見て次数を推測し、いくつかの異なる次数で近似してみて、プロットや残差を比較するという方法から始めるのが良いでしょう。

5. polyfit の応用的な使い方と引数

polyfit 関数には、より高度な用途に対応するためのオプション引数がいくつか用意されています。これらを使うことで、より柔軟で詳細な分析が可能になります。

5.1. 重み付き最小二乗法 (w 引数)

デフォルトでは、polyfit はすべてのデータ点に対して等しく「フィット」しようとします。しかし、特定のデータ点の信頼性が他の点よりも高かったり、逆に低かったりする場合があります。例えば、測定器の精度が異なる場合や、一部のデータが他の要因によって影響を受けている可能性が高い場合などです。

このような場合に、信頼性の高い点にはより強くフィットさせ、信頼性の低い点にはあまり強くフィットさせないようにしたいことがあります。これを実現するのが重み付き最小二乗法です。polyfitw 引数に、各データ点に対応する非負の重みの配列を指定します。重みが大きい点ほど、近似曲線はその点に近づくように(その点での誤差の二乗が小さくなるように)係数が調整されます。

重み付き最小二乗法では、以下の誤差関数 $E_w$ を最小にする係数を求めます。

$E_w = \sum_{i=1}^{N} w_i (y_i – P(x_i))^2$

ここで、$w_i$ はデータ点 $(x_i, y_i)$ の重みです。

例として、データの端の方にある点の信頼性が低いと仮定して、中央付近の点に重みをつけて近似してみましょう。

“`python

例: 重み付き最小二乗法

先ほどの1次関数のデータ (x, y) を使用

データ点数 N = 50

N = len(x)

中央の点に重みをつけるための重み配列を作成

例えば、中央に近づくほど重みが大きくなるようにする

簡単のため、ガウス関数のような重み分布を生成してみる

mean_x = np.mean(x)

xの値と平均値の差が小さいほど重みが大きくなるように exp(-(x-mean)^2 / (2*sigma^2))

sigma を調整して重みの広がりを制御

sigma = 2.5 # 標準偏差。小さいほど中央に集中、大きいほど平坦に
weights = np.exp(-((x – mean_x)2) / (2 * sigma2))

重みを正規化して、最大値が1になるように調整(必須ではないが分かりやすい)

weights = weights / np.max(weights)

重み付きpolyfitを実行

coeffs_weighted = np.polyfit(x, y, 1, w=weights)

print(“Calculated coefficients (weighted):”, coeffs_weighted)

通常のpolyfitの結果と比較するために再計算

coeffs_ordinary = np.polyfit(x, y, 1)
print(“Calculated coefficients (ordinary):”, coeffs_ordinary)

重み付き近似曲線の描画

poly_function_weighted = np.poly1d(coeffs_weighted)
y_fit_weighted = poly_function_weighted(x_fit) # x_fit は最初の1次近似で使ったもの

通常の近似曲線の描画

poly_function_ordinary = np.poly1d(coeffs_ordinary)
y_fit_ordinary = poly_function_ordinary(x_fit)

plt.figure(figsize=(10, 7))
plt.scatter(x, y, label=’Generated Data’, s=20, alpha=0.6)

重みの大きさを点のサイズや色で表現するのも分かりやすい

plt.scatter(x, y, s=weights*50 + 10, c=weights, cmap=’viridis’, label=’Weighted Data Points’) # 例:重みに応じて点のサイズを変える

plt.plot(x_fit, y_fit_weighted, color=’red’, linestyle=’-‘, label=f’Weighted Polyfit (deg=1): y = {coeffs_weighted[0]:.2f}x + {coeffs_weighted[1]:.2f}’)
plt.plot(x_fit, y_fit_ordinary, color=’blue’, linestyle=’–‘, label=f’Ordinary Polyfit (deg=1): y = {coeffs_ordinary[0]:.2f}x + {coeffs_ordinary[1]:.2f}’)

plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Data and Weighted vs Ordinary 1st Degree Polynomial Fit’)
plt.grid(True)
plt.legend()
plt.show()
“`

この例では、データの中心に近い点に重みを与えたため、重み付き近似直線は、通常の近似直線と比較して、中央の点により「引っ張られる」ように変化する可能性があります。重み付き近似は、データの信頼性にばらつきがある場合に非常に有効な手法です。

5.2. 詳細情報の取得 (full=True 引数)

polyfitfull 引数を True に設定すると、多項式の係数だけでなく、最小二乗計算に関する詳細な情報を含むタプルが返されます。これは、近似がどれだけデータに「フィット」しているかを定量的に評価したり、計算上の問題がないかを確認したりするのに役立ちます。

full=True を指定した場合、polyfit は以下の要素を含む5つのタプルを返します。

python
(coeffs, residuals, rank, sv, rcond)

それぞれの意味は以下の通りです。

  • coeffs: 通常の返り値と同じ、計算された多項式の係数配列。
  • residuals: 残差平方和。これは、データ点と近似曲線との間のy方向の距離の二乗の合計です。多項式近似(最小二乗法)はこの値を最小にするように係数を求めているため、この値が小さいほど、近似曲線がデータ点に「近く」フィットしていることを意味します。ただし、次数を上げればこの値は必ず小さくなるため、単独で近似の良し悪しを判断するのではなく、異なる次数のモデルや他の指標と組み合わせて評価することが重要です。w 引数を使った場合は、重み付き残差平方和になります。データ点が少なかったり、次数が高すぎて完全にフィットしてしまったり(例: N個の点に対してN-1次以上の多項式)した場合は [] のように空の配列が返されることがあります(残差が0になるか、一意に定まらない場合)。
  • rank: 設計行列(多項式の係数を求めるための連立一次方程式を表現する行列)のランク。ランクが次数+1(求めたい係数の数)と等しい場合、係数は一意に定まります。ランクが不足している場合は、データ点が少なすぎるか、データが線形従属になっている可能性があり、係数が安定しない可能性があります。
  • sv: 設計行列の特異値 (Singular Values)。特異値分解 (SVD) は polyfit が内部的に係数を計算するために使う手法の一つです。これらの値を見ることで、計算の安定性などを評価できます。小さな特異値がある場合は、データに数値的な問題(多重共線性など)がある可能性を示唆します。
  • rcond: 実際に使用された rcond の値。

最も頻繁に利用するのは residuals でしょう。これを使って、異なる次数のモデルのフィットの良さを比較する例を示します。

“`python

例: full=True を使って残差を確認し、異なる次数を比較

先ほどの2次関数のデータ (x_quad, y_quad) を使用

1次、2次、3次で近似し、残差を取得

coeffs1, residuals1, rank1, sv1, rcond1 = np.polyfit(x_quad, y_quad, 1, full=True)
coeffs2, residuals2, rank2, sv2, rcond2 = np.polyfit(x_quad, y_quad, 2, full=True)
coeffs3, residuals3, rank3, sv3, rcond3 = np.polyfit(x_quad, y_quad, 3, full=True)
coeffs10, residuals10, rank10, sv10, rcond10 = np.polyfit(x_quad, y_quad, 10, full=True)

print(f”Degree 1 – Residuals: {residuals1[0]:.2f}”)
print(f”Degree 2 – Residuals: {residuals2[0]:.2f}”) # residualsは要素1つの配列なので [0] で値を取り出す
print(f”Degree 3 – Residuals: {residuals3[0]:.2f}”)
print(f”Degree 10 – Residuals: {residuals10[0]:.2f}”) # または [] になる可能性あり
“`

この例の出力を見ると、次数を上げるごとに残差平方和が減少していくことが分かります。特に1次から2次への減少幅は大きいでしょう。これは、データが2次的な傾向を持っていることを示唆します。2次から3次への減少幅は小さくなるかもしれません。そして10次ではさらに小さくなるか、あるいはデータ点数によっては完全にゼロ(または非常に小さい値)になる可能性もありますが、これはオーバーフィッティングの兆候である可能性があります。

残差平方和だけでなく、それぞれの近似曲線をプロットして視覚的に確認することと合わせて判断することが重要です。

5.3. 係数の共分散行列 (cov=True 引数)

統計的な分析において、求めた係数の「不確実性」や「推定精度」を知りたい場合があります。polyfitcov 引数を True に設定すると、計算された係数の推定共分散行列が返り値に含まれます。

cov=True を指定した場合、返り値は以下の2つのタプルになります。

python
(coeffs, cov_matrix)

  • coeffs: 通常の返り値と同じ、計算された多項式の係数配列。
  • cov_matrix: 係数の推定共分散行列。これは、係数配列 coeffs の各要素に対応する正方行列です。行列の対角成分は、それぞれの係数の推定分散を表します。推定分散の平方根は、その係数の標準誤差 (Standard Error) と呼ばれ、係数推定値のばらつきの目安となります。非対角成分は、異なる係数間の共分散を表します。

共分散行列を使って、係数の標準誤差を計算してみましょう。

“`python

例: cov=True を使って係数の標準誤差を計算

先ほどの1次関数のデータ (x, y) を使用

1次多項式で近似し、係数と共分散行列を取得

coeffs1_cov, cov_matrix1 = np.polyfit(x, y, 1, cov=True)

print(“Calculated coefficients (deg=1):”, coeffs1_cov)
print(“Covariance Matrix (deg=1):\n”, cov_matrix1)

共分散行列の対角成分が分散

variances1 = np.diag(cov_matrix1)

標準誤差は分散の平方根

std_errors1 = np.sqrt(variances1)

print(“Standard Errors of coefficients (deg=1):”, std_errors1)

結果の解釈: 標準誤差は、それぞれの係数(傾き、切片)の推定値が、真の値に対してどの程度ばらつくと考えられるかを示します。

標準誤差が小さいほど、係数の推定精度が高いと言えます。

“`

この例では、傾きと切片それぞれの標準誤差が表示されます。これらの値は、推定された傾き(coeffs1_cov[0])や切片(coeffs1_cov[1])が、真の値の周りにどの程度ばらつくかを示唆します。例えば、推定された傾きが1.9で標準誤差が0.1だった場合、真の傾きは1.9 ± 0.1 の範囲(より正確には、統計的な信頼区間を計算する必要がありますが、標準誤差はその計算に使われます)に存在する可能性が高い、といった解釈ができます。

cov=True は、単にデータにフィットするだけでなく、そのフィットの統計的な信頼性まで評価したい場合に有用です。ただし、共分散行列や標準誤差の解釈には統計学の知識が必要になります。

5.4. 数値安定性 (rcond 引数)

rcond 引数は、polyfit が内部的に使用する特異値分解 (SVD) における閾値を制御します。非常に小さい特異値は、計算の不安定性(例えば、データが線形従属に近い、つまり多重共線性がある場合など)を示唆します。rcond は、これ以上小さい特異値はゼロとみなすための閾値です。

通常、rcond=None としておけば、Numpyが適切なデフォルト値(浮動小数点数の精度に基づく値)を自動的に設定してくれます。ほとんどの場合はこのデフォルトで問題ありません。

しかし、データに強い相関があったり、次数が高すぎたり、データ点数が少なすぎたりして、polyfit が計算上の警告 (RankWarning など) を出す場合があります。これは、係数が安定して計算できていない可能性があることを示しています。このような警告が出た場合に、rcond の値を調整することが考えられますが、これは数値計算の専門的な知識を要するため、初心者のうちは基本的にデフォルトのままで良いでしょう。警告が出た場合は、次数を見直す、データ点数を増やす、データの前処理を行うなどの対策を検討する方が現実的です。

6. numpy.poly1d の活用

セクション3で簡単に触れましたが、numpy.poly1dpolyfit で得られた係数を非常に扱いやすい「多項式オブジェクト」に変換してくれる便利な関数です。これを使うことで、多項式の値の計算や、微積分といった操作を直感的に行うことができます。

6.1. poly1d オブジェクトの生成

poly1d は、係数配列を引数に取って多項式オブジェクトを生成します。デフォルトでは、係数は次数の高い順に並んでいると解釈されます。

“`python

例: poly1d オブジェクトの生成

係数配列 (例: 2次多項式 ax^2 + bx + c の係数 [a, b, c])

coeffs = np.array([-1, 10, 5])

係数から多項式オブジェクトを作成

poly_obj = np.poly1d(coeffs)

print(“Polynomial object:”)
print(poly_obj) # 多項式のきれいな表現が出力される

係数が次数の低い順の場合 (例: [c, b, a]) は order=’asc’ を指定

coeffs_asc = np.array([5, 10, -1])
poly_obj_asc = np.poly1d(coeffs_asc, variable=’z’, order=’asc’) # 変数名も指定可能
print(“\nPolynomial object (ascending order):”)
print(poly_obj_asc)
“`

出力を見ると、poly1d が係数から多項式を構築し、人間が読める形式で表示してくれることが分かります。order 引数は、係数配列が次数の高い順 (デフォルト) か低い順かを示します。polyfit の返り値は次数の高い順なので、通常は order 引数は不要です。variable 引数で変数名を変更することもできます。

6.2. poly1d オブジェクトを使った値の計算

poly1d オブジェクトは、まるで普通の関数のように使うことができます。多項式に特定のx値を代入して、対応するyの値(関数値)を計算できます。これは単一の値でも、Numpy配列でも可能です。

“`python

例: poly1d オブジェクトを使った値の計算

先ほどの 2次多項式オブジェクト poly_obj を使用

poly_obj = -1 x^2 + 10 x + 5

特定のx値での値を計算

y_at_x0 = poly_obj(0) # x=0 での値
y_at_x5 = poly_obj(5) # x=5 での値

print(f”\nValue at x=0: {y_at_x0}”)
print(f”Value at x=5: {y_at_x5}”)

Numpy配列のx値に対するy値をまとめて計算

x_values = np.array([1, 2, 3, 4])
y_values = poly_obj(x_values)

print(f”Values at x={x_values}: {y_values}”)

近似曲線の描画に使った np.linspace と組み合わせて使うのが一般的

x_for_plot = np.linspace(0, 10, 100)
y_for_plot = poly_obj(x_for_plot) # これで近似曲線のy座標の配列が一気に計算できる
“`

このように、poly1d オブジェクトに関数を呼び出すようにxの値を渡すだけで、簡単に対応するyの値が得られます。これは、近似曲線を生成したり、特定の点での予測値を計算したりする際に非常に便利です。

6.3. poly1d オブジェクトを使った微積分

多項式の優れた性質の一つは、微分や積分が非常に容易なことです。poly1d オブジェクトは、この性質を利用するためのメソッド .deriv().integ() を提供しています。

  • .deriv(m=1): 多項式を m 階微分した新しい poly1d オブジェクトを返します(デフォルトは1階微分)。
  • .integ(m=1, k=0): 多項式を m 階積分した新しい poly1d オブジェクトを返します(デフォルトは1階積分)。k は積分定数です。

“`python

例: poly1d オブジェクトを使った微積分

先ほどの 2次多項式オブジェクト poly_obj を使用

poly_obj = -1 x^2 + 10 x + 5

print(“Original Polynomial:”)
print(poly_obj)

1階微分

(-1 * 2)x^(2-1) + (10 * 1)x^(1-1) + 0 = -2x + 10

poly_deriv1 = poly_obj.deriv()
print(“\n1st Derivative:”)
print(poly_deriv1)

2階微分

(-2 * 1)x^(1-1) + 0 = -2

poly_deriv2 = poly_obj.deriv(m=2)
print(“\n2nd Derivative:”)
print(poly_deriv2)

1階積分 (積分定数 k=0)

(-1/3)x^3 + (10/2)x^2 + 5x + 0 = -0.333x^3 + 5x^2 + 5x

poly_integ1 = poly_obj.integ()
print(“\n1st Integral (k=0):”)
print(poly_integ1)

1階積分 (積分定数 k=10)

poly_integ1_k10 = poly_obj.integ(k=10)
print(“\n1st Integral (k=10):”)
print(poly_integ1_k10)
“`

poly1d オブジェクトは微積分にも対応しているため、多項式で近似した関数を使って、元のデータの変化率(微分)や累積値(積分)といった量を計算することも可能です。これは、物理現象の分析などで役立つことがあります。

7. 多項式近似の注意点

多項式近似は強力なツールですが、使用する際にはいくつかの注意点があります。

7.1. 外挿 (Extrapolation) の危険性

外挿とは、データが存在する既知の範囲(データが得られたx軸の範囲)を超えた領域で、求めた近似モデルを使って値を予測することです。多項式近似で得られた曲線は、データの範囲内ではうまくフィットしていても、データ範囲外では全く異なる振る舞いをすることが非常に多いです。特に高次の多項式は、データ範囲の端を過ぎると急激に上昇または下降したり、不自然な振動を始めたりする傾向があります。

“`python

例: 外挿の危険性を示す

先ほどの2次関数のデータ (x_quad, y_quad) と 2次近似結果 (poly_function_quad) を使用

データ範囲は x=[0, 10]

外挿を含む広い範囲で近似曲線を描画

x_extrapolate = np.linspace(-5, 15, 100) # データ範囲外 (-5から0, 10から15) を含む
y_fit_extrapolate = poly_function_quad(x_extrapolate)

plt.figure(figsize=(8, 6))
plt.scatter(x_quad, y_quad, label=’Generated Data (Quadratic)’, s=20, alpha=0.6)
plt.plot(x_extrapolate, y_fit_extrapolate, color=’red’, label=’2nd Degree Polyfit’)

元のデータ範囲を垂直線で示す

plt.axvline(x_quad.min(), color=’gray’, linestyle=’–‘, label=’Data Range’)
plt.axvline(x_quad.max(), color=’gray’, linestyle=’–‘)

plt.xlabel(‘X-axis’)
plt.ylabel(‘Y-axis’)
plt.title(‘Polynomial Fit and Extrapolation’)
plt.grid(True)
plt.legend()
plt.show()
“`

このプロットを見ると、データが存在する [0, 10] の範囲では近似曲線はデータの傾向をよく捉えていますが、範囲外(特に x<0x>10)では曲線が急激に元のトレンドから外れていく様子が分かります。

したがって、多項式近似で得られたモデルをデータが存在しない領域での予測に使うのは、非常にリスクが高いということを常に意識しておく必要があります。外挿が必要な場合は、その挙動がデータの性質やドメイン知識と矛盾しないか慎重に検討するか、他の予測手法を検討すべきです。

7.2. データの前処理

多項式近似に限らず、データ分析全般に言えることですが、適切なデータ前処理が結果の質に大きく影響することがあります。

  • 外れ値 (Outliers): 最小二乗法は外れ値に非常に敏感です。データの中に他の点から大きく離れた点(外れ値)があると、近似曲線がその外れ値に強く引っ張られ、他の大部分のデータ点からのフィットが悪くなることがあります。必要に応じて、外れ値を特定し、除去または適切に処理することを検討する必要があります。
  • スケーリング: xの値の範囲が非常に広い場合、高次の $x^n$ の項が非常に大きな値になり、数値計算上の精度問題を引き起こす可能性があります。このような場合、x値を正規化(例えば、平均0、標準偏差1になるようにスケーリング)してから polyfit を行い、得られた係数を元のスケールに戻す、といった前処理が有効な場合があります。

7.3. モデルの解釈性

次数が低い多項式(1次や2次)は、その係数が「傾き」「曲がり具合」といった直感的な意味を持つ場合があり、モデルを解釈しやすいです。しかし、次数が高くなるにつれて、個々の係数が持つ物理的・意味的な解釈は難しくなっていきます。単純にデータによくフィットさせるためだけに高次多項式を使うと、得られたモデルがどのような「法則」や「関係性」を表しているのかが不明瞭になることがあります。

モデルを使って何らかの洞察を得たい、背後にあるメカニズムを理解したいという場合は、フィットの良さだけでなく、モデルの解釈性も考慮して次数を選ぶことが重要です。

8. まとめと次のステップ

この記事では、Numpyの polyfit 関数を使った多項式近似について、その基本的な考え方から具体的な使い方、応用的な機能、そして注意点までを詳しく見てきました。

  • 多項式近似は、データ点群を多項式で近似する強力な手法であり、特に最小二乗法によって最適な係数を求めます。
  • numpy.polyfit(x, y, deg) は、指定した次数 deg の多項式の係数を返します。
  • 得られた係数は numpy.poly1d を使って多項式オブジェクトに変換すると、値の計算や微積分が容易になります。
  • 多項式の次数の選択は、アンダーフィッティングやオーバーフィッティングを防ぐ上で非常に重要です。データの視覚化、残差分析、必要に応じて検証データなどを用いて慎重に判断する必要があります。
  • polyfitw 引数で重み付き最小二乗法、full=True で残差などの詳細情報、cov=True で係数の共分散行列を取得できます。
  • 多項式近似で得られたモデルは、データ範囲外での外挿に使うと危険性が高いことを理解しておく必要があります。また、適切なデータ前処理やモデルの解釈性も考慮すべき点です。

polyfit は非常に手軽に多項式近似を実行できる素晴らしいツールです。しかし、ただ関数を使うだけでなく、その背後にある最小二乗法の考え方、次数選択の重要性、そして結果の解釈と注意点を理解することが、データ分析においては不可欠です。

次のステップ

この記事で多項式近似の基本は理解できたはずです。さらにデータ分析のスキルを深めたい場合は、以下のトピックに進んでみるのが良いでしょう。

  • 他の回帰手法: 多項式近似はシンプルなモデルですが、データの関係性が多項式ではうまく表現できない場合もあります。線形回帰モデルを一般化した一般化線形モデル (GLM) や、より柔軟なスプライン補間、あるいは機械学習で使われる決定木サポートベクターマシンニューラルネットワークといったより複雑なモデルを使った回帰手法を学ぶことで、様々なデータに対応できるようになります。
  • モデル評価の指標と手法: 残差平方和だけでなく、R^2、調整済みR^2、AIC、BICといったモデル評価指標について詳しく学ぶことや、交差検証 (Cross-validation) といった評価手法を習得することは、適切なモデル選択に非常に役立ちます。
  • 正則化 (Regularization): 高次のモデルでオーバーフィッティングを防ぐための手法として、Ridge回帰やLasso回帰といった正則化された回帰手法があります。これらは、Numpy単独ではなく、Scikit-learnのような機械学習ライブラリで提供されています。
  • Numpy以外のライブラリ: より高度な回帰分析や統計モデリングには、Scikit-learn(機械学習)、Statsmodels(統計モデリング)といったライブラリが非常によく使われます。これらのライブラリでは、多項式近似だけでなく、より多様なモデルや統計的な検定などが提供されています。

多項式近似と numpy.polyfit は、データ分析への入り口として非常に良いテーマです。この記事が、あなたのデータ分析スキル習得の一助となれば幸いです。


コメントする

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

上部へスクロール