Поиск коэффициентов уравнения

Question

В уравнении вида y = alpha + k*x^n необходимо найти коэффициенты alpha, k, n. Показатель степени n от 0 до 1.

При использовании функции curve_fit из scipy получаются совсем некорректные значения.

Тестовые данные x = [1021.998 10.21998 5.10999] y = [39.8112 15.312 1.5312]. Должно было получится alpha = 215.4611 n = 0.32 k = 4.3831. если использовать curve_fit, то на выходе alpha = -77523.501; k=77518, 284; n = 0.

from scipy.optimize import curve_fit 

def calc_func(x, alpha, k, n): 
    return alpha + k * x ** n

coeff, pcov = curve_fit(calc_func, x[np.isfinite(y)], y[np.isfinite(y)], maxfev=10000)

Выяснила, что для такого типа функции такой метод не подходит, нужно решать численными методами минимизации.

Есть ли таких случаев метод решения средствами scipy? Может есть другой способ решения

Answer 1

У вас очень неудачный набор данных. Две точки находятся почти рядом, а третья точка - очень далеко. В результате вычислительная задача получается крайне неустойчивой и curve_fit не сходится.

Для решения можно использовать тот факт, что n меняется в небольшом диапазоне. Тогда вместо решения задачи с тремя параметрами можно итерировать n и вручную выбрать решение.

import numpy as np
import scipy.optimize as spo
import matplotlib.pyplot as plt

x = np.array([1021.998, 10.21998, 5.10999])
y = np.array([39.8112, 15.312, 1.5312])

def fn(x, a,k,n):
    return a + k*(x**n)

for n in np.linspace(0.01,0.3,20):
    p_1, p_cov_1 = spo.curve_fit(lambda x,a,k: fn(x,a,k, n), x, y )
    err = np.sqrt(np.diag(p_cov_1))
    print("n={0:.3f}, (a,k)={1}±{2}, относительная ошибка {3}".format(n, p_1, err, err/p_1))

Получается вот какой результат:

n=0.010, (a,k)=[-629.213  624.831]±[162.505 156.627], относительная ошибка [-0.258  0.251]
n=0.025, (a,k)=[-234.577  230.849]±[64.832 58.94 ], относительная ошибка [-0.276  0.255]
n=0.041, (a,k)=[-137.351  134.248]±[40.795 34.886], относительная ошибка [-0.297  0.26 ]
n=0.056, (a,k)=[-93.434  90.929]±[29.96  24.033], относительная ошибка [-0.321  0.264]
n=0.071, (a,k)=[-68.471  66.538]±[23.823 17.874], относительная ошибка [-0.348  0.269]
n=0.086, (a,k)=[-52.407  51.019]±[19.894 13.921], относительная ошибка [-0.38   0.273]
n=0.102, (a,k)=[-41.23   40.362]±[17.178 11.18 ], относительная ошибка [-0.417  0.277]
n=0.117, (a,k)=[-33.025  32.653]±[15.202  9.175], относительная ошибка [-0.46   0.281]
n=0.132, (a,k)=[-26.76   26.862]±[13.708  7.653], относительная ошибка [-0.512  0.285]
n=0.147, (a,k)=[-21.834  22.386]±[12.548  6.463], относительная ошибка [-0.575  0.289]
n=0.163, (a,k)=[-17.868  18.85 ]±[11.628  5.511], относительная ошибка [-0.651  0.292]
n=0.178, (a,k)=[-14.616  16.005]±[10.885  4.737], относительная ошибка [-0.745  0.296]
n=0.193, (a,k)=[-11.906  13.684]±[10.278  4.098], относительная ошибка [-0.863  0.299]
n=0.208, (a,k)=[-9.621 11.768]±[9.775 3.564], относительная ошибка [-1.016  0.303]
n=0.224, (a,k)=[-7.672 10.17 ]±[9.356 3.113], относительная ошибка [-1.219  0.306]
n=0.239, (a,k)=[-5.995  8.826]±[9.004 2.729], относительная ошибка [-1.502  0.309]
n=0.254, (a,k)=[-4.541  7.688]±[8.706 2.401], относительная ошибка [-1.917  0.312]
n=0.269, (a,k)=[-3.27   6.718]±[8.452 2.118], относительная ошибка [-2.585  0.315]
n=0.285, (a,k)=[-2.154  5.887]±[8.236 1.873], относительная ошибка [-3.823  0.318]
n=0.300, (a,k)=[-1.169  5.171]±[8.05  1.659], относительная ошибка [-6.889  0.321]

Лично мне нравится n=0.147, alpha = -21.834, k = 22.386. Абсолютная ошибка и относительная приемлемые. Конечно, не идеальные, но терпимые.

Видите на графике как сильно изгибается кривая? Это от того, что первые две точки находятся практически на одной вертикали, если сравнивать с правой точкой. Для более качественной аппроксимации желательно, чтобы вторая точка была где-то посередине между первой и третьей.

Так как первые две точки почти на одной вертикали, то n должно быть около нуля, чтобы кривая решения быстро росла.

Но в этом случае очень быстро нарастает абсолютная ошибка для alpha и k. Как ни крути, очень неудачные данные.

На хороших данных curve_fit работает вполне себе успешно.

a = -10; k = 3; n = 0.1

u = np.array([1,500,100])
v = fn(u, a,k,n)

p, p_cov = spo.curve_fit(fn, u,v)

Результат p = array([-10. , 3. , 0.1]), ковариационная матрица p_cov не определена, так как всего три точки для трёх параметров. Но главное в другом - все три параметра найдены абсолютно точно.

БЛОГ НА HUSL

Поиск коэффициентов уравнения

Ответы (1 шт):