Python入門トップページ

« 戻る 次へ »

目次

  1. データを眺める
  2. モデル1：単回帰分析
  3. モデル2：2次式を加える
  4. モデル3：3次式を加える
  5. モデル4：sinカーブを加える
  6. モデル5：cosカーブを加える

線形回帰モデル

モデル1：単回帰分析

前のページで眺めたデータについて，単回帰分析を実施してみよう．ここで使ったように StatsModels 等のパッケージを使って回帰分析を行うことも可能であるが，次のモデルに発展させるために，ここでは残差2乗和を自身で計算し，これを最小化するパラメータをネルダーミード法やBFGS法などの非線形最適化手法を用いて求めてみる．

ライブラリのインポート

必要なライブラリを最初に読み込んでおく．Pandas, Matplotlib 以外に，NumPy と最適化を行うための SciPy ライブラリをインポートする． 

import numpy as np
import pandas as pd
import scipy.optimize as optimize # 最適化
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize   # contour の色指定

# 高解像度ディスプレイ用
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')

CSVファイルから散布図を描く

CSVファイルを読み込んで，NumPy配列に変換し，これから散布図を描いてみる． 

# CSVファイルを読み込む
url = "https://github.com/rinsaka/sample-data-sets/blob/master/lr.csv?raw=true"
df = pd.read_csv(url)

# NumPy 配列に変換する
x_data = df.loc[:, 'x'].values
y_data = df.loc[:, 'y'].values

# 散布図を作成する
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x_data, y_data)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(0,10)
ax.set_ylim(0,10)
# plt.savefig('lr-data.png', dpi=300, facecolor='white')
plt.show()

lr-data.png

目次に戻る

回帰式を定義する

線形回帰を行いたいので，その回帰式を定義する．このページでは最も単純な一次関数を利用する．

\begin{equation} y = a + bx \end{equation}

上の式で \(a\) は切片，\(b\) は直線の傾きである．ここでは，x の値（またはリスト）とパラメータ w = [a, b] を与えれば，y の値（またはリスト）を返す関数を次のように定義する． 

関数の定義"""
関数の定義
"""
def my_func(w, x):
    y = w[0] + w[1] * x
    return y

定義した関数 my_func に \(a = 1\)，\(b = 2\)，\(x = 3\) を代入して \(y\) の値（つまり，\(1 + 2 \times 3\) ）を求めてみる． 

my_func([1, 2], 3)

7

定義した関数の引数 x に NumPy 配列を与えれば，計算結果が NumPy 配列で得られることに注意する．（つまり，繰り返して計算する必要がない．） 

my_func([1, 2], np.array([0, 1, 2]))

array([1, 3, 5])

回帰直線のグラフを描きたいので，まず x の値を格納する NumPy 配列 x_plotを準備する．このモデルの回帰式は直線であるので，始点と終点の2点のデータがあれば十分であるが，次のページ以降では曲線となる回帰式も描きたいので，等間隔で100点のデータを準備しておく． 

x のNumPy配列 x_plot を準備するx_plot = np.linspace(0, 10, 100)
print(x_plot)

[ 0.          0.1010101   0.2020202   0.3030303   0.4040404   0.50505051
  0.60606061  0.70707071  0.80808081  0.90909091  1.01010101  1.11111111
  1.21212121  1.31313131  1.41414141  1.51515152  1.61616162  1.71717172
  1.81818182  1.91919192  2.02020202  2.12121212  2.22222222  2.32323232
  2.42424242  2.52525253  2.62626263  2.72727273  2.82828283  2.92929293
  3.03030303  3.13131313  3.23232323  3.33333333  3.43434343  3.53535354
  3.63636364  3.73737374  3.83838384  3.93939394  4.04040404  4.14141414
  4.24242424  4.34343434  4.44444444  4.54545455  4.64646465  4.74747475
  4.84848485  4.94949495  5.05050505  5.15151515  5.25252525  5.35353535
  5.45454545  5.55555556  5.65656566  5.75757576  5.85858586  5.95959596
  6.06060606  6.16161616  6.26262626  6.36363636  6.46464646  6.56565657
  6.66666667  6.76767677  6.86868687  6.96969697  7.07070707  7.17171717
  7.27272727  7.37373737  7.47474747  7.57575758  7.67676768  7.77777778
  7.87878788  7.97979798  8.08080808  8.18181818  8.28282828  8.38383838
  8.48484848  8.58585859  8.68686869  8.78787879  8.88888889  8.98989899
  9.09090909  9.19191919  9.29292929  9.39393939  9.49494949  9.5959596
  9.6969697   9.7979798   9.8989899  10.        ]

