Excelのグラフが描く三次近似式を表現したい!!

Question

解決したいこと

Excelのグラフで算出される三次近似式と同等の結果を算出したいです。これまでPythonのnumpyやR言語のライブラリで算出してみたのですが、データのコンディションの問題等でエラーで返されてしまい、上手く算出できませんでした。

Pythonの方では二次近似ではExcelのグラフとほぼ同等の結果を得ることができたのですが、三次近似式になると上手く算出できませんでした。

同じくPythonのnumpyにて下記のような行列計算を行い、三次近似式の最小二乗法を行ってみたのですが、上手く算出できず。整数部が同じ値なので、小数部のみで計算を行ってみても上手く算出できずでした。

該当するソースコード

from django.shortcuts import render
from django.http import HttpResponse
from django.db.models import Max
from io import TextIOWrapper, BytesIO
import csv, base64, math, logging
import pandas as pd, numpy as np, matplotlib.pyplot as plt, japanize_matplotlib
from sklearn.metrics import r2_score
from statistics import mean
import warnings
import scipy as sp
from scipy import optimize
import math
from scipy import interpolate
import decimal
# warnings.simplefilter('ignore', np.RankWarning)

"""
検証
"""
pd.options.display.float_format = '{:.30f}'.format
np.set_printoptions(suppress=True)

a = [34.786724,34.786701,34.786678,34.786656,34.786633,34.786611,34.786588,34.786566,34.786543,34.786521,34.786498,34.786476,34.786453,34.78643,34.786408,34.786387,34.786365,34.786344,34.786323,34.786306,34.786289,34.786272,34.786255,34.786238,34.786221,34.786207,34.786196,34.786185,34.786174,34.786163,34.786151,34.78614,34.786129,34.786123,34.786117,34.786111,34.786105]
b = [133.875015,133.874911,133.874807,133.874704,133.874600,133.874496,133.874392,133.874288,133.874184,133.874080,133.873976,133.873872,133.873768,133.873664,133.873560,133.873455,133.873351,133.873246,133.873142,133.873037,133.872931,133.872826,133.872720,133.872614,133.872509,133.872403,133.872296,133.872190,133.872083,133.871977,133.871870,133.871763,133.871657,133.871550,133.871443,133.871336,133.871229]

"""
二次近似(polyfit)
"""
print("\n>>>>>>>>>>polyfitによる回帰分析")
p_1 = np.polyfit(a, b, 2)
p_2 = np.polyfit(b, a, 2)
print("3次近似式係数(経度/緯度) : %s"%(p_1))
print("3次近似式係数(緯度/経度) : %s"%(p_2))

yfit_1 = np.polyval(p_1, a)
r2_1 = r2_score(b, yfit_1)
r2_1 = round(r2_1, 5)
print("決定係数(経度/緯度) : %s"%(r2_1))

yfit_2 = np.polyval(p_2, b)
r2_2 = r2_score(a, yfit_2)
r2_2 = round(r2_2, 5)
print("決定係数(緯度/経度) : %s"%(r2_2))

"""
三次近似(polyfit)
"""
print("\n>>>>>>>>>>polyfitによる回帰分析")
p_1 = np.polyfit(a, b, 3)
p_2 = np.polyfit(b, a, 3)
print("3次近似式係数(経度/緯度) : %s"%(p_1))
print("3次近似式係数(緯度/経度) : %s"%(p_2))

yfit_1 = np.polyval(p_1, a)
r2_1 = r2_score(b, yfit_1)
r2_1 = round(r2_1, 5)
print("決定係数(経度/緯度) : %s"%(r2_1))

yfit_2 = np.polyval(p_2, b)
r2_2 = r2_score(a, yfit_2)
r2_2 = round(r2_2, 5)
print("決定係数(緯度/経度) : %s"%(r2_2))

"""
三次近似(curve_fit)
"""
def func(x, a, b, c, d):
    return a * pow(x, 3) + b * pow(x, 2) + c * x + d

