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$
1
  • $\begingroup$ You've three disjoint problems in there, one for each sublist in coeffs $\endgroup$ Jan 28, 2014 at 20:30

1 Answer 1

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[[email protected].{k1, k2, k3, k4} == [email protected][[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$
2
  • $\begingroup$ What a nice, complete answer! +1 $\endgroup$
    – ciao
    Jan 29, 2014 at 0:37
  • $\begingroup$ Thank you. This is a much better explanation that the Help pages on Ma. $\endgroup$
    – Pam
    Jan 29, 2014 at 2:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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