python 最小二乗法によるデータの多項式近似の計算エラー

Question

実現したいこと

多項式の係数を正しく導出する．

前提

pythonでcsvからデータを読み込み，
lu分解，前進代入，後退代入を用いて最小二乗法を行い，
得られた係数から曲線をプロットするプログラムを作りました．
しかし，検算のためにライブラリ(polyfit)を使用し求めたところ
自作したものから導出された係数と違う値が出ました．
1次式，2次式は正しくもとまりましたが，
3次式以降の係数が間違って導出されてしまいます．

発生している問題・エラーメッセージ

3次式近似の時の自作プログラムから導出された値
-62289409924.034035
93004892.05532321
-46288.47113927059
7.679190499441964

polyfit関数から求められた値
-3.12020048e+10
4.65818816e+07
-2.31807582e+04
3.84516090e+00    

該当のソースコード

import matplotlib.pyplot as plt
import numpy as np
import csv

def read(filename):
   csv_file = open(filename,'r')
   x = []
   y = []
   for row in csv.reader(csv_file):
      x.append(float(row[0]))
      y.append(float(row[1]))
   return x,y
      
def a_list(x,n):
   m = len(x) #データの個数
   a = [[float(0.0) for i in range(n + 1)] for j in range(n + 1)]
   for i in range(n + 1):
      for j in range(n + 1): #n
         sum_ = 0.0
         for k in range(m):
            sum_ += pow(x[k],i + j)
         a[i][j] = sum_
   return a

def b_list(x,y,n):
   n += 1
   b = [float(0.0) for i in range(n)]
   m = len(x)
   for i in range(n):
      sum_ = 0.0
      for k in range(m):
         sum_ += y[k] * pow(x[k],i)
      b[i] = sum_
   return b

#a行列をLU分解
def luBunkai(data,n):
   n += 1
   a = [[float(0.0) for i in range(n)] for j in range(n)]
   b = [[float(0.0) for i in range(n)] for j in range(n)]
   for i in range(n):
      a[i][i] = 1.0
   for j in range(n):
      for i in range(j + 1): #j+1
         sum_ = 0.0
         for k in range(i): #i
            sum_ += a[i][k] * b[k][j]
         b[i][j] = data[i][j] - sum_

      for i in range(j + 1,n): #j #j+1
         sum_ = 0.0
         for k in range(j): #j
            sum_ += a[i][k] * b[k][j]
         a[i][j] = (data[i][j] - sum_) / b[j][j] #j
   return a,b

#前進代入
def zensin(l,b,n):
   n += 1
   y = [float(0.0) for i in range(n)]
   y[0] = b[0] / l[0][0]
   for i in range(1,n):
      sum_ = 0.0
      for j in range(i): #i
         sum_ += l[i][j] * y[j]
      y[i] = (b[i] - sum_) / l[i][i]
   return y

#後退代入
def koutai(y,u,n):
   n += 1
   x = [float(0.0) for i in range(n)]
   x[n - 1] = y[n - 1] / u[n - 1][n - 1]
   for i in range(n - 2,-1,-1): #n-1
      sum_ = 0.0
      for j in range(i + 1,n): #n
         sum_ += u[i][j] * x[j]
      x[i] = (y[i] - sum_) / u[i][i]
   
   print("係数")
   for _ in x:
      print(_)

   return x

#得られた係数でプロットする関数
def fit_plot(c,x,y,n):
   n += 1
   fit = 0.0
   p = np.linspace(1996,2022,1000)
   for i in range(n):
      fit += c[i] * pow(p,i)
   plt.plot(p,fit,label = "n =" + str(n - 1))
   plt.xlabel("year")
   plt.ylabel("y")
   plt.legend()
   return 0