print("\n>>>>>>>>>>curve_fitによる回帰分析")
# 近似式の作成
popt_1, pcov_1 = optimize.curve_fit(func, a, b)
popt_2, pcov_2 = optimize.curve_fit(func, b, a)
print("3次近似式係数(経度/緯度) : %s"%(popt_1))
print("3次近似式係数(緯度/経度) : %s"%(popt_2))

df_a = pd.DataFrame(a)
df_b = pd.DataFrame(b)

r2_5 = r2_score(df_b[0], func(df_a[0], *popt_1))
r2_5 = round(r2_5, 5)
print("決定係数(経度/緯度) : %s"%(r2_5))

r2_6 = r2_score(df_a[0], func(df_b[0], *popt_2))
r2_6 = round(r2_6, 5)
print("決定係数(緯度/経度) : %s"%(r2_6))

上記の結果

>>>>>>>>>>polyfitによる回帰分析
3次近似式係数(経度/緯度) : [   -4910.5526179    341646.57037832 -5942291.96429974]
3次近似式係数(緯度/経度) : [    25.91131111  -6937.48621553 464394.87552564]
決定係数(経度/緯度) : 0.99421
決定係数(緯度/経度) : 0.9989

>>>>>>>>>>polyfitによる回帰分析
C:\Users\nishi\anaconda3\lib\runpy.py:87: RankWarning: Polyfit may be poorly conditioned
  exec(code, run_globals)
C:\Users\nishi\anaconda3\lib\runpy.py:87: RankWarning: Polyfit may be poorly conditioned
  exec(code, run_globals)
3次近似式係数(経度/緯度) : [     -70.58284125     2455.38573984    85413.15287949 -2971158.21083739]
3次近似式係数(緯度/経度) : [     0.09677439    -12.95517485  -1734.30636101 232206.14080975]
決定係数(経度/緯度) : 0.99421
決定係数(緯度/経度) : 0.9989

>>>>>>>>>>curve_fitによる回帰分析
3次近似式係数(経度/緯度) : [   0.081146     -3.83261409  -22.20142109 2128.17822341]
3次近似式係数(緯度/経度) : [   -0.00623885     1.24087537     3.36861239 -7686.44757283]
決定係数(経度/緯度) : 0.97607
決定係数(緯度/経度) : 0.97383

###Excelグラフの算出結果
下記がExcelの算出結果です。下のグラフの方が二次近似式で表現されており、値は"y = -4,910.5671997070 x2 + 341,647.5848742150 x - 5,942,309.6096078100"と、Pythonのpolyfitの二次近似で算出した値とほぼ同等の結果が得られますが、上のグラフの方の三次近似式で表現したグラフの近似式は同等の結果を得られません。Excelの三次近似式の値は”y = 15,771,219.9687500000 x3 - 1,645,876,935.7577200000 x2 + 57,254,308,654.6703000000 x - 663,892,459,167.0570000000”です。

Ptyhon(numpy)で発生している問題・エラー

また、polyfitで三次近似式を算出すると下記のようなエラーが表示されてしまいます。データのコンディションが不適切だとされてしまいます。

 RankWarning: Polyfit may be poorly conditioned
  exec(code, run_globals)

何とかExcelのグラフの三次近似式と同等の結果を得られるようにしたく、皆様のお知恵をお借りしたいです。お忙しいところ申し訳ありませんが何卒宜しくお願い致します。

Answer 1

from numpy.polynomial import Polynomial as P
p3 = P.fit(a, b, 3)
print(p3.convert().coef[::-1])
print(r2_score(b, p3(np.array(a))))

aa = np.linspace(np.min(a), np.max(a), 100)
plt.scatter(a, b)
plt.plot(aa, p3(aa), 'r-')
plt.show()

Comments

[NS78]

ご回答、誠にありがとうございます。

おぉ！
回答いただいたコードで上手くいきました！

自分なりに色々本も読んで調べたつもりだったのですが、このような関数がnumpyにあったとは知らず、長い間迷走しておりました。

jbpb0様、本当に助かりました！
ありがとうございました！