Progress Report

This blog post is concerned with the milestones achieved and some upcoming features in statsmodels

NONLINEAR MODELS

The model implemented is

y = f(x,θ) + e

where the y is one dimensional endogenous data matrix, f is the nonlinear function, x is exogenous data  matrix, θ is parameter matrix and 'e' is the noise in data. 

Estimation

The estimation of parameters is done using the 'leastq' method from scipy.optimize which minimizes the sum of squares of residuals. We subclass the model class 'NonlinearLS' and provide the 'expr' function which calculates 'f' in the above expression using the parameter values and exogenous data provided to it. It is encouraged that the user provides the analytical derivative of the given function by defining 'jacobian' function in the similar way as 'expr'.

Testing
For testing purposes we used the 'Misra1a' model from NIST data. Details regarding this given in previous post. In summary, we obtained satisfactory results as compared to 'Gretl' which uses the same minpack module used by scipy.

Miscellaneous Features

• Parameters calculated at each iteration by the algorithm can be viewed using view_iter() method
• Prediction table with confidence intervals for each predicted value of endogenous data using prediction_table(alpha) method

Example
A complete example can be viewed here https://github.com/divyanshubandil/statsmodels/commit/db2e388232303323cc9bb36e0fe9f682892993ba

ROBUST NONLINEAR MODELS
I have been working on the M-estimation of nonlinear models for some time now. The best research paper I found having all the tests, computational algorithm and simulation data is here. http://www.tandfonline.com/doi/abs/10.1080/03610920802074836
Recently, I have been able to implement the algorithm in my first commit regarding this topic here https://github.com/divyanshubandil/statsmodels/commit/1745e02b45ebe3f83a8e0d55f477fcef33621d6f
Now I am working with testing this model for the 'Numerical example' given in the paper.

Nonlinear models - Results and Tests

After computing parameters for nonlinear models using least squares (LQ) method of estimation, next step was to add result statistics and writing tests comparing results from other statistical package to those from statsmodels.

Model and Results
The model used for analysis was the 'Misra1a' model selected from NIST nonlinear models. Dataset misra was addded to datasets directory. Parameter estimates and standard error of residuals were provided in the NIST data file. Other results of nonlinear regression analysis of the given model were obtained from Gretl. Reason for choosing gretl was that it uses the same Minpack module for LQ estimation that is used in scipy.optimize. All the results were stored in python class Misra1a.

Testing
Four sets of 18 tests ( 2 start values with and without the jacobian ) provided in test_Misra1a.py were run and we obtained exact matches upto 4 decimal places.

Weighted NLS
Weights (1 for the first seven observations and 0.49 for the next seven observations) were provided while fitting the data.
Gretl does not provide the feature of weighted nls and hence we had to use results from R which has the provision to do so. Test class for weighted NLS was added to test_Misra1a.py and 14 tests were added. Only one test regarding confidence intervals of parameters failed which is only a precision issue (No digits after the decimal match).

NIST Tests - Part 1

As this is my first commit after getting selected for GSoC 2012, I would summarise the status of the project so far followed by progress made this week.

Project Status

• Module declaring Nonlinear Model class using scipy.leastsq() was provided by my mentor, which implements nonlinear least squares method for data fitting. It also included tests against linear regression results.
• During GSoC prelims, I added a mechanism to provide an external jacobian calculating function, both analytically and numerically, which helped to save the parameter values at each iteration of minimizaion algorithm.
• The module test_nonlinearls.py was also updated, including a test for the above feature.
• Now, we are testing the NonlinearLS class derived models against NIST models.

Progress This Week
At present, I am writing tests using NIST data for a corresponding model and comparing the results obtained against the certified values.
A module test_generator.py is used for reading NIST data files and create a corresponding test file. The test_generator module reads the model name, its function expression and analytical jacobian expression from a configuration file model_list.cfg and then creates a test class for the corresponding model.
The objective is to first check the resulting parameter values, which are calculated using the LM algorithm, followed by other result statistics.

Results
Models labeled as of lower difficulty are implemented at present. Most of them them give satisfactory matches for parameter results upto 6-7 decimal places. The code can be viewed here. Also, the tests can be run using nosetests. As I have not uploaded the NIST data files, the test_generator.py module will not be able to give the output test files.

I will continue to develop test code for all the models. NIST results give us only parameter values and standard deviations. After we compare both of them, I will add more result statistics for nonlinear models.

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.