$\chi ^2$ Fitting in Mathematica

We begin by creating arrays x_i, y_i and $\sigma _i$: 

The syntax of the Mathematica command Array is, in this context, identical to the syntax of Maxima's command array. The only difference is that whereas the array x_i created by Maxima would have entries from x₀ through x₃₀₀, in Mathematica the numbering goes from x₁ through x₃₀₀. Unfortunately we also get a screenful of gibberish after the array has been successfully allocated. Recall that in Maxima we could terminate a command with a dollar sign instead of a semicolon, and in this way tell Maxima to shut up. In Mathematica we accomplish the same task by terminating Mathematica's command with a semicolon, instead of just a newline. 

And so, we can now safely proceed to allocate arrays y_i and $\sigma _i$:

Now let us proceed with the definition of $\chi^2 = \sum_{i=1}^N \left(y_i - a - b x_i\right)^2/\sigma_i^2$, and its derivatives with respect to a and b:

Now we are in a position, finally, to formulate our equations and request their solutions:

$\ldots$ and, much to our disappointment, we get a similar answer to the one we got from Maxima. Mathematica returned an empty list, because it doesn't know how to solve these equations. And the reason for the latter is as follows:

The derivatives inside the sums remain unexpanded in the output, even though Mathematica knows how to do them. This can be easily seen by using Mathematica's equivalent of Maxima's pickapart. All Mathematica expressions are stored as lists, much as in Lisp. But they are printed on Mathematica's output in a way, to which Mathematica users are more accustomed.  You can ask Mathematica to print an expression the way it is stored by invoking a command FullForm:

You can get at individual elements of the list as follows: 

Observe that the derivative has been automatically evaluated when we have asked Mathematica to print part 1 of dchi2b[a,b], but that part is not stored in an evaluated form. We can easily create a new expression, in which the derivative is fully evaluated by calling the ReplacePart command: 

The command ReplacePart takes an expression as its first argument. The number of the expression part to be replaced is given in the third argument. The second argument contains the replacement. In our case the replacement is the expression, which Mathematica labelled as Out[32]. Recall that we had the same facility available to us in Maxima: all inputs and outputs were automatically labelled, and could be reused while conversing with Maxima.

We can now repeat this procedure for dchi2a[a,b]:

Let us attempt a solution now:

Well, we've moved forward just a little, but Mathematica's Solve still can't cope with the equations. As in Maxima we need to disentangle a and b explicitly from the sum, so that Mathematica's Solve would notice that the resulting equations are really quite simple and linear in a and b.

Let us begin by expanding the fractions in the sums:

The next step is to distribute  the sums over the expressions as follows:

These equations are still not tractable by Solve. We must factor out a and b. We do that again by manipulating the Sums manually, since tricks such as:

accomplish nothing. Only when a specific numerical substitution is made, we can eventually get somewhere:

This is too little and too late.

So here is how we rewrite these expressions manually:

$\ldots$ and similarly for the second equation:

Make no mistake, this is a butchering job that would not stand up to a scrutiny of any artificial reasoning system , such as HOL  or Isabelle . We are basically creating completely new formulas, by manually selecting parts from the old ones - and such manipulations are prone to mistakes. We had to do considerably more work in Mathematica than even in Maxima. The reason for this is that Mathematica knows even less about manipulating sums than Maxima does. Function Sum in Mathematica is designed specifically for working with power series rather than with arrays.

On the other hand, you must have noticed that Mathematica's vocabulary is rich enough to let us design more sophisticated rewriting rules for sums and derivatives, and we will explore this possibility further down the road.

Now, at long last, we can solve the equations:

We make final assignments to a and b by extracting appropriate expressions from the list returned by Solve:

As was the case in Maxima and in Maple, Mathematica evaluates $\partial a / \partial y_k$ and $\partial b/\partial y_k$ incorrectly too, returning zero in both cases:

In order to perform this differentiation, we must make specific substitutions both for N and for k, e.g.,

and similarly for $\partial b/\partial y_k$:

Compare these results with Section 2.1.8.