# Fit model with two dependent variables

Say I have an experiment that produces two values, for a given input. For example you get a complex response value.

Now I have a model that - of course - also produces two dependent values for one independent value, and some other parameters and I want to find the set of parameters that describe this model best (least-squares).

I basically need two equations for this model then (simple example, sorry, I don't know TeX)

Y1(x) = A * x + B
Y2(x) = (A / 2) * x + B


I want to find A and B, such that the overall sum of squares is minimal

Sum_over_i[(Y1_i - Y1(x_i))^2 + (Y2_i - Y2(x_i))^2] ==> min


and I also want to calculate the R² value for this fit.

How can I do this in Mathematica? I tried using the NonlinearModelFit function like this (data is in the form {X, Y1, Y2}):

data = {{1, 2, 1}, {2, 6, 7}, {3, 8, 9}, {4, 6, 5}, {5, 10, 12}, {6,
23, 18}};
NonlinearModelFit[data, {A*x + B, (A/2)*x + B}, {A, B}, x]


However I get the (not completely unexpected) error

Number of coordinates (2) is not equal to the number of variables (1).


The result I got from another program is A = 2.70824951933389 and B = 1.80751165636783 with R² = 0.52.

I just tried supplying two model functions instead of one. I couldn't find any examples for my case either, since you mostly find multiple-regression problems, with more than one independent variable instead, like in this question How can I use nonlinear model fit with multiple variables? .

• How is the correlation coefficient defined in that "another program"? (And what "another program" are you using, BTW?) – Dr. belisarius Jul 10 '15 at 17:00
• @belisarius I used Origin. From what I found by calculating it by hand they calculate the coefficient by dividing the sum of squared residuals of the fit by the sum of squared differences of each Y point (both Y1 and Y2) from the overall Y-mean (over all Y1 and Y2 values). I don't know if that is the correct approach though and I want to cross-check the method. Please see my other question on stats-overflow stats.stackexchange.com/questions/160326/… – Jens Jul 10 '15 at 17:07
• For complex data specifically, I wrote a package to do the process explained in the answers automatically. – Oleksandr R. Jul 10 '15 at 20:55
• @OleksandrR. Nice work! I will give it a try – Jens Jul 11 '15 at 6:07

There is no built-in way to accomplish what you want with NonlinearModelFit, or with the other fitting functions, unfortunately. Of course, you could write your own target function and NMinimize that, as belisarius has mentioned already in his answer. This is very general and not particularly difficult to do, but it does not give you access to the wealth of built-in fit description parameters that are accessible through NonlinearModelFit. Again, you could recalculate those manually, but it would be useful to be able to leverage the existing functionality instead.

In order to "trick" NonlinearModelFit into fitting your multivariate model, you will want to build an object function that includes both models and can "switch between them", based on a dummy selector variable that you will introduce in your data set and treat as an independent variable.

Let's go through a simple example. I generated some noisy data as a function of a common $x$ range which is the actual independent variable in your physical model. The two functions I used to generate the dependent variables are $y_1=x^2-5x+3$ and $y_2=0.1\ x^3-3$. I then added some random noise to the $y$ variables just to make this slightly more realistic. Here is the data set:

originaldata = {{-10, 168, -83}, {-9, 150.1, -86.9}, {-8, 129.6, -36.2}, {-7, 101.5, -34.3}, {-6, 74.8, -14.6}, {-5, 52.5, -32.5}, {-4, 48.6, -12.4}, {-3, 21.1, -23.7}, {-2, 22, -7.8}, {-1, 15.3, 5.9}, {0, 5, 8}, {1, -3.9, -7.9}, {2, -11.4, -7.2}, {3, -11.5, -11.3}, {4, 10.8, 22.4}, {5, 4.5, 8.5}, {6, 14.6, 29.6}, {7, 28.1, 31.3}, {8, 33, 31.2}, {9, 61.3, 68.9}, {10, 76, 96}};

ListPlot[originaldata[[All, {1, #}]] & /@ {2, 3}, PlotRange -> All]


You will have to reorganize this data set so that each point has the following form: {xValue, datasetIndentifier, yValue}, where we will use the numbers $1$ or $2$ as dataset indentifiers.

Transpose@originaldata;
Insert[%[[{1, #}]], ConstantArray[# - 1, Length[originaldata]], 2] & /@ {2, 3};
reorganizeddata = Flatten[Transpose /@ %, 1]


{{-10, 1, 153}, {-9, 1, 131}, {-8, 1, 107}, {-7, 1, 80}, {-6, 1, 59}, {-5, 1, 59}, {-4, 1, 31}, {-3, 1, 27}, {-2, 1, 21}, {-1, 1, 14}, {0, 1, 5}, {1, 1, 9}, {2, 1, -9}, {3, 1, -10}, {4, 1, 7}, {5, 1, 2}, {6, 1, 8}, {7, 1, 25}, {8, 1, 26}, {9, 1, 33}, {10, 1, 45}, {-10, 2, 139.36}, {-9, 2, 79.6883}, {-8, 2, 0.47189}, {-7, 2, -35.8288}, {-6, 2, -31.6333}, {-5, 2, -47.4649}, {-4, 2, -15.9864}, {-3, 2, 3.2233}, {-2, 2, -4.73885}, {-1, 2, -3.77884}, {0, 2, -20}, {1, 2, -0.221163}, {2, 2, -13.2612}, {3, 2, 13.7767}, {4, 2, 24.9864}, {5, 2, 52.4649}, {6, 2, 54.6333}, {7, 2, 12.8288}, {8, 2, -6.47189}, {9, 2, -80.6883}, {10, 2, -154.36}}

Now let's create an appropriate model function.

Clear[modelfunction]
modelfunction[indepvar_?NumericQ, datasetselector_?NumericQ, p11_?NumericQ, p12_?NumericQ, p13_?NumericQ, p21_?NumericQ, p22_?NumericQ] :=
Piecewise[{
{p11 indepvar^2 + p12 indepvar + p13,
datasetselector == 1},
{p21 indepvar^3 + p22,
datasetselector == 2}
}]


... and use it with NonlinearModelFit:

nlm = NonlinearModelFit[
reorganizeddata,
modelfunction[x, selector, a, b, c, d, e],
{a, b, c, d, e},
{x, selector}
]


We can then very conveniently extract any of the fit descriptors available through NonlinearModelFit:

nlm["ParameterTable"]
nlm["RSquared"]


Of course you can also plot the experimental points and the fit:

Plot[
nlm["BestFit"] /. selector -> # & /@ {1, 2}, {x, -10, 10},
Epilog -> {PointSize[0.015], Gray, Point[reorganizeddata[[All, {1, 3}]]]},
Evaluated -> True
]


• That is a very nice way (goes in line with the answer to a different question i found a few minutes ago, mentioned in my comment to belisarius answer). What I find curious, is that the R-Squared value I get from the nlm result is a lot larger than what Origin gives me (0.84 for my example data, as compared to 0.51 in Origin). Is the value, that comes out of the nlm correct, or is it distorted somehow by the workaround? – Jens Jul 10 '15 at 18:06
• I found the difference in R²: groups.google.com/forum/#!searchin/… Seems like it's a matter of definition and Mathematica actually has a more general description of R² for nonlinear models. How are statisticians not going completely insane after 3 days? – Jens Jul 10 '15 at 18:28
• @Jens That's an interesting point! I can't seem to open the link you provided though: it leads to a banned content warning page from Google. Do you have another way to access the content? Also, yeah, I have come to the conclusion that statisticians have a higher pain tolerance than I do, or otherwise they enjoy the pain they inflict upon others that try to use statistics... :-) – MarcoB Jul 10 '15 at 18:37
• goo.gl/eEdcbJ Let's try that again, I couldn't either, even though I copy pasted the link. Ok that seems to work. Stackexchange added some part's to the URL for whatever reason that caused it to not work. – Jens Jul 10 '15 at 18:38
• @Jens Your most recent link works for me. That's indeed an interesting point about Rsquared. – MarcoB Jul 10 '15 at 18:57

If you know the form you want to minimize, Mathematica provides several ways to do that. For example:

data = {{1, 2, 1}, {2, 6, 7}, {3, 8, 9}, {4, 6, 5}, {5, 10, 12}, {6, 23, 18}};
y1[A_, B_, x_] := A*x + B
y2[A_, B_, x_] := (A/2)*x + B
s = Sum[(i[[2]] - y1[A, B, i[[1]]])^2 + (i[[3]] - y2[A, B, i[[1]]])^2, {i, data}];
NMinimize[s, {A, B}]
(* {211.085, {A -> 2.70825, B -> 1.80751}} *)


Which seems what Origin is doing

Edit

You may try to define something like a Pearson's coefficient:

m = Mean /@ Transpose@data;
mx = First@m;
my = m[[2 ;; 3]];
rr = {Sum[(i[[1]] - mx) (i[[2 ;; 3]] - my), {i, data}],
Sqrt[(Sum[(i[[1]] - mx)^2, {i, data}] Sum[(i[[2 ;; 3]] - my)^2, {i, data}])]};
Divide @@ (Norm /@ rr) // N
(* 0.855367 *)

• Thank you for the answer. I just found a different question about the same topic here: mathematica.stackexchange.com/questions/15905/… Getting the fitted values works well in both ways (yours and their answers). However the R² value I get from the NonlinearModelFit when using @b.gatessucks answer is different than the Origin answer. Can you show me a way to calculate R² using your approach? – Jens Jul 10 '15 at 17:53
• @Jens See edit, please – Dr. belisarius Jul 10 '15 at 18:17
• That's a good idea! Due to the additional benefits of the NonlinearModelFit results I think I will use that one, as explained in the other answer, though. As you have pointed out correctly, there seem to be quite a lot of ways to do this. – Jens Jul 10 '15 at 18:31
• @Jens I believe the other way has its shortcomings too.Just think about a case where due to multiplication factors the y2 series is orders of magnitude larger than the other. – Dr. belisarius Jul 10 '15 at 18:34
• The error function to minimize should be the same in both approaches, isn't it? – Jens Jul 10 '15 at 18:37