Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am CS major, taking Computational Numerical Analysis course. Instructor gave us freedom of choice, we were allowed to use anything or any computer language we picked, I picked Mathematica.This is my 3rd day with the program.

I am writing a simple procedure to calculate the Least-Square $m^\text{th}$ Degree Polynomials based on data sample. Algorithm is somewhat simple: http://www.efunda.com/math/leastsquares/lstsqrmdcurve.cfm
This is how far I got:

y_sample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
x_sample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}(*Some sample data*)

f[x_] = Sum[Subscript[a, i] x^i, {i, 0, 4}](*generating polynomial, in this case, degree of 4*)    
Sum[(y_sample[[i]] - f[x_sample[[i]]])^2, {i, 0, 4}] 

Last line of code should roughly correspond to

$$ \Pi = \sum^n_{i = 1}\left[y_i - f(x_i)\right]^2 = \sum^n_{i=1}\left[y_i- (a_0 + a_1 x_i + a_2 x_i^2 + \ldots + a_m x_i^m)\right]^2 = \text{min} $$

List of problems and questions, I am facing:

  1. Executing code above, produces number of errors that I don't mean much to me.
  2. How do I differentiate in respect to all $\mathbb{a}_n$
  3. Differentiation should produce number of linear equations, how do I set the results of differentiation to zero and solve for unknown coefficients.
  4. Outputs are enormously long, is there a way to suppress some outputs !
share|improve this question
6  
Please DON'T USE underscoring in vriable's names! –  belisarius Jun 11 '12 at 15:49
    
Why ? Does it have something to do with how Mathematica interprets things ? –  newprint Jun 11 '12 at 15:54
    
@user465292 Yes, see e.g. Blank, BlankSequence, BlankNullSequence! –  Artes Jun 11 '12 at 15:56
8  
I really recommend you to grab some basic language totorial and go through it, will save you a lot of trouble. Have a look at this discussion for a comprehensive list of available resources. –  Leonid Shifrin Jun 11 '12 at 16:12

3 Answers 3

up vote 12 down vote accepted

You have some errors in your syntax:

  • you name your lists x_sample and y_sample, but in Mathematica, an underscore is not allowed in names (as it is reserved to patterns).
  • your last sum runs from 0, but in Mathematica, the first element in a List has index 1
  • your last sum should run until the number of data points, not 4
  • furthermore, I would advise you to use SetDelayed(:=) in your function definition, as it prevents possible naming conflicts
  • also in your function definition, letting the sum run from 0 will give problems for x=0, as it produces 0^0, which is Indeterminate.

This gives:

ysample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
xsample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}

f[x_] := Subscript[a,0] + Sum[Subscript[a, i] x^i, {i, 1, 4}]
result = Sum[(ysample[[i]] - f[xsample[[i]]])^2, {i, 1, Length[ysample]}]

Next, to calculate the derivative with respect to, say, Subscript[a,0], use D[result,Subscript[a,0]].

Last, to get a solution to the set of linear equations, use Solve:

sols = Solve[Table[D[result, Subscript[a, i]] == 0, {i, 0, 4}], Table[Subscript[a, i], {i, 0, 4}]]

To put the values of the coefficients in your polynomial, use Replace(/.):

final[x_] := (f[x] /. sols)[[1]]
final[x]
final[2]

EDIT

As you can see in the result, there are some values which are almost zero (10^-11 or smaller). These are simply the result of some internal round-off errors, and should of course be equal to zero. You can easily add a piece of code to discard all coefficients which are a factor smaller than the largest coefficient.

The largest coefficient is given by:

