找到的散點線性擬合方法都是基于最小二乘法的（numpy.polyfit()，scipy.optimize())

以下是根據 GB/T 18459-2001中附錄 A2 提供的獨立線性度擬合方法，求得的最佳擬合直線

import mathdef find_line(x0, y0):'''根據散點求得端基直線k,b,并得到每點對端基直線的偏差dy:param x0: x坐標數組:param y0: y坐標數組:return: 每點的偏差dy'''# 求得端基直線的k,bk = (y0[-1] - y0[0]) / (x0[-1] - x0[0])d = y0[-1] - (k * x0[-1])# 根據端基直線的k,b,求得每點對端基直線的偏差dydy = [(y0[i] - ((k * x0[i]) + d)) for i in range(len(x0))]return dydef find_poly(x_lst, y_lst):'''找到凸多邊形,可以包含全部偏差點dy:param x_lst:x數組:param y_lst:偏差點dy數組:return:最佳凸多邊形各個點'''g_x = []g_y = []# 從起始點開始，向右（x增大方向），找到最大斜率點，再從最大斜率點開始，向右繼續尋找start_index = 0x_tmp = x_lst[1:]y_tmp = y_lst[1:]while len(x_tmp) > 0:k_lst = [((j - y_lst[start_index]) / (i - x_lst[start_index])) for i, j in zip(x_tmp, y_tmp)]start_index = k_lst.index(max(k_lst)) + 1 + start_indexg_x.append(x_lst[start_index])g_y.append(y_lst[start_index])x_tmp = x_lst[start_index + 1:]y_tmp = y_lst[start_index + 1:]# 從最末點開始，向左（x減小方向），也找到最大斜率點（有人說找最小斜率點，但是結果算出來不對），再從最大斜率點開始，向左繼續尋找start_index = len(x_lst) - 1x_tmp = x_lst[:start_index]y_tmp = y_lst[:start_index]while len(x_tmp) > 0:k_lst = [((j - y_lst[start_index]) / (i - x_lst[start_index])) for i, j in zip(x_tmp, y_tmp)]start_index = k_lst.index(max(k_lst))g_x.append(x_lst[start_index])g_y.append(y_lst[start_index])x_tmp = x_lst[:start_index]y_tmp = y_lst[:start_index]return g_x, g_ydef find_cross(p1, p2, p3):'''根據一個點p1和一條直線(p2和p3的連線)，求得該點對直線的鉛垂線（平行于縱軸坐標的直線）和直線的交點:param p1: 點:param p2: 線上一點:param p3: 線上另一點:return: 鉛垂線與線（p2-p3）的交點'''k = (p2[1] - p3[1]) / (p2[0] - p3[0])b = p2[1] - k * p2[0]p4 = [p1[0], (k * p1[0]) + b]return p4def find_perfect_line(a, b, x0, y0):'''根據凸多邊形的每個點，找到凸多邊形內最長的一根鉛垂線，與最長垂線相交的直線l1的斜率就是最佳直線的斜率過最長鉛垂線的中點作直線l2平行于直線l1，直線l2為最佳直線:param a:凸多邊形的x軸數組:param b:凸多邊形的y軸數組:param x0: 實際x軸數組:param y0: 實際y軸數組:return:最佳直線的斜率k, 截距b'''dic = {}point_lst = [[i, j] for i, j in zip(a, b)]for i in range(len(a)):p1 = (a[i], b[i]) # 凸多邊形的每個點坐標rp1 = (x0[x0.index(p1[0])], y0[x0.index(p1[0])])for j in range(len(a)):p2 = (a[j], b[j]) # 凸多邊形的每個點坐標rp2 = (x0[x0.index(p2[0])], y0[x0.index(p2[0])])if j == len(a) - 1:p3 = (a[0], b[0])else:p3 = (a[j + 1], b[j + 1])rp3 = (x0[x0.index(p3[0])], y0[x0.index(p3[0])])if p1 == p2 or p1 == p3:passelse:# print('====+++===', p1, p2, p3)# print('====+++===', rp1, rp2, rp3)# 鉛垂線長度s_length = abs((p2[1] * (p1[0] - p3[0])) + (p1[1] * (p3[0] - p2[0]) + (p3[1] * (p2[0] - p1[0]))) / (p3[0] - p2[0]))s = find_cross(p1, p2, p3)if isPointinPolygon(s, point_lst):dic[s_length] = [rp1, rp2, rp3]print('鉛垂線長度', s_length, s, dic[s_length])# print(dic)max_s_length = dic[max(dic)]rp1 = max_s_length[0]rp2 = max_s_length[1]rp3 = max_s_length[2]# 點和直線兩端點中點連線的斜率k = (((rp1[1] + rp2[1]) / 2) - ((rp1[1] + rp3[1]) / 2)) / (((rp1[0] + rp2[0]) / 2) - ((rp1[0] + rp3[0]) / 2))b = ((rp1[1] + rp2[1]) / 2) - (k * ((rp1[0] + rp2[0]) / 2))return k, bdef isPointinPolygon(point, rangelist): # [[0,0],[1,1],[0,1],[0,0]] [1,0.8]print(point)# 判斷是否在外包矩形內，如果不在，直接返回falselnglist = []latlist = []for i in range(len(rangelist) - 1):lnglist.append(rangelist[i][0])latlist.append(rangelist[i][1])# print(lnglist, latlist)maxlng = max(lnglist)minlng = min(lnglist)maxlat = max(latlist)minlat = min(latlist)# print(maxlng, minlng, maxlat, minlat)if (point[0] > maxlng or point[0] < minlng orpoint[1] > maxlat or point[1] < minlat):return Falsecount = 0point1 = rangelist[0]for i in range(1, len(rangelist)):point2 = rangelist[i]# 點與多邊形頂點重合if (point[0] == point1[0] and point[1] == point1[1]) or (point[0] == point2[0] and point[1] == point2[1]):print("在頂點上")return False# 判斷線段兩端點是否在射線兩側 不在肯定不相交 射線（-∞，lat）（lng,lat）if (point1[1] < point[1] and point2[1] >= point[1]) or (point1[1] >= point[1] and point2[1] < point[1]):# 求線段與射線交點 再和lat比較point12lng = point2[0] - (point2[1] - point[1]) * (point2[0] - point1[0]) / (point2[1] - point1[1])# print(point12lng)# 點在多邊形邊上if (point12lng == point[0]):print("點在多邊形邊上")return Trueif (point12lng < point[0]):count += 1point1 = point2print(count)if count % 2 == 0:return Trueelse:return Trueif __name__ == '__main__':x0 = [1.00, 2.00, 3.00, 4.00, 5.00, 6.00]y0 = [2.02, 4.00, 5.98, 7.90, 10.10, 12.05]dy = find_line(x0, y0)print(dy)a0, b0 = find_poly(x0, dy)print(a0, b0)print(find_perfect_line(a0, b0, x0, y0))

轉載于:https://www.cnblogs.com/sunqim16/p/9121011.html

總結

以上是生活随笔為你收集整理的独立线性度 最佳直线的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。