Using lmfit minimize for the first time to fit z=f(x,y), it rune but coefficient always end at 0

144 Views Asked by At

I would like to fit z = f(x,y) using an objective function. I plan to fit more parameters later on, and lmfit sounded a nice abstraction to try.

For the sake of testing I created a controlled data set. The data is an array of coordinate X, coordinate Y, Vector X, Vector Y

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import gridspec

from scipy.optimize import leastsq
from lmfit import Parameters, fit_report, minimize

#Creat sample
xs = 10
ys = 10
s = 11
coo_x = np.linspace(-xs, xs,s)
coo_y = np.linspace(-ys, ys,s)
#Get all permutations of X,Ycoordinates
mesh = np.array(np.meshgrid(coo_x, coo_y))
coo = mesh.T.reshape(-1, 2)

header =  ["W_DesignPosX","W_DesignPosY","W_Registration_X_A","W_Registration_Y_A"]
transX = 3
transY = 0
angle = 0
magX = 0
magY = 0
orthX = 0

trans = np.linspace((transX,transY),(transX,transY),s*s)
rot = np.flip(coo, axis=1)*np.array ([-angle,angle])
mag = np.array([magX,magY])
orth = np.flip(coo, axis=1)*orthX/2

np.random.seed(seed=30)
random = np.random.normal(0,0.1, (s*s,2))
#random =  np.zeros((s*s,2))

#Compute data
test= np.concatenate((coo, trans+coo*mag+rot+orth+random), axis=1)
test_df = pd.DataFrame(data=test, columns=header)

In the test case above TransX = 3, all the other input are = 0 Running the minimize it should fit to the following A=3, B=0, C=0, but all end at 0 :(

def residual_x(param, x, y, data):
       
    A=params['A']
    B=params['B']
    C=params['C']
    
    model = A + B*x +C*y
    return (model-data)

params = Parameters()
params.add('A', value=0.0)
params.add('B', value=0.0)
params.add('C', value=0.0)


x,y =test[:,:2].T
reg_x = test[:,2]
out = minimize(residual_x,params, args = (x,y,reg_x))
print(fit_report(out))
print()
print(out.params.pretty_print())

I did eyeball the array and the quiver chart. The data has a horizontal vector.

def vector_summary(df,Design_x,Design_y,Reg_x,Reg_y,s=1):
    c = 'g'
    
    fig = plt.figure(figsize=(8, 4))
    grid = plt.GridSpec(2, 3,width_ratios=[1.5, 0.25,1])

    #Vector map
    ###########
    ax_q = fig.add_subplot(grid[:,0])
    X = list(df[Design_x])
    Y = list(df[Design_y])
    U = list(df[Reg_x])
    V = list(df[Reg_y])
    
    ax_q.quiver(X,Y,U,V,scale=0.04/s,color=c)
    ax_q.set_title("Vector map",fontsize=20)
    ax_q.set_xlabel('W_DesignPosX')
    ax_q.set_ylabel('W_DesignPosY')

    #ax_q.set_ylim([-20000,20000])

  
    #X_registration
    ###############
    ax_x= fig.add_subplot(grid[0,2])
    sns.histplot(df, x=Reg_x,ax=ax_x,color=c)
    ax_x.set_title("Reg_X",fontsize=20)
    
    #Y_registration
    ###############
    ax_y= fig.add_subplot(grid[1,2])
    sns.histplot(df, x=Reg_y,ax=ax_y,color=c)
    ax_y.set_title("Reg_Y",fontsize=20)
    plt.tight_layout()
    plt.show()

vector_summary(test_df,'W_DesignPosX','W_DesignPosY','W_Registration_X_A','W_Registration_Y_A',0.0005)

I am not a computer scientist and only have some instinct that my issue lies in the objective function. but I cannot point my finger on the issue.

Any advises would be appreciated! I am eager to learn. It is about the journey right ;-)

1

There are 1 best solutions below

0
On

You have a simple typo in your residuals function

def residual_x(param, x, y, data):

needs to be params and not param

def residual_x(params, x, y, data):
     

Hence instead of accessing the updated params from residuals (that exists in local namespace), it was just looking at your original params (that existed in global name space). There was no error raised because the minimizer doesn't check if the special keyword 'params' is passed or not, instead Python goes from local namespace to global namespace outside residuals and minimizer, and of course, that variable doesn't change.

Later when you were trying to access params you would get the original one you have created.