#1~nまでプロットする
#最終的に1~nとするためfor
for n in range(3,4):
   x,y = read("data.csv")
   a = a_list(x,n)
   b = b_list(x,y,n)
   l,u = luBunkai(a,n)
   lu_y = zensin(l,b,n)
   c = koutai(lu_y,u,n)
   fit_plot(c,x,y,n)
plt.xlim(1995,2022)
plt.ylim(3000,20000)
plt.scatter(x,y,label = "data")
plt.show()

#検算
for i in range(3,4):
   p4 = np.polyfit(x, y, i)
   print('ライブラリの値{0}'.format(p4))
   plt.plot(x, np.poly1d(np.polyfit(x, y, i))(x))
plt.xlim(1995,2022)
plt.ylim(3000,20000)
plt.scatter(x,y,label = "data")
plt.show()

"""data.datの内容
1996,6216.833333333333
1997,10050.819444444443
1998,12947.568287037036
1999,11282.380690586418
2000,12065.031724215534
2001,9488.919310351295
2002,11673.326609195943
2003,14021.860550766329
2004,16618.821712563862
2005,16300.151809380322
2006,17646.09598411503
2007,21248.091332009586
2008,13241.674277667466
2009,12615.556189805624
2010,16353.713015817135
2011,9480.726084651427
2012,5124.893840387619
2013,4430.157820032302
2014,4850.513151669358
2015,14092.876095972446
2016,18467.406341331036
2017,20329.700528444253
2018,17423.975044037023
2019,17081.74792033642
2020,14347.562326694702
2021,18582.130193891226
2022,11001.927516157602
"""

試したこと

プログラム中のa,b,l,u,y,xをトレースして，lu分解，前進代入，後退代入が間違っていないことを確認しました．
data.datを用いて最小二乗法を行う外部のツールで計算したところ，polyfit関数と同じ値になりました．

追記

2023/12/27

Answer 1

演算をdecimal型にしたところ，導出された値がライブラリに近い値となったため解決とします．（値に誤差はあります）
原因は精度の問題だったと思います．

Comments

[winterboum]

おめでとう。
その精度の問題が私の回答です。
もう解決したの( と codeをきちんと読み解く時間ないの）で簡単に言うと、
X 1996〜2022、 Y 4000〜20000 の世界を 座標移動して X -13〜13、Y -12000 〜 12000 の世界にして計算すると精度を得するのです。
例えば精度８桁の計算機だったとすると
2022を3乗すると 8000000000位 になりますが13だと2197。8000000000 の上の桁は意味がなくて最後の4桁が重要ですが、8000000000 のままだと精度が足りなくなって下４桁部分が吐き出されてしまいます。

Answer 2

気になることがあります。
データが 6216.833333333333 と16桁もあります。有効桁数そんなに無いでしょう？
通常の人が手にできる有効桁って4桁、精密な秤で5桁が良いところ。16桁というのは生の測定値ではなく計算結果ですよね？ 計算の元になってる生データの有効桁数から計算結果の有効桁数を求めて、その桁で丸めるのが望ましいです。

で、プログラム
分散を求めるのに Σ (X**2) - 1/n (ΣX)**2 してますね？
筆算ならこれですが、計算機使うなら分散の定義通りに Σ(( ｘ - xの平均) **2)でやる方が精度が出ます。
それが原因かどうかは試していただかないとなんとも、ですが。

Comments

[meg_]

重複投稿の様です。

[Saku]

データに関してはご指摘の通りで計算結果です．
11 y.append(round(float(row[1])))
とround関数で大まかに小数点以下を処理しましたが，
ライブラリで求めた値には近づきませんでした．
（追記　違いました．計算結果を丸めるということなので，丸めたものを計算することではなかったです）

プログラムの方ですが，分散を求めるというのは具体的に何行目のプログラムでしょうか？
a_listとb_listというのは，xとyの二乗誤差を最小にするために変微分して0とする
という式から持ってきたものですが（質問に貼っておきます），そのあたりでしょうか．
理解不足で申し訳ありません．