いま準備した NumPy 配列 x_plot に対応する y の値（NumPy配列）を計算する．まずは，繰り返しを用いた C言語風の記述方法を使ってみる．

C言語風の書き方# 傾きと切片のリストを準備
w = np.array([4, 0.25])
# y_pred のリストを準備（ゼロで初期化）
y_pred = np.zeros(len(x_plot))
# 配列x_plotから一つずつ取り出して，y_predを計算する
for i, x in enumerate(x_plot):
    y_pred[i] = my_func(w, x)
print(y_pred)

[4.         4.02525253 4.05050505 4.07575758 4.1010101  4.12626263
 4.15151515 4.17676768 4.2020202  4.22727273 4.25252525 4.27777778
 4.3030303  4.32828283 4.35353535 4.37878788 4.4040404  4.42929293
 4.45454545 4.47979798 4.50505051 4.53030303 4.55555556 4.58080808
 4.60606061 4.63131313 4.65656566 4.68181818 4.70707071 4.73232323
 4.75757576 4.78282828 4.80808081 4.83333333 4.85858586 4.88383838
 4.90909091 4.93434343 4.95959596 4.98484848 5.01010101 5.03535354
 5.06060606 5.08585859 5.11111111 5.13636364 5.16161616 5.18686869
 5.21212121 5.23737374 5.26262626 5.28787879 5.31313131 5.33838384
 5.36363636 5.38888889 5.41414141 5.43939394 5.46464646 5.48989899
 5.51515152 5.54040404 5.56565657 5.59090909 5.61616162 5.64141414
 5.66666667 5.69191919 5.71717172 5.74242424 5.76767677 5.79292929
 5.81818182 5.84343434 5.86868687 5.89393939 5.91919192 5.94444444
 5.96969697 5.99494949 6.02020202 6.04545455 6.07070707 6.0959596
 6.12121212 6.14646465 6.17171717 6.1969697  6.22222222 6.24747475
 6.27272727 6.2979798  6.32323232 6.34848485 6.37373737 6.3989899
 6.42424242 6.44949495 6.47474747 6.5       ]

しかしながら，定義した関数 my_func に NumPy 配列を与えれば，すべての結果を計算できるので，Python では次のような書き方ができる．まだ試していませんが，この書き方の方が実行速度の面でも有利なはずです． 

切片 4, 傾き 0.25 として，回帰直線のデータを作成する# 傾きと切片のリストを準備
w = np.array([4, 0.25])
# NumPy配列では一気に計算できる
y_pred = my_func(w, x_plot)
print(y_pred)

[4.         4.02525253 4.05050505 4.07575758 4.1010101  4.12626263
 4.15151515 4.17676768 4.2020202  4.22727273 4.25252525 4.27777778
 4.3030303  4.32828283 4.35353535 4.37878788 4.4040404  4.42929293
 4.45454545 4.47979798 4.50505051 4.53030303 4.55555556 4.58080808
 4.60606061 4.63131313 4.65656566 4.68181818 4.70707071 4.73232323
 4.75757576 4.78282828 4.80808081 4.83333333 4.85858586 4.88383838
 4.90909091 4.93434343 4.95959596 4.98484848 5.01010101 5.03535354
 5.06060606 5.08585859 5.11111111 5.13636364 5.16161616 5.18686869
 5.21212121 5.23737374 5.26262626 5.28787879 5.31313131 5.33838384
 5.36363636 5.38888889 5.41414141 5.43939394 5.46464646 5.48989899
 5.51515152 5.54040404 5.56565657 5.59090909 5.61616162 5.64141414
 5.66666667 5.69191919 5.71717172 5.74242424 5.76767677 5.79292929
 5.81818182 5.84343434 5.86868687 5.89393939 5.91919192 5.94444444
 5.96969697 5.99494949 6.02020202 6.04545455 6.07070707 6.0959596
 6.12121212 6.14646465 6.17171717 6.1969697  6.22222222 6.24747475
 6.27272727 6.2979798  6.32323232 6.34848485 6.37373737 6.3989899
 6.42424242 6.44949495 6.47474747 6.5       ]

回帰直線のデータが作成できたので，散布図に重ねてみよう． 

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x_data, y_data)
ax.plot(x_plot, y_pred)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(0,10)
ax.set_ylim(0,10)
# plt.savefig('lr-model1-1.png', dpi=300, facecolor='white')
plt.show()

  lr-model1-1.png

