これなら分かる応用数学教室 最小二乗法 p4 サンプル

lsq.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# 参考: http://kaiseki-web.lhd.nifs.ac.jp/documents/Python/leastmean.htm

from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt

# これなら分かる応用数学教室 p.4
# 【例1.2】 4点(4,-17),(15,-4),(30,-7),(100,70)に直線を当てはめよ

# numpy.array
# -------------
# list(Python組み込み) に比べて下記の利点
# - 多次元 array 用の Python 拡張パッケージ
# - ハードウェア寄り(効率)
# - 科学技術計算のために設計(便利さ)

x = np.array([4, 15, 30, 100])
y = np.array([-17, -4, -7, 70])


# 最小二乗法
# --------------
# 一次式の場合: A=[[x0,1],[x1,1],...,[xm,1]] の配列を作成する
A1 = np.array([x, np.ones(len(x))])
A2 = np.array([x**2,  x, np.ones(len(x))])
A3 = np.array([x**3, x**2,  x, np.ones(len(x))])

A1 = A1.T  # 行と列の入れ替え
A2 = A2.T  # 行と列の入れ替え
A3 = A3.T  # 行と列の入れ替え
print("A1=", A1)


a1, b1 = np.linalg.lstsq(A1, y)[0]
a2, b2, c2 = np.linalg.lstsq(A2, y)[0]
a3, b3, c3, d3 = np.linalg.lstsq(A3, y)[0]

# データをプロット
plt.plot(x, y, "o")

# 求めた N次線を上書き
px = np.arange(-20, 120, 0.1)
plt.plot(px, (a1 * px + b1), label="n=1")
plt.plot(px, (a2*px**2 + b2 * px + c2), label="n=2")
plt.plot(px, (a3*px**3 + b3*px**2 + c3*px + d3), label="n=3")
plt.legend()
plt.show()