Next: Interpolation, Previous: Standard Nonlinear Models, Up: Curve Fitting

Calc's internal least-squares fitter can only handle multilinear
models. More precisely, it can handle any model of the form
‘`a f(x,y,z) + b g(x,y,z) + c h(x,y,z)`’, where ‘`a,b,c`’
are the parameters and ‘`x,y,z`’ are the independent variables
(of course there can be any number of each, not just three).

In a simple multilinear or polynomial fit, it is easy to see how
to convert the model into this form. For example, if the model
is ‘`a + b x + c x^2`’, then ‘`f(x) = 1`’, ‘`g(x) = x`’,
and ‘`h(x) = x^2`’ are suitable functions.

For most other models, Calc uses a variety of algebraic manipulations to try to put the problem into the form

Y(x,y,z) = A(a,b,c) F(x,y,z) + B(a,b,c) G(x,y,z) + C(a,b,c) H(x,y,z)

where ‘`Y,A,B,C,F,G,H`’ are arbitrary functions. It computes
‘`Y`’, ‘`F`’, ‘`G`’, and ‘`H`’ for all the data points,
does a standard linear fit to find the values of ‘`A`’, ‘`B`’,
and ‘`C`’, then uses the equation solver to solve for ‘`a,b,c`’
in terms of ‘`A,B,C`’.

A remarkable number of models can be cast into this general form.
We'll look at two examples here to see how it works. The power-law
model ‘`y = a x^b`’ with two independent variables and two parameters
can be rewritten as follows:

y = a x^b y = a exp(b ln(x)) y = exp(ln(a) + b ln(x)) ln(y) = ln(a) + b ln(x)

which matches the desired form with
‘`Y = ln(y)`’,
‘`A = ln(a)`’,
‘`F = 1`’, ‘`B = b`’, and
‘`G = ln(x)`’.
Calc thus computes the logarithms of your ‘`y`’ and ‘`x`’ values,
does a linear fit for ‘`A`’ and ‘`B`’, then solves to get
‘`a = exp(A)`’
and ‘`b = B`’.

Another interesting example is the “quadratic” model, which can be handled by expanding according to the distributive law.

y = a + b*(x - c)^2 y = a + b c^2 - 2 b c x + b x^2

which matches with ‘`Y = y`’, ‘`A = a + b c^2`’, ‘`F = 1`’,
‘`B = -2 b c`’, ‘`G = x`’ (the *-2* factor could just as easily
have been put into ‘`G`’ instead of ‘`B`’), ‘`C = b`’, and
‘`H = x^2`’.

The Gaussian model looks quite complicated, but a closer examination
shows that it's actually similar to the quadratic model but with an
exponential that can be brought to the top and moved into ‘`Y`’.

The logistic models cannot be put into general linear form. For these models, and the Hubbert linearization, Calc computes a rough approximation for the parameters, then uses the Levenberg-Marquardt iterative method to refine the approximations.

Another model that cannot be put into general linear
form is a Gaussian with a constant background added on, i.e.,
‘`d`’ + the regular Gaussian formula. If you have a model like
this, your best bet is to replace enough of your parameters with
constants to make the model linearizable, then adjust the constants
manually by doing a series of fits. You can compare the fits by
graphing them, by examining the goodness-of-fit measures returned by
`I a F`, or by some other method suitable to your application.
Note that some models can be linearized in several ways. The
Gaussian-plus-`d` model can be linearized by setting ‘`d`’
(the background) to a constant, or by setting ‘`b`’ (the standard
deviation) and ‘`c`’ (the mean) to constants.

To fit a model with constants substituted for some parameters, just store suitable values in those parameter variables, then omit them from the list of parameters when you answer the variables prompt.

A last desperate step would be to use the general-purpose
`minimize`

function rather than `fit`

. After all, both
functions solve the problem of minimizing an expression (the
‘`chi^2`’
sum) by adjusting certain parameters in the expression. The `a F`
command is able to use a vastly more efficient algorithm due to its
special knowledge about linear chi-square sums, but the `a N`
command can do the same thing by brute force.

A compromise would be to pick out a few parameters without which the
fit is linearizable, and use `minimize`

