折れ線回帰をpythonで実施したいです。

python

1import numpy as np 2import matplotlib.pyplot as plt 3from scipy.optimize import curve_fit 4 5data = [6 [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17],7 [-1.040122, -1.563332, -1.792177, -0.926176, -0.746776, -1.260266, -1.046665, -0.058491, 0.999898, 1.642013, 3.016936, 2.671379, 0.865827, 0.802594, -0.738069, -0.666671, -0.159904]8]9 10x = np.array(data[0], dtype=float)11y = np.array(data[1], dtype=float)12n = x.shape[0]13 14def func1(x, a1, b1):15 return a1 * x + b1 16 17def func2(x, a2, x0, a1, b1):18 return a2 * (x - x0) + (a1 * x0 + b1)19 20best_error = float('inf')21best_i = None22best_params1 = None23best_params2 = None24 25for i in range(2, n - 2):26 seg1_x = x[:i + 1]27 seg1_y = y[:i + 1]28 seg2_x = x[i:]29 seg2_y = y[i:]30 31 params1, _ = curve_fit(func1, seg1_x, seg1_y)32 a1, b1 = params1 33 34 params2, _ = curve_fit(func2, seg2_x, seg2_y, p0=[1, x[i], a1, b1])35 a2, x0, _, _ = params2 36 37 ans1 = func1(seg1_x, a1, b1)38 ans2 = func2(seg2_x, a2, x0, a1, b1)39 40 error = np.sum((ans1 - seg1_y) ** 2) + np.sum((ans2 - seg2_y) ** 2)41 42 if error < best_error:43 best_error = error 44 best_i = i 45 best_params1 = params1 46 best_params2 = params2 47 48seg1_x = x[:best_i + 1]49seg1_y = y[:best_i + 1]50seg2_x = x[best_i:]51seg2_y = y[best_i:]52 53a1, b1 = best_params1 54a2, x0, _, _ = best_params2 55 56ans1 = func1(seg1_x, a1, b1)57ans2 = func2(seg2_x, a2, x0, a1, b1)58 59plt.plot(x, y, '-o', color='red', label="data")60plt.plot(seg1_x, ans1, '--', color='blue')61plt.plot(seg2_x, ans2, '--', color='green')62plt.show()

コメントを投稿

0 コメント