目次に戻る

誤差を計算する

あるパラメータ w = [a, b] を与えて計算した回帰直線と一つのデータとの2乗誤差を計算してみよう．まずは2乗誤差を求める関数 get_square_error を定義する． 

def get_square_error(w, x, y):
    y_pred = my_func(w, x)
    error = (y - y_pred)**2
    return error

上と同じパラメータの組み合わせ（w = [4, 2.5]）のときに，先頭のデータ x = 0.1, y = 4.0 に対する予測値は 4.025 となり，その2乗誤差は 0.0006250000000000178 となることがわかる． 

w = np.array([4, 0.25])
print(x_data[0], y_data[0], my_func(w, x_data[0]))
print(get_square_error(w, x_data[0], y_data[0]))

0.1 4.0 4.025
0.0006250000000000178

次のデータ x = 0.2, y = 4.8 に対する予測値は 4.05 となるので，その2乗誤差は 0.5625 である． 

print(x_data[1], y_data[1], my_func(w, x_data[1]))
print(get_square_error(w, x_data[1], y_data[1]))

0.2 4.8 4.05
0.5625

さらに別のデータ x = 7.2, y = 1.0 に対する予測値は 5.8 となるので，その2乗誤差は大きくなり 23.04 である． 

print(x_data[34], y_data[34], my_func(w, x_data[34]))
print(get_square_error(w, x_data[34], y_data[34]))

7.2 1.0 5.8
23.04

目次に戻る

残差2乗和を計算する

上の手順では，1つのデータに対する2乗誤差を計算した．次は，すべてのデータについての2乗誤差の総和，つまり，残差2乗和 (rss: residual sum-of-squares) を計算する関数を定義しよう．Python では NumPy を使うと非常に短いコードで計算が可能になる．つまり for や while 等の繰り返し処理を必要としない． 

def get_rss(w, x, y):
    y_pred = my_func(w, x)
    error = (y - y_pred)**2
    return np.sum(error)

上と同じパラメータの組み合わせ（w = [4, 2.5]）のときの残差2乗和を計算すると rss = 331.03 となった． 

w = np.array([4, 0.25])
rss = get_rss(w, x_data, y_data)
print(rss)

331.031875

別のパラメータの組み合わせ（w = [5, -0.25]）に対しては rss = 261.04 となり，残差2乗和が減少した． 

w = np.array([5, -0.25])
rss = get_rss(w, x_data, y_data)
print(rss)

261.041875

上の2つの組み合わせでの回帰直線と散布図を重ね合わせると次のようになる． 

y_pred1 = my_func([4, 0.25], x_plot)
y_pred2 = my_func([5, -0.25], x_plot)
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x_data, y_data, label="data")
ax.plot(x_plot, y_pred1, label="y = 4 + 0.25x")
ax.plot(x_plot, y_pred2, label="y = 5 - 0.25x")
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(0,10)
ax.set_ylim(0,10)
ax.legend()
# plt.savefig('lr-model1-2.png', dpi=300, facecolor='white')
plt.show()

  lr-1-1

つまり，問題は残差2乗和 (rss) が最小になるようなパラメータ（切片，傾き）の組み合わせを求めることである．単回帰分析では連立方程式を解析的に解くことによって最適なパラメータを求めることが可能ですが，次ページ以降のモデルに対応するために，次のステップで最適化のイメージを掴んだ後に，さらに次のステップでは代表的な非線形最適化手法（ネルダーミード法，BFGS法）を使って求めてみます．

目次に戻る

最適化のイメージを掴む

上の例では，\(y = 4 + 0.25x\) と \(y = 5 -0.25x\) での残差2乗和を比較した．問題は \(y = a + bx\) について，残差2乗和が最も小さくなるような \(a\) と \(b\) を求めることである．

最適化のイメージを掴むために，まずは傾きを \(b = 0.25\) に固定した状態で，様々な \(a\) について rss を求めたものが下のグラフである．もしも \(b=0.25\) ならば，図から \(a=3.4\) あたりで rss が最小になることがわかる．これは上の図で傾きが 0.25 の直線を上下に平行移動したときに，切片が 3.4 になる付近で rss が最小になることを意味する．つまり，上の図で青の直線は（傾きを変化させないのであれば）もう少し下に移動したほうが良いということである． 

