# Using FindFit to fit $a\,b^t$: how to avoid introducing complex numbers?

I'd like to find a model $f(t)=a\,b^t$ which matches the following data

data = {{0.0, 100.0}, {0.02, 81.87}, {0.04, 67.03}, {0.06, 54.88},
{0.08, 44.93}, {0.10, 36.76}};


I tried this:

FindFit[data, a*b^t, {a, b}, t]


But that results in an error:

FindFit::nrlnum: The function value {-4.38571+0. I,<<4>>,136.386+<<1>>} is not
a list of real numbers with dimensions {6} at {a,b} = {95.6143,-626.386}. >>


The answer to the problem is $a \approx 100.01244$ and $b \approx 0.000045146$.

Is there a way to use FindFit[] to get the correct answer? Should I be using some other function?

This is from problem 3.4.85 in Calculus/Stewart 6E stated as follows:

The flash unit on a camera operates by storing charge on a capacitor and releasing it suddenly when the flash is set of f. The following data describe the charge remaining on the capacitor (measured in microcoulombs, µC) at time (measured in seconds):

...

Use a graphing calculator or computer to find an exponential model for the charge.

-
In newer versions of Mathematica you can also use NonlinearModelFit[]. If you're fitting an exponential, consider taking the logarithm of values. – Szabolcs Feb 17 '12 at 9:47
@Szabolcs: taking the logarithm is good for getting initial values for subsequent polishing with a nonlinear regression method, but is poor for actually determining regression parameters, since the transformation distorts the error in the data as well. – J. M. Feb 17 '12 at 10:23
@J.M. That is true. Which way is better depends on the application. For certain applications it may be preferable to have the error transformed as well, for others it may be harmful. – Szabolcs Feb 17 '12 at 10:26
can't you try to simplify your code to the minimum that gives the error, then post that code? as it is one would have to guess what the problem is – acl Aug 25 '12 at 15:21
Ignoring the \[VeryThinSpace] issue, I think the issue lies in the message itself: your data isn't real (as in it contains complex numbers). Try LinearModelFit, instead. – rcollyer Aug 25 '12 at 16:53

## The solution

It's usually a very very very very good idea to plot your data to see what shape it roughly has in the beginning. In your case, it looks like this:

... wait, that's not an exponentially rising function! Of course $a\,b^t$ won't fit that, but $a\,b^{-t}$ might work, let's try:

FindFit[data, a*b^(-t), {a, b}, t]

{a -> 100.004, b -> 22099.8}


## Some further remarks

I have of course used the fact that $b^t=\left(\frac1b\right)^{-t}$ here. You could have fitted with your function as well, only that Mathematica chose inappropriate starting parameters. I do not know how Mathematica determines these, but I think assuming it uses something around 1 is reasonable. My impression is always that it favors large to small numbers, but I wouldn't count on that.

You can specify starting parameters like this:

FindFit[data, a*b^t, {{a, aStart}, {b, bStart}}, t]


This will tell Mathematica to use aStart and bStart initially; it will then continue to slighly alter these values to see whether the fit gets better.

"Wait, but if I don't know these parameters at all, how can I find good starting values out" - here's one of the great powers of Mathematica: easy interactivity. Simply code a Manipulate environment where you can set the parameters yourself, like this:

Manipulate[
Show[
ListLinePlot[data],
Plot[a b^(-t), {t, 0, .10}, PlotStyle -> Red]
],
{a, 0, 200},
{b, 0, 50000}
]


Here you can play around with the parameters, until your function (in red) roughly matches the data. Read the values of the parameters, use them as starting values, get delicious fit.

Let's say you found out your starting parameters using Manipulate or from some solution manual (cough). Say your approximate values are $a\approx 100$, $b\approx10^{-5}$:

FindFit[data, a*b^t, {{a, 100}, {b, 10^(-5)}}, t]

{a -> 100.004, b -> 0.0000452493}


That's your textbook solution right there.

Another thing worth mentioning is NonlinearModelFit, the big brother of FindFit. It features not only fitting, it also generates a huge amount of extra data for statistical analysis, for example error estimates of the parameters and plotting confidence bands. Here's an example where I fitted a Gaussian point set with the theoretical curve:

-
Hi David, thanks for the help! When I plot the curve: Plot[a*b^t /. {a -> 100.01244, b -> 0.000045146}, {t, 0, 0.1}] using the non-negative exponent with the values of a and b given, it appears to match the data pretty well. :-) – dharmatech Feb 17 '12 at 5:49
I added some paragraphs on this. Basically you've had the problem that Mathematica chose inappropriate starting parameters for the fit. It's usually a good idea to always give this information to the program. – David Feb 17 '12 at 6:19
Excellent info and answer @David! Accepted! – dharmatech Feb 17 '12 at 7:07
Just to add that NonlinearModelFit's parameter determination can be made substantially more robust in version 8 using the NMinimize method (this should have been in version 7, but did not work due to a bug). For example, I have quite successfully used Method -> {NMinimize, Method -> {"DifferentialEvolution", "PostProcess" -> {FindMinimum, Method -> "QuasiNewton"}} to fit a (more complicated than shown here) exponential model to experimental data. Differential evolution is a particularly robust minimization method. – Oleksandr R. Feb 17 '12 at 10:20
I agree with Oleksandr; for particularly complicated models, it's a good idea to use a global optimization strategy like differential evolution and then a subsequent polishing with a more traditional unconstrained optimization method like quasi-Newton. This is particularly beneficial when the model is complicated enough that one suspects that there are a lot of local optima (e.g. fitting to a combination of trigonometric functions with different frequencies)... – J. M. Feb 17 '12 at 10:49

