Степень регрессии в степенном законе между curve_fit, python и excel - PullRequest
1 голос
/ 18 февраля 2020

Извините, я новичок в python и стека потока. Поэтому я не могу опубликовать изображение.

Я хочу сделать степенную регрессию с помощью функции curve_fit в python. Но результат довольно странный для меня. Я использую Excel, чтобы проверить это дальше. Это выглядит большой разницей между этими двумя. Черная линия - результат кривой_fit, а параметр красной линии - Excel. Может ли кто-нибудь дать мне знать разницу? Спасибо!

x=[164000,400,13000,700,57000,108,12000]
y=[0.011970,0.000098,0.066100,0.004300,0.042600,0.000061,0.002858 ]

def f(x,a,b):
    return a*x**b

popt,pocv=curve_fit(f,x,y)

ax.set_xscale("log")
ax.set_yscale("log")
ax.set_ylim(0.00001,0.1)
ax.set_xlim(10,1000000)

ax.scatter(x,y)

px=np.linspace(10,1000000,1000)

#parameter form curve_fit
py=a*px**b
[enter image description here][1]
#parameter from excel
pyy=3E-6*px**0.8305

ax.loglog(px,pyy,color="red")
ax.loglog(px,py,color="k")

Ответы [ 2 ]

0 голосов
/ 19 февраля 2020

Тот факт, что вы строите данные в пространстве журнала, должен дать вам хороший совет для размещения в пространстве журнала . То есть подойдет np.log(a*x**b) к np.log(y). Модификация вашего сценария, которая на самом деле выполняется и получает хорошее соответствие, будет выглядеть следующим образом:

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

x=[164000,400,13000,700,57000,108,12000]
y=[0.011970,0.000098,0.066100,0.004300,0.042600,0.000061,0.002858 ]

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

popt,pcov=curve_fit(f, x, np.log(y), [1.e-6, 0.9])

ax = plt.gca()

ax.set_xscale("log")
ax.set_yscale("log")
ax.set_ylim(0.00001,0.1)
ax.set_xlim(10,1000000)

ax.scatter(x,y)

px = np.linspace(10,1000000,1000)
a, b = popt
print("Parameters: a=%g,  b=%g" % (a, b))

#parameter form curve_fit
py=a*px**b

#parameter from excel
pyy=3e-6*px**0.8305

ax.loglog(px,pyy, color="red")
ax.loglog(px,py,  color="k")
plt.show()

Всегда обязательно указывайте начальные значения для параметров и обязательно распечатывайте результаты. В качестве примера, при выполнении этого будет напечатано Parameters: a=2.78612e-06, b=0.829763 и показано, что обе прогнозируемые линии находятся почти друг на друге.

Для лучшего подбора кривой вы можете найти lmfit (https://lmfit.github.io/lmfit-py/) полезно (да, я ведущий автор и предвзятый). С lmfit ваша подгонка может быть:

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

x=[164000,400,13000,700,57000,108,12000]
y=[0.011970,0.000098,0.066100,0.004300,0.042600,0.000061,0.002858 ]

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

model = Model(f)
params = model.make_params(a=1.e-6, b=0.9)
result = model.fit(np.log(y), params, x=x)

print(result.fit_report())

px = np.linspace(10,1000000,1000)
plt.scatter(x,y)
plt.loglog(px, np.exp(result.eval(x=px)),  color="k")
plt.show()

Обратите внимание, что с lmfit параметры имеют имя , используя имена в вашей функции модели f(). Это распечатает отчет о подгонке, который включает предполагаемые неопределенности:

[[Model]]
    Model(f)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 16
    # data points      = 7
    # variables        = 2
    chi-square         = 14.7591170
    reduced chi-square = 2.95182340
    Akaike info crit   = 9.22165592
    Bayesian info crit = 9.11347621
[[Variables]]
    a:  2.7861e-06 +/- 6.3053e-06 (226.31%) (init = 1e-06)
    b:  0.82976271 +/- 0.25700150 (30.97%) (init = 0.9)
[[Correlations]] (unreported correlations are < 0.100)
    C(a, b) = -0.958

и даст график enter image description here

0 голосов
/ 18 февраля 2020

Бревно с least_squares проблематично c, потому что различия между точками и выровненной линией увеличиваются с увеличением значения x. Вот почему подходит не так, как вы ожидали. Чтобы исправить это, вы можете добавить границы параметров, как показано ниже:

popt,pcov=curve_fit(f,x,y, bounds=(0, [3*10**(-6), 0.9]))

Суммарная ошибка параметров, которую вы можете получить из pcov.

error = np.sum(np.sqrt(np.diag(pcov)))

Вывод:

enter image description here

или изменить метод подгонки.

...