a_list = np.linspace(1,8,500)  # 傾き a を変化させる
b = 0.25 # 切片 b は固定
rss = np.zeros(a_list.shape[0]) # ゼロで初期化
for i, a in enumerate(a_list):
    rss[i] = get_rss([a, b], x_data, y_data)

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(a_list, rss, label="b = 0.25")
ax.set_xlabel('a')
ax.set_ylabel('rss')
ax.set_xlim(1,8)
ax.set_ylim(200,500)
ax.legend()
# plt.savefig('lr-rss-1.png', dpi=300, facecolor='white')
plt.show()

  lr-rss-1

さらに，4通りの傾きについて同じ様に rss を求めたものが下のグラフである．この図から \(b = -0.25\) 付近で rss が最も小さくなりそうであることが伺える．このときの切片は \(a = 5.6\) 程度になることもわかる． 

a_list = np.linspace(1,8,500)  # 傾き a を変化させる
b_list = np.array([-0.50, -0.25, 0.00, 0.25]) # 切片 b も変化させる
rss = np.zeros((a_list.shape[0], b_list.shape[0])) # ゼロで初期化
for j, b in enumerate(b_list):
    for i, a in enumerate(a_list):
        rss[i, j] = get_rss([a, b], x_data, y_data)

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
for j, b in enumerate(b_list):
    ax.plot(a_list, rss[:,j], label=f"b = {b}")
ax.set_xlabel('a')
ax.set_ylabel('rss')
ax.set_xlim(1,8)
ax.set_ylim(200,500)
ax.legend()
# plt.savefig('lr-rss-2.png', dpi=300, facecolor='white')
plt.show()

  lr-rss-2

今度は，切片を \(a = 6\) に固定したときに，最適な傾きがどの程度になるかを考える．下の図から傾きが \(b = -0.3\) 付近で rss が最小化される． 

a = 6
b_list = np.linspace(-0.8, 0.5, 500)
rss = np.zeros(b_list.shape[0]) # ゼロで初期化
for j, b in enumerate(b_list):
    rss[j] = get_rss([a, b], x_data, y_data)

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.plot(b_list, rss, label="a = 6")
ax.set_xlabel('b')
ax.set_ylabel('rss')
ax.set_xlim(-.8, .5)
ax.set_ylim(200,500)
ax.legend()
# plt.savefig('lr-rss-3.png', dpi=300, facecolor='white')
plt.show()

  lr-rss-3

さらに，様々な切片 \(a\) について同様の計算を行うと，次のようなグラフが得られる．やはりこのグラフからも \(a = 6\)，\(b = -0.3\) の周辺に rss を最小化する最適解がありそうであることが伺える． 

a_list = np.array([2, 4, 6, 8])
b_list = np.linspace(-0.8, 0.5, 500)
rss = np.zeros((b_list.shape[0], a_list.shape[0])) # ゼロで初期化
for i, a in enumerate(a_list):
    for j, b in enumerate(b_list):
        rss[j, i] = get_rss([a, b], x_data, y_data)

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
for j, a in enumerate(a_list):
    ax.plot(b_list, rss[:,j], label=f"a = {a}")
ax.set_xlabel('b')
ax.set_ylabel('rss')
ax.set_xlim(-.8, .5)
ax.set_ylim(200,500)
ax.legend()
# plt.savefig('lr-rss-4.png', dpi=300, facecolor='white')
plt.show()

  lr-rss-4

最後に，横軸に切片 \(a\)，縦軸に傾き \(b\) をとり，rss の大小を色の変化によってプロットしたものが下のグラフである．やはりこの図からも \(a = 5.5\)，\(b = -0.25\) の付近に最適解が存在しそうなことが読み取れる．このように，未知パラメータが1個や2個の場合には，グラフを描くことで最適解を探索するイメージを掴むことができる．未知パラメータが3個以上になると図では表現できず，最適解を見つけることは難しくなる．しかしながら，ネルダーミード法や準ニュートン法（BFGS法）などのよく知られた最適化アルゴリズムを用いれば未知パラメータの数が多くなっても最適解を比較的簡単に発見可能である．特にPythonではSciPy ライブラリなどこのような最適化が可能なライブラリが複数存在し，無料で利用できます．（さらに，有料の商用ソルバーを使えばより大規模な問題を高速・高精度に解くこともできる（らしい）．） 

aa, bb = np.meshgrid(np.linspace(4.0, 7.0, 300), np.linspace(-0.5, 0.0, 300))
rss = np.zeros(aa.ravel().shape[0])
for i in range(aa.ravel().shape[0]):
    rss[i] = get_rss([aa.ravel()[i], bb.ravel()[i]], x_data, y_data)
