Решите для констант 4-параметра (уравнение Родбарда) с использованием Python

Я новичок в python, и я пытаюсь написать алгоритм для решения 4 неизвестных параметров в уравнении Родбарда, где мы связываем значение оттенков серого, измеренное с помощью ImageJ, с калибровочными дисками оптической плотности. Это уравнение является нелинейным и записывается как y = c*((x-a)/(d-x))^(1/b), где a, b, c и d неизвестны. У меня есть значения x и y для четырех точек (176,5, 0), (161,333, 0,1), (66,1667, 0,9) и (40,833, 2,5). Ниже я разместил свою попытку решить эти 4 неизвестных. Любая помощь, чтобы указать мне в правильном направлении, будет принята с благодарностью!

    import scipy.optimize as opt

    def f(a, b, c, d):
         0 == [c * ((176.5 - a)/(d - 176.5))**(1/b)]
         0.1 == [c * ((161.333 - a)/(d - 161.333))**(1/b)]
         0.9 == [c * ((66.1667 - a)/(d - 66.1667))**(1/b)]
         2.5 == [c * ((40.833 - a)/(d - 40.833))**(1/b)]
    return f

    opt.curve_fit(a, b, c, d)

    print a
    print b
    print c
    print d

person A.LeBrun    schedule 17.02.2016    source источник
comment
Я попытался отредактировать его, следуя примеру, но он все еще не работает для меня. Извините, я не смог заставить его выглядеть как код, но вот он: импортируйте numpy как np из scipy.optimize )/(d - x))**(1/b)] xdata = np.linspace(0, 1, 255) y = func(xdata, 0, 0,1, 0,9, 2,5) ydata = y + 0,2 * np. random.normal(size=len(xdata)) popt, pcov = curve_fit(func, xdata, ydata)   -  person A.LeBrun    schedule 18.02.2016


Ответы (2)


Если вы хотите использовать curve_fit, вы должны сделать следующее:

def f(x1, a1, b1, c1, d1):
    return c1 * (((x1 - a1)/(d1 - x1))**1/b1)

x_data = np.array([176.5, 161.333, 66.1667, 40.833])
y_data = np.array([0., 0.1, 0.9, 2.5])
p0 = np.array([168., -0.01, -7.4, 35000.])

popt, pcov = opt.curve_fit(f, x_data, y_data, p0, None, False, True, ftol = 0.00001)

p0 - это начальное предположение, и если вы его не сообщите, будет принят массив из всех единиц.

С предоставленными данными я пробовал разные параметры, но не смог найти решение.

Надеюсь, это поможет. Удачи!

person Francesc Torradeflot    schedule 17.02.2016
comment
Я пробовал это, но продолжаю получать ошибки, это то, что у меня есть сейчас, когда p0 близок к фактическому значению: return c1 * (((x1 - a1)/(d1 - x1))**(1/b1)) x_data = np.array([169,85, 172,15, 42,68, 65,13]) y_data = np.array([0,03, 0., 2.5, 0.9]) p0 = np.array([34., -1.5, 0.5, 170.]) popt, pcov = opt.curve_fit(f, x_data, y_data, p0, None, False, True, ftol = 0,00001) - person A.LeBrun; 18.02.2016
comment
Даже знаю, что я знаю, что эти значения производят значение r-квадрата 1, как говорится... Оптимальные параметры не найдены: количество вызовов функции достигло maxfev = 1000. Нужно ли мне менять ftol? - person A.LeBrun; 18.02.2016
comment
Я думаю, что оптимизация этой функции проблематична, потому что основание возведения в степень может иметь отрицательные значения. Если вы используете p0=np.array([15., -1./1.5, 0.5, 180.]), база положительна для всех наблюдаемых значений x, и curve_fit сходится (с предупреждениями) к довольно хорошему решению: [ 21.63157788 - 0,77630377 0,53499971 176,50006451] - person Francesc Torradeflot; 19.02.2016
comment
Большое спасибо! Это на самом деле производит то, что я ожидал! - person A.LeBrun; 19.02.2016

Это похоже на этот вопрос: Как решить пару нелинейных уравнений с помощью Python?

Но это, кажется, не дает мне правильного ответа. Попробуйте скопировать это и добавить больше очков?

from scipy.optimize import fsolve

def equations(p):
    a,b,c,d = p
    return ((c * ((176.5-a)/(d-176.5))**(1/b)),
            (c * ((161.333-a)/(d-161.333))**(1/b))-0.1,
            (c * ((66.1667-a)/(d-66.1667))**(1/b))-0.9,
            (c * ((40.833-a)/(d-40.833))**(1/b))-2.5)

a,b,c,d = fsolve(equations, (1,1,1,1))
print a,b,c,d

вывод:

1.0 1.0 1.0 1.0
person steven    schedule 17.02.2016
comment
Я добавил больше, но смог заставить работать нижнюю часть: (c * ((131.5-a)/(d-131.5))**(1/b))-0.3, (c * ((115.1667-a)/ (д-115,1667))**(1/б))-0,4, (с*((91,333-а)/(д-91,333))**(1/б))-0,6, (с*((175,333 -а)/(д-175,333))**(1/б))-0,03, (в * ((174,333-а)/(д-174,333))**(1/б))-0,05, (в * ((154,333-а)/(д-154,333))**(1/б))-0,15, - person A.LeBrun; 18.02.2016
comment
Было бы проще использовать SymPy? - person A.LeBrun; 18.02.2016
comment
a,b,c,d = fsolve(уравнения, (1,1,1,1)) SyntaxError: недопустимый синтаксис - person A.LeBrun; 18.02.2016
comment
Я вижу, да, похоже, что вы не можете переопределить эту функцию fsolve. Тогда попробуйте этот? (из того же потока, что и раньше) stackoverflow.com/a/8954984/5889975 - person steven; 18.02.2016
comment
Я пробовал это: from sympy import nsolve nsolve([c * ((176.5-a)/(d-176.5))**(1/b), c * ((66.1667-a)/(d-66.1667))* *(1/b)-0,9, c * ((40,833-a)/(d-40,833))**(1/b)-2,5, c * ((115,1667-a)/(d-115,1667))* *(1/b)-0,4], [a, b, c, d], [22, -1,5, 0,5, 175]), теперь я получаю сообщение об ошибке, говорящее, что a, b, c, d не определены, хотя Я знаю, что значения [22, -1,5, 0,5, 175] близки к тому, что должно быть. - person A.LeBrun; 18.02.2016