#13273: extension of sage.numerical.optimize.find_fit
-----------------------------------------------------------+----------------
Reporter: hackstein | Owner:
hackstein
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:
-----------------------------------------------------------+----------------
Comment (by hackstein):
I fixed some errors (see the updated code below), but it still don't fit
together with the leastsq routine. Thus any comments and improvement would
be highly appreciated as it is my first ticket.
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:
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)
for row in xrange(len(x1_data)):
fparams = [x1_data[row], x2_data[row]]+ params.tolist()
result[row]= (func[0](*fparams), func[1](*fparams))
return numpy.array([result[:,-2] - y1_data, result[:,-1] -
y2_data])
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:6>
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.
