$\chi ^2$ Fitting in Maxima

Within Maxima we have to declare three arrays: x_i, y_i, and $\sigma _i$: 

Here we have used the command array, to declare simultaneously three arrays , whose names are given in a list, which has been inserted into the first slot of the command. The arrays are all of the same shape and size, 300 entries long and are numbered from 0 through 299, although we won't make any use of the latter.

Now let us define the figure of merit, $\chi ^2$ :

We locate the minimum of $\chi ^2$ by equating its partial derivatives with respect to a and b to zero, and solving the resulting two equations for a and b:

So far so good. We have just evaluated that

$$\frac{\partial\chi^2} {\partial a} = - 2 \sum_{i=1}^N \frac{... ...m_{i=1}^N \frac{x_i\left(y_i - b x_i - a\right)} {\sigma_i^2}$$

In principle all that we need to do now, should be to invoke solve as follows:

But solve responds with an empty list, because it cannot untangle a and b from the sums and solve this simple linear system of equations. Maxima's solve, sadly, is not very clever. The above, as you will see soon, will work without a glitch in Maple and in Mathematica.

In Maxima we have to help it a little. We do that by extracting a and b from the sum manually:

Maxima gets infuriatingly stubborn at this stage and refuses to split sums eda and edb into three separate sums each, which is what we need to do in order to persuade solve that the equations are easily tractable.

So, at this stage we tell Maxima what to do ourselves. The command pickapart(eda, 4) splits the whole expression , eda into its components down to the fourth nesting level. The components are referred to as E21, E22, and E23. Using these we assemble a new expression manually by saying:

eeda : sum(e21, i, 1, n) + sum(e22, i, 1, n) + sum(e23, i, 1, n);

Here is how it all happens:

and now similarly for edb:

At this stage a and b have been taken out of the sums explicitly, and the resulting equations are easily seen to be linear in aand b. Consequently solve can be invoked again, this time on a simplified problem:

Translating this into a normal mathematical notation we obtain the following results:

$$a = - \frac{\sum_{i=1}^N\frac{x_i^2}{\sigma_i^2} \sum_{i=1}^... ...}^N\frac{1}{\sigma_i^2} \sum_{i=1}^N\frac{x_i^2}{\sigma_i^2}}$$

and

$$b = \frac{\sum_{i=1}^N\frac{x_i}{\sigma_i^2} \sum_{i=1}^N\fr... ...}^N\frac{1}{\sigma_i^2} \sum_{i=1}^N\frac{x_i^2}{\sigma_i^2}}$$

But our job is still not finished. We have to find expressions for standard deviations in a and b. We extract the solution from the list of lists returned by solve, again, by using pickapart:

This way we can make new assignments to a and b:

Placing a dollar  instead of a semicolon at the end of an expression supresses Maxima chatter mode, so that we don't see the same large expressions being printed all over again.

a and b so defined are functions of y_i. Errors in y_i, the sigmas, $\sigma _i$, propagate into errors in a and in b  following the error propagation law:

$$\sigma_a = \sqrt{\sum_{k=1}^N \sigma_k^2 \left(\frac{\partial a}{\partial y_k}\right)^2}$$

and

$$\sigma_b = \sqrt{\sum_{k=1}^N \sigma_k^2 \left(\frac{\partial b}{\partial y_k}\right)^2}$$

But here we are going to hit a problem. The problem is that both Maxima and Maple will evaluate quite incorrectly the expression

$$\frac{\partial a}{\partial y_k}$$

namely:

This is a bug! In this context all three, Maxima, Maple, and Mathematica, think that y_k is an altogether different beast from y_i, and since they don't find any explicit occurrence of y_k in a they return $\partial a/\partial y_k = 0$.

But we can work around it.

Consider a simpler expression:

$$\sum_{i=1}^N x_i y_i$$

and its derivative with respect to say, y₃:

$$\frac{\partial}{\partial y_3}\sum_{i=1}^N x_i y_i = \frac{\p... ...( x_1 y_1 + x_2 y_2 + x_3 y_3 + x_4 y_4 + \ldots \right) = x_3$$

Now let us try to reproduce this result in Maxima.

Observe that Maxima is reluctant to evaluate the derivative, even when it has been told about the values of N and k (until then Maxima would be quite correct to say that it could be zero, in case k > N). Only when we have forced Maxima to expand the sum into a finite expression, by invoking ev with the simpsum flag, did the differentiation finally proceed and the correct result, x₃, was returned.

So, let us use this trick to evaluate $\partial a / \partial y_k$ for some specific value of k, say, 3, in Maxima:

This looks quite good, but it is a touch too long to digest. We need to rewrite expression D70 in a more compact form, using again capital sigmas for sums. We use our faithful pickapart to do that:

I have performed this calculation earlier (as you can see from the C number that jumped suddenly to 200), which is why Maxima does not print the values of E172, E150, and so on. But you can always ask Maxima what they are. For example:

To save you the effort:

$$\begin{eqnarray*}\frac{1}{\texttt{E150}} &=& \frac{1}{\sigma_3^2} \\ \texttt{E1... ... S_{xx} \sum_{i=1}^5\frac{1}{\sigma_i^2} = S_{xx} S_{\sigma} \\ \end{eqnarray*}$$

Making the following substitutions: $5 \rightarrow N$, and $3 \rightarrow k$, we can finally write down the following result:

$$\frac{\partial a}{\partial y_k} = - \frac{\frac{1}{\sigma_k... ...sigma_k^2 \left(S_{xx} S_{\sigma} - \left(S_x\right)^2\right)}$$

We will now apply the same methodology to evaluate $\partial b/\partial y_k$.

where, as before,

$$\begin{eqnarray*}\texttt{E156} &=& S_x \\ \texttt{E174} &=& S_{xx} S_\sigma \\ ... ...a_3^2} \\ \texttt{E215} &=& \frac{x_3}{\sigma_3^2} S_\sigma \\ \end{eqnarray*}$$

Again replacing $5 \rightarrow N$ and $3 \rightarrow k$ we obtain the following result:

$$\frac{\partial b}{\partial y_k} = \frac{\frac{1}{\sigma_k^2}... ...\sigma_k^2 \left(S_{xx} S_\sigma - \left(S_x\right)^2 \right)}$$

Those quantities, $S_\sigma$, S_xx, etc., are clearly very useful in $\chi ^2$-fitting. By introducing S_y and S_xywe can rewrite expressions for a and b in a more compact form:

$$\begin{eqnarray*}a &=& - \frac{\sum_{i=1}^N\frac{x_i^2}{\sigma_i^2} \sum_{i=1}^... ...ma S_{xy} - S_x S_y} {S_\sigma S_{xx} - \left(S_x\right)^2} \\ \end{eqnarray*}$$