M = Max[Apply[#2 &, sols, {2}][[1]]

(don't mind the syntax too much yet, it is just a way to select the values from sols in order to be able to use Maxon it)

Next we make the replacement for every x which is at least 10^6 times smaller than the maximum:

sols = sols/. x_ /; x < M / 10^6 :> 0
final[x_] := (f[x] /. sols)[[1]]
final[x]
final[2]

EDIT 2

It appears Mathematica has a built-in function to deal with near-zero numerics: Chop. It's al lot more easy to use than the method I proposed in my former edit; just write Chop[sols] and it will automatically thread over the rules. Thanks J.M. for pointing that out!

Then the final code would be:

ysample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
xsample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}

f[x_] := Subscript[a,0] + Sum[Subscript[a, i] x^i, {i, 1, 4}]
result = Sum[(ysample[[i]] - f[xsample[[i]]])^2, {i, 1, Length[ysample]}]

sols = Solve[Table[D[result, Subscript[a, i]] == 0, {i, 0, 4}], Table[Subscript[a, i], {i, 0, 4}]]

final[x_] := (f[x] /. Chop[sols])[[1]]
final[x]
final[2]
share|improve this answer
    
Could you elaborate on the last 3 lines of code, it looks like definition is recursive. Thanks ! –  newprint Jun 11 '12 at 16:43
    
The first line, final[x_] := (f[x] /. sols)[[1]], just defines the final function to be the function f (the generating polynomial) with its coefficients replaced by the solutions. I added the [[1]] to get rid of the extra {}, because the result of this replacement is a List. The two other lines are just examples; final[x] gives the polynomial in function of x and final[2] gives the polynomial for x=2 (the result is 20, as the function is y=10 x). –  freddieknets Jun 11 '12 at 16:59
    
Thank you for taking time to write this extensive answer! –  newprint Jun 11 '12 at 19:28
2  
"You can easily add a piece of code to discard all coefficients which are a factor smaller than the largest coefficient." - I'd have told OP about Chop[]... –  J. M. Jun 12 '12 at 8:24
    
Wow, thanks! I didn't know this function (I'm not too much used to work with floating points). –  freddieknets Jun 12 '12 at 9:42

freddieknets already told you what you're supposed to be doing. Here's how I'd have done your approach to least squares:

n = 4;
f[x_] = Map[C, Range[0, n]].x^Range[0, n];


Solve[Thread[
    D[Total[(ysample - f /@ xsample)^2], {Map[C, Range[0, n]]}] == 0],
    Map[C, Range[0, n]]] // First // Chop
{C[0] -> 0, C[1] -> 10., C[2] -> 0, C[3] -> 0, C[4] -> 0}


f[x] /. %
10. x

Notes:

  • C was always intended for indexed "arbitrary constant" use; this is one case where it certainly is useful.

  • The dot product is useful for assembling a polynomial from its coefficients. An alternative is to use a "Horner approach": f[x_] = Expand[Fold[(#1 x + #2) &, 0, Map[C, Range[n, 0, -1]]]]

  • The least-squares objective function is compactly implemented here as Total[(ysample - f /@ xsample)^2]; Sum[] is useful notationally, but in Mathematica, there are always other methods.

  • Remember that in an optimization such as this, one is supposed to be solving for the zeroes of the gradient; in this case, note that D[f[x, y, (* other variables *)], {{x, y, (* other variables *)}}] is the way to quickly generate a gradient.

  • Thread[list == 0] equates all the components of list to 0.

  • As I told freddieknets, Chop[] is very useful for problems such as this, when you want to treat tiny terms as actually being zero.


For giggles, here's one way I'd have done, if none of Fit[], FindFit[], LeastSquares[], or PseudoInverse[] were available:

n = 4;
vm = Take[LinearAlgebra`VandermondeMatrix[xsample], n + 1];
{qm, rm} = QRDecomposition[Transpose[vm]];
LinearSolve[rm, qm.ysample].x^Range[0, n] // Chop

I could have used CholeskyDecomposition[] (the normal equations route), but QRDecomposition[] is often more numerically sound for applications like these.

There's also SingularValueDecomposition[], but that requires slightly more machinery to discuss properly.

Finally, I gave a rundown of most of the usual methods for linear least squares in this stats.SE answer; you might be interested in it.

share|improve this answer

Also, remember you can check your results with:

ysample = {0, 1, 2, 3, 4, 5, 6, 7, 8}
xsample = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8}(*Some sample data*)
model = a x^3 + b x^2 + c x + d;
fit = FindFit[data = Transpose@{xsample, ysample}, model, {a, b, c, d}, x]
modelf = Function[{x}, Evaluate[model /. fit]]
Plot[modelf[x], {x, 0, 1}, Epilog -> Map[Point, data]]
share|improve this answer
2  
Also of interest is the LeastSquares[] function. –  J. M. Jun 11 '12 at 17:52
    
@J.M. ... or PseudoInverse –  Sjoerd C. de Vries Jun 11 '12 at 19:41

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.