3
$\begingroup$

I have some composition data for a geological sample and its constituent minerals. I want to estimate the proportions of the minerals that make up the bulk sample. This is essentially a mass balance problem. Each variable is pseudo-independent (but constrained by constant sum).

I'm starting with the matrix form $a.x=b$, where, $a$ is the array $elemental~composition~vector \times mineral ~species$, $x$ is the column vector of mineral proportions, and $b$ is the bulk composition. The key constraint with this kind of problem is that

$0 \leq x_{j} \leq 1~~~~~(j=1,...m)$

and

$\sum_{j=1}^m x_{j}=1$

here is some example data

a={{63.1545, 64.3049, 100., 37.1417, 32.4026, 30.0382, 30.7033, 0., 36.5444, 0., 0., 0., 0.0034},
   {0.0016, 0.01, 0., 0.6641, 2.35946, 26.91, 0., 0., 0.297125, 0., 0., 0., 51.3026}, 
   {17.8683, 21.096, 0., 11.2387, 14.3019, 6.46115, 0.34805, 0., 14.368, 0., 0., 0., 0.}, 
   {0.0223, 0.1191, 0., 30.6859, 32.0779, 1.5471, 1.5014, 0., 11.5299, 0., 0., 0., 43.3621}, 
   {0., 0., 0., 0.94408, 0.6131, 0.0986, 0., 0., 0., 0., 0., 0., 2.1039}, 
   {0., 0., 0., 1.14598, 1.82486, 0.03695, 0.00935, 0., 0.2312, 0., 0., 0., 0.0166}, 
   {0.625667, 9.77073, 0., 1.35584, 0.0558, 0.05635, 0., 0., 0.079275, 0., 0., 0., 0.018}, 
   {0., 2.80603, 0., 10.7388, 0.0245429, 27.7454, 2.19443, 55.08, 10.5478, 56.03, 78.13, 10.44, 0.0262}, 
   {15.7365, 0.0874667, 0., 1.89652, 8.54851, 0.00475, 0., 0., 0.09665, 0., 0., 0., 0.0115}, 
   {0., 0., 0., 0.16636, 0., 0., 0.0293, 42.4, 0.0114, 0., 0., 0., 0.0617}}

b={65.67, 0.52, 14.77, 5.418, 0.13, 0.3, 3.22, 2.05, 6.07, 0.13}

If I use the LeastSquares function I get

lsq = LeastSquares[a, b]

{0.331681, 0.283916, 0.166439, 0.172393, 0.0586919, 0.0795275, 0.00802007, 0.00244529, -0.030593, -0.0157832, -0.0220087, -0.00294087, -0.0363984}

This result looks OK, except negative proportions are not allowed and the sum does not equal 1.

Total@lsq

0.99539

How do I implement linear least squares fitting for multivariate data under the constraint that values lie between 0 to 1 and sum to 1 in Mathematica?

Also how can I derive the residuals for the fit?

$\endgroup$
  • 2
    $\begingroup$ My usual soapbox speech: Why not consult a statistician first and then figure out how to implement an appropriate process in Mathematica? Consider stats.stackexchange.com/questions/267014/…. $\endgroup$ – JimB Mar 20 at 14:53
  • $\begingroup$ @JimB precisely because someone will tell me that I can't use continuum methods with compositional data. The constant sum problem is considered a vague oddity in geological circles and largely ignored (for better or worse). $\endgroup$ – geordie Mar 20 at 15:26
  • $\begingroup$ I'm sure there must be more to it than "They won't let me do what I want!" Compositional data analysis has a pretty standard set of statistical approaches. Chemists use it all the time. These approaches are many times called "mixture designs": itl.nist.gov/div898/handbook/pri/section5/pri54.htm. Maybe you've been talking to the wrong folks. $\endgroup$ – JimB Mar 20 at 15:49
  • $\begingroup$ Like I said, geologists don't. Perhaps it has to do with the general 'noisiness' of geoscience datasets. Most other sciences are able to rely on fairly precise numbers. $\endgroup$ – geordie Mar 20 at 16:22
  • 1
    $\begingroup$ "Most other sciences are able to rely on fairly precise numbers." Ha, ha. You're funny! $\endgroup$ – JimB Mar 20 at 16:30
6
$\begingroup$

One approach is to reframe it as a minimization problem:

xVec = Array[x, 13];
NMinimize[{Total[(a.xVec - b)^2], Total[xVec] == 1, Thread[xVec >= 0]}, xVec]

You can also get an improved residual by relaxing the equality constraint:

NMinimize[{Total[(a.xVec - b)^2] + (Total[xVec] - 1)^2, Thread[xVec >= 0]}, xVec]

Thanks to Roman for some simplifications.

$\endgroup$
  • $\begingroup$ Thread[xVec > 0] would be simpler on the third line and does the same thing. Also, you can constrain the sum by equality: NMinimize[{(a.xVec-b).(a.xVec-b), Total[xVec]==1, Thread[xVec>0]}, xVec]. I think the square in Total[(a.xVec - b)]^2 should be taken inside of the Total, not outside. $\endgroup$ – Roman Mar 20 at 14:15
4
$\begingroup$

This can be done as a linear programming problem although the optimization is not least-squares in that case. The idea is to set up new variables, which are absolute values of discrepancies between a.x. and b. This can be done as below.

xvec = Array[x, Length[a[[1]]]];
dvec = Array[d, Length[a]];
lpolys = a.xvec - b;
ineqs = Flatten[{Thread[xvec >= 0], Thread[dvec - lpolys >= 0], 
    Thread[dvec + lpolys >= 0], Total[xvec] == 1}];

Now use NMinimize.

{min, vals} = 
 NMinimize[{Total[dvec], ineqs}, Join[xvec, dvec], PrecisionGoal -> 10]
Total[xvec /. vals]

(* Out[90]= {0.627674050324, {x[1] -> 0.339337833732, 
  x[2] -> 0.292918812222, x[3] -> 0.195908896153, 
  x[4] -> 0.10491683141, x[5] -> 0.0591148361052, x[6] -> 0., 
  x[7] -> 0., x[8] -> 0.00181140763364, x[9] -> 0., x[10] -> 0., 
  x[11] -> 0., x[12] -> 0., x[13] -> 0.00599138274489, d[1] -> 0., 
  d[2] -> 0., d[3] -> 0.502611255236, d[4] -> 0., 
  d[5] -> 0.0178984583702, d[6] -> 0.0717916527728, d[7] -> 0., 
  d[8] -> 0., d[9] -> 0., d[10] -> 0.0353726839451}}

Out[91]= 1. *)
$\endgroup$

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.