#13273: extension of sage.numerical.optimize.find_fit
-----------------------------------------------------------+----------------
Reporter: hackstein | Owner:
Type: enhancement | Status:
needs_review
Priority: major | Milestone:
sage-5.3
Component: numerical | Resolution:
Keywords: fitting, interpolation, two-dimensional | Work issues:
enhancement
Report Upstream: N/A | Reviewers:
evanandel, was, ncohen
Authors: Urs Hackstein | Merged in:
Dependencies: | Stopgaps:
-----------------------------------------------------------+----------------
Changes (by hackstein):
* owner: hackstein =>
Comment:
It runs through now, but it still needs careful review. Any comments and
improvements are very welcome.
def find_fit2dim(data, model, initial_guess=None, parameters=None,
variables=None, solution_dict=False):
r"""
Finds numerical estimates for the parameters of the two-dimensional
function model to
give a best fit to data.
INPUT:
- ``data`` -- A two dimensional table of floating point numbers of the
form `[[x_{1,1}, x_{1,2}, f_{1,1}, f_{1,2}],
[x_{2,1}, x_{2,2}, f_{2,1}, f_{2,2}],
\ldots,
[x_{n,1}, x_{n,2}, f_{n,1}, f_{n,2}]]` given as either a
numpy array.
- ``model`` -- Either a pair of symbolic expressions, of symbolic
functions, or
of Python function. ``model`` has to be a function of the variables
`(x_1, x_2, \ldots, x_k)` and free parameters
`(a_1, a_2, \ldots, a_l)`.
- ``initial_guess`` -- (default: ``None``) Initial estimate for the
parameters `(a_1, a_2, \ldots, a_l)`, given as either a list, tuple,
vector or numpy array. If ``None``, the default estimate for each
parameter is `1`.
- ``parameters`` -- (default: ``None``) A list of the parameters
`(a_1, a_2, \ldots, a_l)`. If model is a symbolic function it is
ignored, and the free parameters of the symbolic function are used.
- ``variables`` -- (default: ``None``) A list of the variables
`(x_1, x_2, \ldots, x_k)`. If model is a symbolic function it is
ignored, and the variables of the symbolic function are used.
- ``solution_dict`` -- (default: ``False``) if ``True``, return the
solution as a dictionary rather than an equation.
EXAMPLES:
First we create some data points of a function f
\colon\R^{2}\to\R^{2}. Let f(x,y)=(x-1,3*y+2)
sage: data=np.array([[1,0,0,2],[1,1,0,5],[2,1,1,5],[2,3,1,11]])
We define a function with free parameters "a" and "b":
sage: model(x,y)=(a+x,b+3*y)
We search for the parameters that give the best fit to the data:
find_fit2dim(data,model)
[a == -1.0000000064095422, b == 1.9999999999999998]
ALGORITHM:
Uses ``scipy.optimize.leastsq`` which in turn uses MINPACK's lmdif and
lmder algorithms.
"""
import numpy
if not isinstance(data, numpy.ndarray):
try:
data = numpy.array(data, dtype = float)
except (ValueError, TypeError):
raise TypeError, "data has to be a numpy array"
elif data.dtype == object:
raise ValueError, "the entries of data have to be of type float"
if data.ndim != 2:
raise ValueError, "data has to be a two dimensional table of
floating point numbers"
from sage.symbolic.expression import Expression
if isinstance(model[0], Expression) and
isinstance(model[1],Expression):
if variables is None:
variables = list(model[0].arguments())
if parameters is None:
L=list(model[0].variables())
M=list(model[1].variables())
for v in L:
if v in M:
M.remove(v)
parameters=L+M
for v in variables:
parameters.remove(v)
if data.shape[1] != len(variables) + 2:
raise ValueError, "each row of data needs %d entries, only %d
entries given" % (len(variables) + 2, data.shape[1])
if parameters is None or len(parameters) == 0 or \
variables is None or len(variables) == 0:
raise ValueError, "no variables given"
if initial_guess == None:
initial_guess = len(parameters) * [1]
if not isinstance(initial_guess, numpy.ndarray):
try:
initial_guess = numpy.array(initial_guess, dtype = float)
except (ValueError, TypeError):
raise TypeError, "initial_guess has to be a list, tuple, or
numpy array"
elif initial_guess.dtype == object:
raise ValueError, "the entries of initial_guess have to be of type
float"
if len(initial_guess) != len(parameters):
raise ValueError, "length of initial_guess does not coincide with
the number of parameters"
if isinstance(model[0], Expression) and
isinstance(model[1],Expression):
var_list = variables + parameters
var_names = map(str, var_list)
func = (model[0]._fast_float_(*var_names),
model[1]._fast_float_(*var_names))
else:
func = model
def function(x1_data, x2_data, params):
s=(len(x1_data), 2)
result = numpy.zeros(s)
for row in xrange(len(x1_data)):
fparams = numpy.hstack((x1_data[row],
x2_data[row],params)).tolist()
result[row] = (func[0](*fparams),func[1](*fparams))
return result
def error_function(params, x1_data, x2_data, y1_data, y2_data):
t=(len(x1_data),2)
result = numpy.zeros(t)
result2 = numpy.zeros(len(x1_data))
for row in xrange(len(x1_data)):
fparams = [x1_data[row], x2_data[row]]+ params.tolist()
result[row]= (func[0](*fparams), func[1](*fparams))
result2[row]=(result[row,-2] - y1_data[row])*(result[row,-2] -
y1_data[row])+ (result[row,-1] - y2_data[row])*(result[row,-1] -
y2_data[row])
return result2
x1_data = data[:, -4]
x2_data = data[:, -3]
y1_data = data[:, -2]
y2_data = data[:,-1]
from scipy.optimize import leastsq
estimated_params, d = leastsq(error_function, initial_guess, args =
(x1_data, x2_data, y1_data, y2_data))
if isinstance(estimated_params, float):
estimated_params = [estimated_params]
else:
estimated_params = estimated_params.tolist()
if solution_dict:
dict = {}
for item in zip(parameters, estimated_params):
dict[item[0]] = item[1]
return dict
return [item[0] == item[1] for item in zip(parameters,
estimated_params)]
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/13273#comment:7>
Sage <http://www.sagemath.org>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en.
