python：対数近似の求め方

Question

python初心者です．
pythonで作図した時系列データに対し，最小二乗法を用いた曲線近似（対数近似　：　y=a+b*INx）を求めたいと思っています．

scipy.optimize.curve_fit()を使おうと思ったのですが，私の力が足りずできませんでした.....

以下に近似線を引く前のプログラムです．

N1 = np.genfromtxt
(r_id_2020+'/' + date + '/' + 'sokudo'+ date+'/' + 'sokudo_001' +'_'+date+ kakutyo ,
delimiter= "," , skip_header=2, dtype='float')

T = N1[:,0]
x = N1[:,1]
y = N1[:,2]
Vx = N1[:,3]
Vy = N1[:,4]

plt.plot(x,Vx , c="red" ,linewidth = 3.0 ,linestyle = "--" ,label = "experiment")
plt.tight_layout()
plt.grid(True)
plt.legend(loc = 4,fontsize =18)
plt.xlim([0,100])
plt.ylim([0,200])
plt.show()

以下に示すのが読み取った値です．

横軸：x=[0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00
0.00000000e+00 1.90961041e-02 5.24035548e-02 2.85298953e-01
7.21075890e-01 1.27872008e+00 1.93054469e+00 2.74415196e+00
3.68723589e+00 4.75847352e+00 5.87141905e+00 7.04011478e+00
8.36277713e+00 9.70280669e+00 1.11143465e+01 1.25954020e+01
1.41392145e+01 1.57012355e+01 1.73083533e+01 1.89795096e+01
2.06446459e+01 2.24057242e+01 2.41901302e+01 2.60303676e+01
2.78307003e+01 2.95976600e+01 3.13087960e+01 3.30842440e+01
3.47851258e+01 3.65105431e+01 3.82488748e+01 4.00159721e+01
4.17402677e+01 4.34410278e+01 4.52448555e+01 4.69817836e+01
4.87180675e+01 5.04680193e+01 5.22016781e+01 5.39695510e+01
5.56590214e+01 5.73945792e+01 5.92032241e+01 6.09687188e+01
6.27228611e+01 6.45467608e+01 6.63263339e+01 6.81291919e+01
6.99411467e+01 7.17997065e+01 7.36481266e+01 7.54604648e+01
7.72207331e+01 7.90777275e+01 8.08353332e+01 8.26670779e+01
8.46531627e+01]

縦軸：Vx=[ 0. 0.60085865 0. 0. 0.
1.99844704 13.97372388 26.14661623 33.45865115 39.10947701
48.81643602 56.58503559 64.27425794 66.77673165 70.12174381
79.35974097 80.40177396 84.69239012 88.86332873 92.62875202
93.72125846 96.42706474 100.26938091 99.90817593 105.66469813
107.064359 110.41424644 108.019962 106.01758445 102.6681563
106.52688151 102.05290984 103.52503684 104.29990184 106.02583642
103.45773805 102.04560391 108.22966384 104.21568472 104.17703327
104.99710988 104.01952463 106.07237378 101.36822886 104.13346659
108.51869321 105.92968331 105.24853623 109.43398266 106.77438881
108.17147708 108.71728851 111.51358727 110.90520583 108.74029412
105.61609635 111.41966345 105.45634023 109.90468592 119.16508715
116.03951328]

以下の図のような近似線を引きたいです．

Answer 1

scipy.optimize.curve_fit() に最適化する関数、x の点一覧、y の点一覧を渡せば、最適化した結果の係数が返り値として得られます。

詳しい使い方はドキュメントを参照ください。

scipy.optimize.curve_fit — SciPy v1.5.0 Reference Guide

サンプルコード

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

np.random.seed(0)

a, b = 1.2, 2  # 真のパラメータ

def f(x, a, b):
    return a + b * np.log(x)


x = np.linspace(0.001, 10, 100)
y = f(x, a, b) + np.random.normal(0, 0.5, len(x))

a_pred, b_pred = curve_fit(f, x, y)[0]
y_pred = f(x, a_pred, b_pred)


# 描画する。
fig, ax = plt.subplots()

ax.plot(x, y, "gray", label=fr"${a:.2f} + {b:.2f} \log x + \mathcal{{N}}(0, 0.5^2)$")
ax.plot(x, y_pred, "r--", label=fr"${a_pred:.2f} + {b_pred:.2f} \log x$")
ax.legend()
ax.grid()

plt.show()