rss = rss.reshape(aa.shape)

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.contourf(aa, bb, rss)
mappable0 = ax.pcolormesh(aa, bb, rss, cmap='gist_rainbow', norm=Normalize(vmin=246, vmax=255))
fig.colorbar(mappable0, ax=ax, orientation="vertical")
ax.set_xlabel('a')
ax.set_ylabel('b')
ax.set_xlim(4.5,6.5)
ax.set_ylim(-0.5,0.0)
# plt.savefig('lr-rss-5.png', dpi=300, facecolor='white')
plt.show()

  lr-rss-5

目次に戻る

最適化

ここでは代表的な非線形最適化手法（ネルダーミード法，BFGS法）を使ってパラメータの最適化を行います．Python の SciPy ライブラリには強力な最適化関数が実装されているのでこれを用いれば比較的簡単に最適化が行なえます．なおこのページでは最適化のパラメータ数はわずか2個，次のページ以降の他のモデルでも5個程度ですが，ディープラーニングでは数千万や数億を超えるパラメータ数の最適化問題が解かれています．ここで説明している単純な文字認識でもおよそ3万個のパラメータです．

まず，SciPyライブラリがプログラムの先頭でインポートされていることを確認します．インポートされていなければ次の行を追加して実行してください． 

import scipy.optimize as optimize # 最適化

SciPy の optimize.minimze 関数を使ってネルダーミード法による最適化を行います．optimize.minimize 関数の引数には

  1. 最小化したい目的関数
  2. パラメータの初期値
  3. 目的関数に与える引数
  4. 最適化メソッド（使いたいアルゴリズム）

を指定します．今回は残差2乗和を最小化したいので，第1引数には rss を指定します．第2引数には探索したいパラメータに初期値を与える必要があります．初期値が異なれば別の局所解が得られたり，最適解を得られなかったりするので，問題の種類によっては試行錯誤的に適切な初期値を選ぶ必要があります．しかしながら，この問題ではそれほど気にする必要はないと思われます．目的関数である rss には最適化したいパラメータ w 以外に xとy を与える必要があるので，第3引数には args=(x_data,y_data) のようにタプル形式で指定します．第4引数には利用したい最適化アルゴリズムを指定します． 

w = np.array([0.1, 0.1]) # 初期値を設定
results_nm = optimize.minimize(get_rss, w, args=(x_data, y_data), method='Nelder-Mead')
print(results_nm)

 final_simplex: (array([[ 5.59136809, -0.24365817],
       [ 5.59141949, -0.24366314],
       [ 5.59145674, -0.24367417]]), array([244.84285679, 244.8428568 , 244.84285681]))
           fun: 244.84285679273873
       message: 'Optimization terminated successfully.'
          nfev: 121
           nit: 62
        status: 0
       success: True
             x: array([ 5.59136809, -0.24365817])

上のコードを実行した結果，目的関数値 (fun) が 244.84 となり，最適化が成功した (message: 'Optimization terminated successfully' および success: True) ことがわかります．最適解となるパラメータは x: [ 5.59136809, -0.24365817] になります．

最適化を行った結果から，目的関数値やパラメータは次のように取得できる． 

目的関数値print(results_nm["fun"])

244.84285679273873

最適解print(results_nm["x"])

[ 5.59136809 -0.24365817]

切片の最適解print(results_nm["x"][0])

5.591368093953666

傾きの最適解print(results_nm["x"][1])

-0.24365816990138328

次に，ブロイデン・フレッチャー・ゴールドファーブ・シャンノ法（BFGS法）によって同じ問題を最適化してみると，ほぼ同じ結果が得られました． 

w = np.array([0.1, 0.1]) # 初期値を設定
results_bfgs = optimize.minimize(get_rss, w,  args=(x_data, y_data), method='BFGS')
print(results_bfgs)

      fun: 244.8428567809353
 hess_inv: array([[ 0.05204723, -0.00864194],
       [-0.00864194,  0.00186052]])
      jac: array([1.90734863e-06, 3.81469727e-06])
  message: 'Optimization terminated successfully.'
     nfev: 21
      nit: 6
     njev: 7
   status: 0
  success: True
        x: array([ 5.591403  , -0.24366441])

ネルダーミード法によって得られた結果と散布図を重ねてみよう．このようなデータに直線を当てはめることには予想通り無理がありそうです． 

