# Non-linear regression with arbitrary precision arithmetic

I'm looking to do something like this:

R = RealField(1000)
data = [[R(1),R(5)],[R(5),sqrt(R(6))],[R(8),R(9)]]
var('a, b, c, x')
model(x) = a*x*x + b*x + c
find_fit(data, model)


However when I run this a,b and c seem to have been truncated to a double. How would I get arbitrary precision out of find_fit?

edit retag close merge delete

Sort by » oldest newest most voted

By looking at the source code directly in Sage:

sage: find_fit??


or online here one sees that the first thing that happens to your data

sage: R = RealField(1000)
sage: data = [[R(1),R(5)],[R(5),sqrt(R(6))],[R(8),R(9)]]
[[1.000000000000...00, 5.000000000000...00],
[5.000000000000...00, 2.449489742783...84],
[8.000000000000...00, 9.000000000000...00]]


is

sage: import numpy
sage: data = numpy.array(data, dtype = float)
sage: data
array([[ 1.        ,  5.        ],
[ 5.        ,  2.44948974],
[ 8.        ,  9.        ]])


Therefore, all the precision you gave is lost from the start :(

Then, few lines later, your problem is reduced to a least square problem and the function scipy.optimize.leastsq from scipy is used to solve your problem (which in turn uses MINPACK's lmdif and lmder algorithms according to the doc):

from scipy.optimize import leastsq
estimated_params, d = leastsq(error_function, initial_guess, args = (x_data, y_data))


Then, by looking at the documentation of leastsq:

sage: from scipy.optimize import leastsq
sage: leastsq?


You can see that some options that you may be searching for are available:

   ftol : float, optional
Relative error desired in the sum of squares.
xtol : float, optional
Relative error desired in the approximate solution.
gtol : float, optional
Orthogonality desired between the function vector and the
columns of the Jacobian.


The documentation leastsq? also gives you the default values that you may want to change: ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0.

So now the next steps are that you redo on your own the reduction to the least square problem done in find_fit without losing the precision by not casting entries to float. Then, you change the default values for ftol, xtol, gtol in the call to leastsq to see if you are able achieve your goal in increasing the precision of the result.

Then, if you succeed, you may consider proposing an improvement to the Sage function find_fit.

more