Hi Per,

The fitting routines in Statistics all work in the floating point domain only. That means that the result you get will be accurate in a floating point sense - the inputs will be regarded as being specified up to a certain number of digits, and the result will be accurate up to a certain number of digits as well.

As with all other floating point computations, you can set the required accuracy by changing the *Digits* environment variable. If you want to have a reasonable sense of confidence that all *n* digits that you see in the answer are correct (that is, that this *n*-digit number provides a closer fit than any other *n*-digit number - not that it is an *n*-digit approximation to a perfect fit: such a perfect fit may not exist!), you should set *Digits* to a higher value initially (maybe *n* + 3 or so, potentially depending on the number of data points you have), then round the result to *n* digits. You could do that as follows:

oldDigits := Digits:

Digits := Digits + max(3, ilog10(numelems(X)));

answer := PowerFit(X, Y, x);

Digits := oldDigits:

answer := evalindets(answer, 'float', x -> evalf(SFloat(x)));

That last line is complicated by the split personality of floating point computations in Maple. Maple supports floating point numbers with an arbitrary number of digits, which are implemented in software, and hardware-precision floating point numbers which are more efficient but provide limited accuracy (15 decimal digits on all platforms, currently). It tries to use hardware floats in many situations (including in the Statistics package whenever the setting of *Digits* allows it), but to round a number to a given accuracy, we have to force it to be a software floating point number.

Another meaning of your question might be that you're looking for a symbolic answer, not a floating point one. In that case, first of all you would definitely not want to have any floating point numbers in your input. Furthermore, the Statistics package can't help you for this case. Thankfully, it's fairly straightforward to implement the required linear algebra yourself:

exactPowerFit := proc(X :: Vector[column],

Y :: Vector[column],

x :: name,

$)

local

A :: Vector[column],

M :: Matrix,

lnX :: Vector[column],

lnXsum :: algebraic,

lnY :: Vector[column],

n :: nonnegint,

v;

n := numelems(X);

if numelems(Y) <> n then

error "dependent and independent data should have equally many data points, "

"received: %1 and %2, respectively", numelems(Y), n;

elif n < 2 then

error "need at least two data points, received: %1", n;

elif hastype([X, Y], 'float') then

WARNING("this procedure is not meant to be used with floating point values");

end if;

# Precomputation.

lnX := ln~(X);

lnY := ln~(Y);

lnXsum := add(v, v in lnX);

M := <n, lnXsum | lnXsum, add(v^2, v in lnX)>;

A := <add(v, v in lnY), add(lnX[v] * lnY[v], v = 1 .. n)>;

A := LinearAlgebra:-LinearSolve(M, A);

return exp(A[1]) * x^A[2];

end proc;

Now you can do

X := <seq(i, i = 1 .. 6)>;

Y := <seq(i^2, i = 1 .. 6)>;

exactPowerFit(X, Y, x);

and get a symbolic result. Or more interestingly,

Y := <seq(i^2 + 1, i = 1 .. 6)>;

and find a large expression for the coefficient and exponent that give the closest fit in a least-squares sense (in the log domain). It's a little too big to print here comfortably.

Hope this helps,

Erik Postma

Maplesoft.