x_plot = np.linspace(0, 10, 100)
y_pred = my_func(results_nm["x"], x_plot)

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x_data, y_data, label="data")
ax.plot(x_plot, y_pred, label='model 1')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(0,10)
ax.set_ylim(0,10)
ax.legend()
# plt.savefig('lr-model1-3.png', dpi=300, facecolor='white')
plt.show()

  lr-model1-3.png

目次に戻る

StatsModelsとの比較

いま残差2乗和の最小化によって得られた結果と，ここで説明した StatsModels による単回帰分析と同じ方法で得られる結果を比較しておこう．どちらの方法でもほぼ同じ結果（切片と傾き）が得られていることが確認できました． 

import statsmodels.api as sm
x = sm.add_constant(x_data)
model = sm.OLS(y_data, x)
results = model.fit()
print(results.params)

[ 5.59140282 -0.24366437]

当てはまりの良さを R-squared や Adj. R-squared の値から確認する．これらの値が非常に小さいことから，やはりこのデータに対して単回帰分析による直線の当てはめは適切な方法でないことがわかる． 

print(results.summary2())

                 Results: Ordinary least squares
=================================================================
Model:              OLS              Adj. R-squared:     0.038
Dependent Variable: y                AIC:                197.2346
Date:               2021-07-17 16:39 BIC:                200.7099
No. Observations:   42               Log-Likelihood:     -96.617
Df Model:           1                F-statistic:        2.607
Df Residuals:       40               Prob (F-statistic): 0.114
R-squared:          0.061            Scale:              6.1211
-------------------------------------------------------------------
            Coef.    Std.Err.      t      P>|t|     [0.025   0.975]
-------------------------------------------------------------------
const       5.5914     0.7982    7.0049   0.0000    3.9782   7.2047
x1         -0.2437     0.1509   -1.6147   0.1142   -0.5487   0.0613
-----------------------------------------------------------------
Omnibus:              7.545        Durbin-Watson:           0.377
Prob(Omnibus):        0.023        Jarque-Bera (JB):        4.056
Skew:                 0.550        Prob(JB):                0.132
Kurtosis:             1.949        Condition No.:           11
=================================================================

目次に戻る

コード全体

不要な部分を削除した後のコード全体は次のような感じになる． 

import numpy as np
import pandas as pd
import scipy.optimize as optimize # 最適化
import matplotlib.pyplot as plt

# 高解像度ディスプレイ用
from IPython.display import set_matplotlib_formats
# from matplotlib_inline.backend_inline import set_matplotlib_formats # バージョンによってはこちらを有効に
set_matplotlib_formats('retina')

"""
関数の定義
"""
def my_func(w, x):
    y = w[0] + w[1] * x
    return y

"""
残差2乗和を求める関数
Residual Sum-of-Squares
"""
def get_rss(w, x, y):
    y_pred = my_func(w, x)
    error = (y - y_pred)**2
    return np.sum(error)

# CSV ファイルを読み込む
url = "https://github.com/rinsaka/sample-data-sets/blob/master/lr.csv?raw=true"
df = pd.read_csv(url)

# NumPy 配列に変換する
x_data = df.loc[:, 'x'].values
y_data = df.loc[:, 'y'].values

# ネルダーミード法による最適化を行う
w = np.array([0.1, 0.1]) # 初期値を設定
results_nm = optimize.minimize(get_rss, w, args=(x_data, y_data), method='Nelder-Mead')
print(results_nm)

# 最適解を使って回帰直線のデータを作成する
x_plot = np.linspace(0, 10, 100)
y_pred = my_func(results_nm["x"], x_plot)

# グラフを描く
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.scatter(x_data, y_data, label="data")
ax.plot(x_plot, y_pred, label='model 1')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(0,10)
ax.set_ylim(0,10)
ax.legend()
# plt.savefig('lr-model1-3.png', dpi=300, facecolor='white')
plt.show()

 final_simplex: (array([[ 5.59136809, -0.24365817],
       [ 5.59141949, -0.24366314],
       [ 5.59145674, -0.24367417]]), array([244.84285679, 244.8428568 , 244.84285681]))
           fun: 244.84285679273873
       message: 'Optimization terminated successfully.'
          nfev: 121
           nit: 62
        status: 0
       success: True
             x: array([ 5.59136809, -0.24365817])
lr-model1-3.png

次のページからは曲線を当てはめるいくつかの方法を試していこう．

目次に戻る

« 戻る 次へ »