Analytical Jacobian Use

In the previous post I mentioned the issue of user provided analytical derivative not giving correct results. With the help of Josef, the issue is resolved now. It turned out that, before passing on to leastsq the jacobian values had to transformed into those of "whitened" exogenous data. Here is an example of a typical Nonlinear Model Class:

class Myfunc(NonlinearLS):

    def expr(self, params, exog=None):
        if exog is None:
            x = self.exog
        else:
            x = exog

        x0, x1 = x[:,0], x[:,1]
        a, b, c = params
        return a + b*x0 + c*x1

    def jacobian(self,params, exog=None):
        if exog is None:
            x = self.exog
        else:
            x = exog

        x0, x1 = x[:,0], x[:,1]
        a, b, c = params
        a1 = np.ones(len(x0),)
        b1 = x0
        c1 = x1
        jacob = np.column_stack((a1,b1,c1))
        return jacob
 

Some points to note:

• The above class is for linear model y = a + b*x0 + c*x1 + e, where a,b and c are the parameters to be estimated.
• The expr is the 'expression' of the function of exogenous variable. Simply, it is 'f' in y = f(x) + e. Here f is a + b*x0 + c*x1.
• The jacobian is the analytical derivative matrix of 'f' given above with respect to the parameters to estimated.
• The jacobian function is optional. If it is not provided, jacobian is calculated using numerical derivative.

A test has been included in test_nonlinearls.py for above class. The results obtained after nonlinear least squares fitting were:

Next we move on inlcude the result statistics for NonLinearLS class.

Jacobian Using Analytical Derivatives

Testing the numerical derivative jacobian

I updated the code and ran the test_nonlinearls.py. It tests the results obtained from NonlinearLS fitting with WLS fitting results. All the tests gave fine results for all test models using jacobian calculated explicitly using numerical derivative.

For the case of linear model, y = a + b*x0 + c*x1 + e ,where a,b,c are the parameters to be estimated and e is the noise. The output is as given below (.params and .bse used for the following outputs)

leastsq Parameters [ 0.85542169 -0.72289157  1.80722892]
leastsq Standard Errors [ 0.69147062  0.85276594  2.04464611]

WLS Parameters [ 0.85542169 -0.72289157  1.80722892]
WLS Standard Errors [ 0.69147062  0.85276594  2.04464611]
 

Providing the Exact Jacobian
As mentioned in the previous post, I introduced a function 'jacobian' to let the user provide the analytical formula for calculating the jacobian.

I tested it using the linear regression model: y = a + b*x0 + c*x1 + e. The output is as given below (.params and .bse used for the following outputs):

leastsq Parameters [ 0.72754286 -0.81228571  2.15571429]
leastsq Standard Errors [ 0.69737916  0.86005273  2.06211739]

WLS Parameters [ 0.85542169 -0.72289157  1.80722892]
WLS Standard Errors [ 0.69147062  0.85276594  2.04464611]

The output values clearly don't match.

However, when the WLS estimates are provided as starting values for the leastsq method, the output is as follows.

WLS Parameters [ 0.85542169 -0.72289157  1.80722892]
WLS Standard Errors [ 0.69147062  0.85276594  2.04464611]

leastsq Parameters [ 0.85542169 -0.72289157  1.80722892]
leastsq Standard Errors [ 0.69147062  0.85276594  2.04464611]

The output values do match now.

The following points can be inferred:

• leastsq using LM algorithm gives quite different values for parameters than WLS. However the standard errors are close for both methods.This is the case when the jacobian is calculated analytically for leastsq.  In case of numerical derivative jacobian, the results from WLS and leastsq match exactly. This point needs to be tested with some other results to understand this behaviour.
• The start values for fitting using nonlinearls can be provided with WLS estimates. This gives a good result atleast in the case of linear regression model.

Next, we will test the code written now with some 'real' nonlinear models and compare it with results with other statistical packages.

Parameter estimates

Jacobian issue resolved

The discrepancy mentioned in my previous blog post is resolved. The Jacobian calculation formula was working fine. It turned out some code debugging was required. Now the jacobian is calculated explicitly for each step using forward differences method. It is then passed on to leastsq function for the calculation of parameters.
Next I will build some mechanism to let the user pass his own jacobian calculating function.
The code edits can be viewed here.


Viewing parameter estimates

The explicit jacobian approximations helps to save the parameter estimates at each iteration of the LM algorithm. These parameter estimates are passed on to an instance of the results class.
Presently the user can view the estimates using view_iter() function of the result instance. A screenshot of the parameter estimates table is given below.

Starting off with NonLinearModel Class

Here are some of the points I have been working on-

Base Class
The nonlinear regression models can be expressed as

y = f(x,parameters) + e

Here f is the nonlinear function of explanatory variable x and different parameters.

NonLinearModel, the base class for all nonlinear model needs to have a function which contains the expression of 'f' and calculates the response for a set of given data points and parameters.
The expr function in NonLinearModel is introduced for the above reason.An example using NonLinearLS class is given here.
Slightly refactored 'nonlinls.py'  is given here.

Calculation of Jacobian
The leastsq() procedure is based on the minpack/lmder.f and minpack/lmdif.f. These modules use Levenberg-Marquardt algorithm for nonlinear least squares optimization. Both the modules use forward-differences method to calculate the Jacobian. It gives the following output for the above example-

 

However, calculating the jacobian explicitly using the numdiff/approx_fprime1 module in sandbox and passing on to 'leastsq' gave different results as following.

Currently working on it to find a solution.
The reason for explicitly calculating jacobian is to replace the forward differences method by better methods like automatic differentiation or n-point numerical derivative. Jacobian values are needed to keep track and prevent any deviation from the expected descent direction.

Nonlinear Regression Classes

Suggestions about class inheritance to be followed during the course of project -

NonLinearModel - Base Class for all non linear models. Derived from Class Model
Derived Classes-

• NonlinearLS - Nonlinear Least Squares class similar to OLS and WLS. Use Levenberg-Marquardt Algorithm for curve fitting. Currently using scipy.optimize.leastsq() for this purpose.
• NonlinearRobust - Class for calculating nonlinear robust estimators. IRLS algorithm or its variant used for calculation.

NonLinearRegressionResults -Base Class for holding the estimates and predictions. Derived from Regression Results
Derived Classes-

• NonlinearLSResults - For holding results from fitting NonlinearLS
• NonlinearRobustResults - For holding results from NonlinearRobust. Furthur classes may be derived for different estimators.

The user will have to provide an expression in the expr method of a class derived from the respective NonLinear classes for fitting the data.
Optional expressions will be provided.

GSOC Application

Key points of the application

• Implementation of NonLinearModel Class with development of nonlinear curve fitting machinery
• Implementation of NonLinearModel_Results Class for storing results
• Integrating OLS, WLS and RLM model classes with nonlinear framework
• Introducing MM-estimator and tau-estimator for robust linear models
• Writing test suite with suitable datasets for different estimators