2

I am new to python and what I am trying to do is write an algorithm to solve for the 4 unknown parameters in the Rodbard Equation where we are relating a grayscale value measured using ImageJ to optical density calibration discs. This equation is nonlinear and is written as y = c*((x-a)/(d-x))^(1/b) where a, b, c, and d are unknown. I have the values of x and y for four point (176.5, 0), (161.333, 0.1), (66.1667, 0.9), and (40.833, 2.5). Below, I have posted my attempt to solve for these 4 unknowns. Any help to point me in the right direction would be greatly appreciated!

    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
A.LeBrun
  • 161
  • 1
  • 1
  • 11
  • I tried editing it by following the example but it is still not working for me. Sorry I couldn't get it to come up like code but here it is: import numpy as np from scipy.optimize import curve_fit def func(x, a, b, c, d): return [c * ((x - a)/(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) – A.LeBrun Feb 17 '16 at 21:52

2 Answers2

2

If you want to use curve_fit, you should do the following:

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 is the initial guess and if you do not inform it, an array of all ones will be assumed.

With the provided data I have tried with different params but I've not been able to find a solution.

I hope this helps. Good luck!

Francesc Torradeflot
  • 126
  • 5
  • I tried this but I keep get errors, this is what I have now where p0 is close to what the actual value is: import numpy as np import scipy.optimize as opt def f(x1, a1, b1, c1, d1): 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) – A.LeBrun Feb 17 '16 at 23:00
  • Even know I know these values produce an r-squared value of 1 it says... Optimal parameters not found: Number of calls to function has reached maxfev = 1000. Would I need to change ftol? – A.LeBrun Feb 18 '16 at 13:10
  • I think the optimization of this function is problematic because the base of the exponentiation can have negative values. If you use p0=np.array([15., -1./1.5, 0.5, 180.]) the base is positive for all observed x values and curve_fit converges (with warnings) to a quite good solution: [ 21.63157788 -0.77630377 0.53499971 176.50006451] – Francesc Torradeflot Feb 18 '16 at 22:35
  • Thank you so much! This actually produces something I was expecting! – A.LeBrun Feb 19 '16 at 02:25
0

This looks similar to this question: How to solve a pair of nonlinear equations using Python?

But this does not seem to give me the correct answer. Try copying this and adding more points?

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

output:

1.0 1.0 1.0 1.0
Community
  • 1
  • 1
steven
  • 226
  • 1
  • 13
  • I added more but I could get the bottom portion to work: (c * ((131.5-a)/(d-131.5))**(1/b))-0.3, (c * ((115.1667-a)/(d-115.1667))**(1/b))-0.4, (c * ((91.333-a)/(d-91.333))**(1/b))-0.6, (c * ((175.333-a)/(d-175.333))**(1/b))-0.03, (c * ((174.333-a)/(d-174.333))**(1/b))-0.05, (c * ((154.333-a)/(d-154.333))**(1/b))-0.15, – A.LeBrun Feb 17 '16 at 22:18
  • Would it be easier to use SymPy? – A.LeBrun Feb 17 '16 at 22:31
  • a,b,c,d = fsolve(equations, (1,1,1,1)) SyntaxError: invalid syntax – A.LeBrun Feb 17 '16 at 22:32
  • I see yeah it looks like you can't do overdetermined with that fsolve function. Try this one then? (from the same thread as before) http://stackoverflow.com/a/8954984/5889975 – steven Feb 17 '16 at 22:33
  • I tried this: 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]), now I get an error saying a, b, c, d are not defined even though I know the values [22, -1.5, 0.5, 175] are close to what the values should be. – A.LeBrun Feb 17 '16 at 22:47