I have been trying to use python's scipy.optimize
library to estimate the parameters of a model but without success so far. I tried to use scipy.optimize.leastsq
which uses the levenberg-marquardt algorithm. Unfortunately, it always fails to find the minimum of my model function even if I set the initial parameter guess very close to the best fit. Actually, it always returns the initial guessed parameters. So, I assume I am doing something wrong. My model is a simple circle but in order to make things even simpler only the radius is the actual parameter, the center of the circle in the data is known and hard-coded. The data is a 10x10 pixel float image with a circle with center 5,5 and radius 4. Actually, the data have been generated using the model I am trying to fit. So, a perfect fit exists. Here is my code:
import math
import numpy
import scipy.optimize
# ========================================================================
g_data_width = 10
g_data_height = 10
g_xo = 5.0
g_yo = 5.0
def evaluate_model01(x,y,r):
x2 = x*x
y2 = y*y
r2 = r*r
v = 0.0
if(x2 + y2 <= r2):
v = 20.0
return v
def model01(params,data_o):
data_r = numpy.zeros(g_data_height*g_data_width)
r = params[0]
for y in range(g_data_height):
for x in range(g_data_width):
xnew = x - g_xo
ynew = y - g_yo
data_r[y*g_data_width+x] = math.fabs(data_o[y,x]-evaluate_model01(xnew,ynew,r))
return data_r
# ========================================================================
g_data_o = numpy.array([[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, 0, 0, 0, 20, 0, 0, 0, 0],
[ 0, 0, 0, 20, 20, 20, 20, 20, 0, 0],
[ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0],
[ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0],
[ 0, 20, 20, 20, 20, 20, 20, 20, 20, 20],
[ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0],
[ 0, 0, 20, 20, 20, 20, 20, 20, 20, 0],
[ 0, 0, 0, 20, 20, 20, 20, 20, 0, 0],
[ 0, 0, 0, 0, 0, 20, 0, 0, 0, 0]],dtype=numpy.float32)
g_params = numpy.array([8.0])
print(scipy.optimize.leastsq(model01,g_params,args=(g_data_o),full_output=1))
# ========================================================================
I hard-coded the data in order to remove any data dependencies and allow the code to run out of the box on any machine with scipy
installed. What I don't quite understand is what the model01
function is supposed to return. According to the documentation it is supposed to return an array. An array of what? In my code I assumed that I have to return an array of the residual for each data point. Is that correct? My data is a 2D array because it is an image but my residual is a flattened 2D array of residual. Is that OK? Can someone tell me what exactly I am doing wrong? Can someone modify and fix the code? As I mentioned above, the code should run out of the box on any machine with scipy
and numpy
installed. In case what I want to achieve is not possible with scipy.optimize.leastsq
, can you suggest some other library which fits models using the levenberg-marquardt algorithm?