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 trying to learn how to use the DesignMatrix function in Mathematica and how to perform a LeastSquares fit. I am at end of my rope and seek expert intervention.

Say, I have some experimental data and some experimental fit coefficients

size = 20; (*no of points *) 
expdata = (Array[#, size] & /@ {x, y, z})~Join~{ConstantArray[1, size]};  
coeffs = {{a1, b1, c1, d1}, {e1, f1, g1, h1}, {i1, j1, k1, l1}}; (*can be larger*) 

My experimental equations would look like:

lhs = coeffs.expdata; (*Left hand side *)
rhs = Array[#, size] & /@ {u, v, w}; (*Right hand side *)
eqs = Flatten[ Table[ lhs[[#1, i]] == rhs[[#1, i]] & /@ Range[3], {i, 1, size}]];

I know I can solve these equations exactly using the Solve function as follows for a small number of data points:

Solve[ eqs, Flatten[coeffs]] 

But assuming I have a large number of noisy data points, how do I setup a LeastSquares problem. In specific, I’d like to understand how to use the DesignMatrix function….

share|improve this question
    
You've three disjoint problems in there, one for each sublist in coeffs –  belisarius Jan 28 at 20:30

1 Answer 1

up vote 8 down vote accepted

Let's say you have a set of data:

data = {{-2., -9.424}, {-1.5, -2.586}, {-1., -3.047}, {-0.5, -1.203}, \
{0., 0.551}, {0.5, 4.566}, {1., 12.077}, {1.5, 21.2118}, {2., 
    44.752}};

and you suspect it can be modeled with a 3rd order polynomial. You can set up the design matrix using the following:

a = DesignMatrix[data, {x, x^2, x^3}, x]

This link provides some useful information on how the least squares approximation works. Simply, if the equation a x = b has no solution, multiply both sides by the Transpose of a and solve the resulting equation. This can be done in Mathematica by:

soln = Solve[Transpose@a.a.{k1, k2, k3, k4} == Transpose@a.data[[All, 2]], 
  {k1,k2, k3, k4}]
(* {{k1 -> 0.368633, k2 -> 3.98313, k3 -> 4.23867, k4 -> 2.3121}} *)

Which gives us a decent fit

Plot[k1 + k2 x + k3 x^2 + k4 x^3 /. soln, {x, -2, 2},
  Epilog -> {Red, PointSize[0.02], Point /@ data}]

Mathematica graphics

which is also identical to doing it the simple way

nlm = NonlinearModelFit[data, 
  k1 + k2 x + k3 x^2 + k4 x^3, {k1, k2, k3, k4}, x]
nlm["BestFitParameters"]
(* {k1 -> 0.368633, k2 -> 3.98313, k3 -> 4.23867, k4 -> 2.3121} *)

Note I use NonlinearModelFit instead of LinearModelFit here only because the output more closely resembles the output of soln

share|improve this answer
    
What a nice, complete answer! +1 –  rasher Jan 29 at 0:37
    
Thank you. This is a much better explanation that the Help pages on Ma. –  Pam Jan 29 at 2:01

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.