One way to address this issue, and one that I think works well, is to fit to an explicitly complex model. While such an approach initially seems unnecessary for this real-valued problem, the crucial point is that, when complex values are encountered, you need a mechanism (that FindFit can handle without any difficulties) for penalizing the unwanted imaginary part that doesn't involve changing the essential meaning of the model. Just adding Abs or Re--a common approach that usually works reasonably well, I admit--will not suffice in more difficult cases. If you're unable to give reasonable initial guesses for the parameters, the fit is then very likely to stray off the real line and into some false minimum somewhere in the complex plane.

It's quite easy to convert a real model to a complex one automatically--and, to a large extent, you don't even have to think about it. Here's a way to do that.

data = {
{0.0, 100.0}, {0.02, 81.87}, {0.04, 67.03},
{0.06, 54.88}, {0.08, 44.93}, {0.10, 36.76}
};
model = a b^t;


First we split the data into separate components. Here we simply take the real and imaginary parts, but you can choose whatever transformation you like as long as it represents a complete (not necessarily orthogonal) basis, and maps the complex numbers onto the reals. You should bear in mind, though, that the residuals will now be calculated in this space, and different choices will lead to different error estimates on the parameters. This may be enough (especially for a nonlinear fit) to change the optimum values of the parameters themselves, so choose your transformation wisely.

It's convenient at this point to work with the abscissae and ordinates separately. (And I should note that we assume here a one-dimensional model, which we are converting into an two-dimensional one. A similar approach can be used to change an n-dimensional complex fit into an n+1-dimensional real one, but the code will differ in minor details.)

transformation = {Re, Im};
{abscissae, ordinates} = Transpose[data];
data = Transpose[{abscissae, #}] & /@ Through@transformation[ordinates];


Now we have two data arrays rather than one, but this isn't what FindFit needs (or can accept). So, we merge these by prepending an index variable to each abscissa:

data = MapIndexed[Prepend[#1, First[#2]] &, data, {2}] ~Flatten~ 1;


Next, the model needs to be extended by a dimension. This can be done using If, Which, Piecewise and so on, but I personally think KroneckerDelta is neatest:

model = Inner[
#1[model] KroneckerDelta[i, #2] &,
transformation, Range@Length[transformation]
];


Now, FindFit gets the right answer straight away without any complaints (but don't forget to add the label for the extra dimension!):

FindFit[data, model, {a, b}, {i, t}]
(* -> {a -> 100.004, b -> 0.0000452493} *)


Obviously, if this model wasn't limited to the real line, one might wish to transform the parameters as well. That can be done, of course, but I think it's best left as a subject for another answer.

-
Useful method, thanks. – Simon Woods Aug 26 '12 at 16:31

This error can occur when the fit function is complex for the parameter values which Mathematica tries. Here's a simple example:

data = N @ Table[{x, Sqrt[100 - x]}, {x, 10}];
FindFit[data, Sqrt[a - x], {a}, x]


Mathematica evaluated the fit function with a=1 at each x value and got complex results, hence the error.

The best way to fix it is to provide starting estimates for the parameter values. The estimates often don't have to be particularly accurate, the important thing is to provide values which make the fit function real for all the x values in the data. So in our example we can choose any starting estimate of a which is greater than 10 and FindFit will work as desired:

FindFit[data, Sqrt[a - x], {{a, 10.1}}, x]


{a -> 100.}

Another possibility for this simple example is to force the fit function to be real with Re. This works here, but will probably not work well for more complicated cases.

FindFit[data, Re@Sqrt[a - x], {a}, x]


{a -> 100.}

In general you will find that the more parameters there are in the fit function, the more important it is to provide good starting estimates.

-

I considered this very regression model in this math.SE answer. Briefly stating the lessons of that answer: one might consider performing an initial linearization of $y=a b^x$ as $\log\,y=\log\,a+x\log\,b$ and then do a linear regression like usual, but as I explained in that answer, a logarithmic transformation will distort the error in your original data, such that the errors in the transformed no longer necessarily follow, say, a normal distribution, even if the error in your original data was normally distributed to begin with. (Remember that the usual assumption, often violated, for performing regression is your data having normally-distributed error.)

So, one should not settle for the results of the linearization. Here in Mathematica, due to the availability of nonlinear regression in FindFit[] and its more elaborate cousin NonlinearModelFit[], one should definitely not settle for the results of a linearization! The right way is to use FindFit[] first for the linear regression on the transformed data, to get initial estimates that can then be subsequently polished by FindFit[], but this time doing nonlinear regression:

data = {{0.0, 100.0}, {0.02, 81.87}, {0.04, 67.03},
{0.06, 54.88}, {0.08, 44.93}, {0.10, 36.76}};

{a, b} = {Exp[u], Exp[v]} /.
FindFit[Transpose[{data[[All, 1]], Log[data[[All, 2]]]}], u + v t, {u, v}, t]
{100.012, 0.0000451495}

{p, q} /. FindFit[data, p q^t, {{p, a}, {q, b}}, t]
{100.004, 0.0000452493}


Similar things can be done for all other regression models that possess a convenient linearization.

-