Fit model with two dependent variables

Question

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)

$$Y_1(x) = A x + B$$ $$Y_2(x) = (A / 2) x + B$$

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

$$\sum_i[(Y_{1,i} - Y_1(x_i))^2 + (Y_{2,i} - Y_2(x_i))^2] \Rightarrow \textrm{min}$$

and I also want to calculate the $R^2$ 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? .

Answer 1

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
]

Answer 2

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 *)

Answer 3

If you simply have two sets of data and no shared parameters, one should just perform two fits.

If you have two separate sets of data with shared parameters and know that the errors are independent with a single common variance, then using MultiNonlinearModelFit by @SjoerdSmit is what you want.

But you appear have "paired" response variables (or "matched" response variables when there are more than 2) which (to me) screams that the error terms are likely correlated and maybe of unequal variance. And you have common parameters. For such cases the structure of the error terms become much more important and ignoring that structure will get you wrong estimates of precision of the non-error terms. (It is true that many times the estimates of the non-error terms won't differ by much even if you ignore the error structure. But one should really want an appropriate set of estimates of precision.)

Below I describe a reasonable starting model from which one should then check on the assumptions I'm going to make explicit.

We assume that the error structure is that of a bivariate normal distribution with that same distribution for all observations and that individual observations are independent of each other. (Again, one should check on those assumptions after fitting.)

Suppose @MarcoB 's models are slightly modified:

f1[x_, a_, b_, c_] := a x^2 - b x + c
f2[x_, c_, d_] := d x^3 - c

(* Generate some data *)
x = Range[-10, 10];
parms = {a -> 1, b -> 5, c -> 3, d -> 0.1, σ1 -> 7, σ2 -> 15, ρ -> 0.5};
(* Random error from bivariate normal *)
SeedRandom[12345];
error = RandomVariate[BinormalDistribution[{0, 0}, {σ1, σ2}, ρ] /. parms, Length[x]];
(* Construct response vectors *)
response = Transpose[{f1[x, a, b, c], f2[x, c, d]} /. parms] + error ;

(* Look at the data *)
ListPlot[{Transpose[{x, response[[All, 1]]}], Transpose[{x, response[[All, 2]]}]}]

(* Get initial estimates of parameters by fitting separately and ignoring the
   correlation structure *)
(* The example functions are linear in the parameters but we'll use
   NonlinearModelFit to keep things more general *)
nlm1 = NonlinearModelFit[Transpose[{x, response[[All, 1]]}], f1[z, a, b, c], {a, b, c}, z];
nlm2 = NonlinearModelFit[Transpose[{x, response[[All, 2]]}], f2[z, c, d], {c, d}, z];
{a0, b0, c10} = {a, b, c} /. nlm1["BestFitParameters"];
σ10 = Sqrt[nlm1["EstimatedVariance"]];
{d0, c20} = {d, c} /. nlm2["BestFitParameters"];
c0 = (c10 + c20)/2;  (* Average the two separate estimates of c *)
σ20 = Sqrt[nlm2["EstimatedVariance"]];
ρ0 = Correlation[Transpose[{nlm1["FitResiduals"], nlm2["FitResiduals"]}]][[1, 2]];

(* Define the log of the likelihood *)
logL = LogLikelihood[BinormalDistribution[{0, 0}, {σ1, σ2}, ρ], 
   response - Transpose[{f1[x, a, b, c], f2[x, c, d]}]];

(* Find maximum likelihood estimates *)
mle = FindMaximum[{logL, σ1 > 0 && σ2 > 0 && -1 < ρ < 1}, 
  {{a, a0}, {b, b0}, {c, c0}, {d, d0}, {σ1, σ10}, {σ2, σ20}, {ρ, ρ0}}];

(* Estimate parameter covariance matrix and associated standard errors *)
cov = -Inverse[(D[logL, {{a, b, c, d, σ1, σ2, ρ}, 2}]) /. mle[[2]]];
se = Sqrt[Diagonal[cov]] ;(* Standard errors *)

(* Create a parameter table *)
t = ({a, b, c, d, σ1, σ2, ρ} /. mle[[2]])/se;
pvalue = 2 (1 - CDF[NormalDistribution[0, 1], Abs[t]]);
TableForm[Transpose[{{a, b, c, d, σ1, σ2, ρ} /. mle[[2]], se, t, pvalue}],
 TableHeadings -> {{"a", "b", "c", "d", "σ1", "σ2", "ρ"}, 
 {"Estimate", "Standard Error", "z-Statistic", "P-value"}}]

We can compare this against the parameter tables of the separate fits:

nlm1["ParameterTable"]

nlm2["ParameterTable"]

We see that the standard error for $c$ is about half of that of the separate fits.

Here is the data and the fit:

Show[ListPlot[{Transpose[{x, response[[All, 1]]}], Transpose[{x, response[[All, 2]]}]}],
 Plot[{f1[z, a, b, c] /. mle[[2]], f2[z, c, d] /. mle[[2]]}, {z, Min[x], Max[x]}]]

With this approach one can also fit relatively simple mixed models.