on a call to `fit`

which efficiently takes care of the rest of the parameters. The thing
to be minimized would be the value of
‘`chi^2`’
returned as the fifth result of the `xfit`

function:

minimize(xfit(gaus(a,b,c,d,x), x, [a,b,c], data)_5, d, guess)

where `gaus`

represents the Gaussian model with background,
`data`

represents the data matrix, and `guess`

represents
the initial guess for ‘`d`’ that `minimize`

requires.
This operation will only be, shall we say, extraordinarily slow
rather than astronomically slow (as would be the case if `minimize`

were used by itself to solve the problem).

The `I a F` [`xfit`

] command is somewhat trickier when
nonlinear models are used. The second item in the result is the
vector of “raw” parameters ‘`A`’, ‘`B`’, ‘`C`’. The
covariance matrix is written in terms of those raw parameters.
The fifth item is a vector of filter expressions. This
is the empty vector ‘`[]`’ if the raw parameters were the same
as the requested parameters, i.e., if ‘`A = a`’, ‘`B = b`’,
and so on (which is always true if the model is already linear
in the parameters as written, e.g., for polynomial fits). If the
parameters had to be rearranged, the fifth item is instead a vector
of one formula per parameter in the original model. The raw
parameters are expressed in these “filter” formulas as
‘`fitdummy(1)`’ for ‘`A`’, ‘`fitdummy(2)`’ for ‘`B`’,
and so on.

When Calc needs to modify the model to return the result, it replaces
‘`fitdummy(1)`’ in all the filters with the first item in the raw
parameters list, and so on for the other raw parameters, then
evaluates the resulting filter formulas to get the actual parameter
values to be substituted into the original model. In the case of
`H a F` and `I a F` where the parameters must be error forms,
Calc uses the square roots of the diagonal entries of the covariance
matrix as error values for the raw parameters, then lets Calc's
standard error-form arithmetic take it from there.

If you use `I a F` with a nonlinear model, be sure to remember
that the covariance matrix is in terms of the raw parameters,
*not* the actual requested parameters. It's up to you to
figure out how to interpret the covariances in the presence of
nontrivial filter functions.

Things are also complicated when the input contains error forms.
Suppose there are three independent and dependent variables, ‘`x`’,
‘`y`’, and ‘`z`’, one or more of which are error forms in the
data. Calc combines all the error values by taking the square root
of the sum of the squares of the errors. It then changes ‘`x`’
and ‘`y`’ to be plain numbers, and makes ‘`z`’ into an error
form with this combined error. The ‘`Y(x,y,z)`’ part of the
linearized model is evaluated, and the result should be an error
form. The error part of that result is used for
‘`sigma_i`’
for the data point. If for some reason ‘`Y(x,y,z)`’ does not return
an error form, the combined error from ‘`z`’ is used directly for
‘`sigma_i`’.
Finally, ‘`z`’ is also stripped of its error
for use in computing ‘`F(x,y,z)`’, ‘`G(x,y,z)`’ and so on;
the righthand side of the linearized model is computed in regular
arithmetic with no error forms.

(While these rules may seem complicated, they are designed to do
the most reasonable thing in the typical case that ‘`Y(x,y,z)`’
depends only on the dependent variable ‘`z`’, and in fact is
often simply equal to ‘`z`’. For common cases like polynomials
and multilinear models, the combined error is simply used as the
‘`sigma`’
for the data point with no further ado.)

It may be the case that the model you wish to use is linearizable,
but Calc's built-in rules are unable to figure it out. Calc uses
its algebraic rewrite mechanism to linearize a model. The rewrite
rules are kept in the variable `FitRules`

. You can edit this
variable using the `s e FitRules` command; in fact, there is
a special `s F` command just for editing `FitRules`

.
See Operations on Variables.

See Rewrite Rules, for a discussion of rewrite rules.

Calc uses `FitRules`

as follows. First, it converts the model
to an equation if necessary and encloses the model equation in a
call to the function `fitmodel`

(which is not actually a defined
function in Calc; it is only used as a placeholder by the rewrite rules).
Parameter variables are renamed to function calls ‘`fitparam(1)`’,
‘`fitparam(2)`’, and so on, and independent variables are renamed
to ‘`fitvar(1)`’, ‘`fitvar(2)`’, etc. The dependent variable
is the highest-numbered `fitvar`

