9
$\begingroup$

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….

$\endgroup$
  • $\begingroup$ You've three disjoint problems in there, one for each sublist in coeffs $\endgroup$ – Dr. belisarius Jan 28 '14 at 20:30
9
$\begingroup$

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

$\endgroup$
  • $\begingroup$ What a nice, complete answer! +1 $\endgroup$ – ciao Jan 29 '14 at 0:37
  • $\begingroup$ Thank you. This is a much better explanation that the Help pages on Ma. $\endgroup$ – Pam Jan 29 '14 at 2:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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