. For example, the power law
model ‘`a x^b`’ is converted to ‘`y = a x^b`’, then to

fitmodel(fitvar(2) = fitparam(1) fitvar(1)^fitparam(2))

Calc then applies the rewrites as if by ‘`C-u 0 a r FitRules`’.
(The zero prefix means that rewriting should continue until no further
changes are possible.)

When rewriting is complete, the `fitmodel`

call should have
been replaced by a `fitsystem`

call that looks like this:

fitsystem(Y,FGH,abc)

where `Y` is a formula that describes the function ‘`Y(x,y,z)`’,
`FGH` is the vector of formulas ‘`[F(x,y,z), G(x,y,z), H(x,y,z)]`’,
and `abc` is the vector of parameter filters which refer to the
raw parameters as ‘`fitdummy(1)`’ for ‘`A`’, ‘`fitdummy(2)`’
for ‘`B`’, etc. While the number of raw parameters (the length of
the `FGH` vector) is usually the same as the number of original
parameters (the length of the `abc` vector), this is not required.

The power law model eventually boils down to

fitsystem(ln(fitvar(2)), [1, ln(fitvar(1))], [exp(fitdummy(1)), fitdummy(2)])

The actual implementation of `FitRules`

is complicated; it
proceeds in four phases. First, common rearrangements are done
to try to bring linear terms together and to isolate functions like
`exp`

and `ln`

either all the way “out” (so that they
can be put into `Y`) or all the way “in” (so that they can
be put into `abc` or `FGH`). In particular, all
non-constant powers are converted to logs-and-exponentials form,
and the distributive law is used to expand products of sums.
Quotients are rewritten to use the ‘`fitinv`’ function, where
‘`fitinv(x)`’ represents ‘`1/x`’ while the `FitRules`

are operating. (The use of `fitinv`

makes recognition of
linear-looking forms easier.) If you modify `FitRules`

, you
will probably only need to modify the rules for this phase.

Phase two, whose rules can actually also apply during phases one
and three, first rewrites `fitmodel`

to a two-argument
form ‘`fitmodel( Y, model)`’, where

`fitmodel`

is changed to a
four-argument `fitsystem`

, where the fourth argument is
Phase three converts a sum of items in the `model` to a sum
of ‘`fitpart( a, b, c)`’ terms which represent
terms ‘

`fitpart`

s
with equal Phase four moves the `fitpart`

terms into the `FGH` and
`ABC` vectors. Also, some of the algebraic expansions that
were done in phase 1 are undone now to make the formulas more
computationally efficient. Finally, it calls the solver one more
time to convert the `ABC` vector to an `abc` vector, and
removes the fourth `model` argument (which by now will be zero)
to obtain the three-argument `fitsystem`

that the linear
least-squares solver wants to see.

Two functions which are useful in connection with `FitRules`

are ‘`hasfitparams(x)`’ and ‘`hasfitvars(x)`’, which check
whether ‘`x`’ refers to any parameters or independent variables,
respectively. Specifically, these functions return “true” if the
argument contains any `fitparam`

(or `fitvar`

) function
calls, and “false” otherwise. (Recall that “true” means a
nonzero number, and “false” means zero. The actual nonzero number
returned is the largest `n` from all the ‘`fitparam( n)`’s
or ‘

The `fit`

function in algebraic notation normally takes four
arguments, ‘`fit( model, vars, params, data)`’,
where

If `params` is omitted, the parameters are all variables in
`model` except those that appear in `vars`. If `vars`
is also omitted, Calc sorts all the variables that appear in
`model` alphabetically and uses the higher ones for `vars`
and the lower ones for `params`.

Alternatively, ‘`fit( modelvec, data)`’ is allowed
where

If Calc is unable to do the fit, the `fit`

function is left
in symbolic form, ordinarily with an explanatory message. The
message will be “Model expression is too complex” if the
linearizer was unable to put the model into the required form.

The `efit`

(corresponding to `H a F`) and `xfit`

(for `I a F`) functions are completely analogous.