3
$\begingroup$

I've got some data with 2 independent variables and 1 dependent variable I'm trying to fit a curve to. I'm really kind of a dummy with Mathematica except for the most basic things, so I looked through the documentation of FindFit and to be honest I'm a little lost.

My data is:

X = {50, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 62, 63, 63, 
 64, 65, 66, 67, 67, 68, 69, 69, 70, 70, 70, 71, 72, 72, 72, 73, 73, 
 73, 73, 73, 73, 74, 74, 74, 75, 76, 76, 76, 76, 76, 76, 76, 76}

Y = {19.9, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 30.4, 50, 
 24.3, 50, 50, 50, 50, 47.1, 50, 50, 32.4, 50, 51.1, 45.3, 50, 50, 51, 
 106.5, 50, 39.7, 71.4, 80.2, 89.2, 160.6, 50, 148.9, 36, 50, 50, 
 71.4, 97.2, 116.8, 207.6, 224.2, 230.1, 232.6, 50}

Z = {0.199, 0.500, 0.495, 0.490, 0.485, 0.481, 0.476, 0.472, 0.467, 
 0.463, 0.459, 0.455, 0.450, 0.272, 0.446, 0.216, 0.442, 0.439, 0.435, 
 0.431, 0.402, 0.427, 0.424, 0.272, 0.420, 0.425, 0.377, 0.417, 0.413, 
 0.418, 0.874, 0.410, 0.322, 0.580, 0.653, 0.725, 1.306, 0.407, 1.202, 
 0.291, 0.403, 0.400, 0.567, 0.771, 0.926, 1.647, 1.780, 1.826, 1.846, 0.397}

I'm looking for a function f such that f(X,Y)=Z.

I know I have to give it a model, but I'm not even sure what a good model to try would be. I know that for constant X, f is almost linear in Y, so maybe having f be quadratic in terms of Y would be good. I'm not sure about f's relationship to X, but I suspect it's a similar situation so quadratic in terms of X maybe as well.

How do I set this up and find a good fit?

Improve this question
$\endgroup$
3
  • $\begingroup$ Sorry about that, it was a partial copy/paste of the X values. They should each be Length 50 now $\endgroup$
    – Zachary Turner
    Commented Dec 30, 2015 at 2:28
  • 2
    $\begingroup$ If I want to find a model I look at a plot first. Graphics3D[Line[Transpose[{X, Y/8, Z*12}]]] I've scaled Y and Z to make the plot approximately a cube. Putting my mouse inside that cube, pressing the left button and dragging seems to hint that roughly x,y,z lie in a plane and I might start looking for a model there. $\endgroup$
    – Bill
    Commented Dec 30, 2015 at 3:11
  • 2
    $\begingroup$ What sort of model would be suited to your data is not a Mathematica issue. It depends on what real world phenomena the data represents and knowing what models experts in that phenomena have developed. $\endgroup$
    – m_goldberg
    Commented Dec 30, 2015 at 5:42

2 Answers 2

Reset to default
10
$\begingroup$

As @Bill points out, plotting the data is essential. (Also "multivariate" is usually reserved for multiple responses rather than multiple predictor variables and what you want to do is called "multiple regression".) Once you have some idea as to what form the candidate models might take, you can fit those using several fitting functions:

x = {50, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 62, 63, 
   63, 64, 65, 66, 67, 67, 68, 69, 69, 70, 70, 70, 71, 72, 72, 72, 73,
    73, 73, 73, 73, 73, 74, 74, 74, 75, 76, 76, 76, 76, 76, 76, 76, - 76};

y = {19.9, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 30.4, 50, 
   24.3, 50, 50, 50, 50, 47.1, 50, 50, 32.4, 50, 51.1, 45.3, 50, 50, 
   51, 106.5, 50, 39.7, 71.4, 80.2, 89.2, 160.6, 50, 148.9, 36, 50, 
   50, 71.4, 97.2, 116.8, 207.6, 224.2, 230.1, 232.6, 50};

z = {0.199, 0.500, 0.495, 0.490, 0.485, 0.481, 0.476, 0.472, 0.467, 
   0.463, 0.459, 0.455, 0.450, 0.272, 0.446, 0.216, 0.442, 0.439, 
   0.435, 0.431, 0.402, 0.427, 0.424, 0.272, 0.420, 0.425, 0.377, 
   0.417, 0.413, 0.418, 0.874, 0.410, 0.322, 0.580, 0.653, 0.725, 
   1.306, 0.407, 1.202, 0.291, 0.403, 0.400, 0.567, 0.771, 0.926, 
   1.647, 1.780, 1.826, 1.846, 0.397};

data = Transpose[{x, y, z}];

FindFit[data, a + b xx + c yy, {a, b, c}, {xx, yy}]
(* {a -> 0.26709594179703233,b -> -0.0035772589306014156,c -> 0.00799624708810228} *)

LinearModelFit[data, {xx, yy}, {xx, yy}]["BestFitParameters"]
(* {0.26709594179703283,-0.0035772589306014243,0.007996247088102279} *)


NonlinearModelFit[data, 
  a + b xx + c yy, {a, b, c}, {xx, yy}]["BestFitParameters"]
(* {a -> 0.26709594179703233,b -> -0.0035772589306014156, c -> 0.00799624708810228} *)

(And using lower case letters for your variables is almost always better in Mathematica.)

When you want to compare several models, I find that NonlinearModelFit even when fitting linear models is more convenient. Here are a few possibilities:

nlm1 = NonlinearModelFit[data, a + b xx + c yy, {a, b, c}, {xx, yy}]
nlm2 = NonlinearModelFit[data, 
  a + b xx + c yy + d xx^2 + e yy^2, {a, b, c, d, e}, {xx, yy}]
nlm3 = NonlinearModelFit[data, 
  a + b xx + c yy + d xx yy, {a, b, c, d}, {xx, yy}]
nlm4 = NonlinearModelFit[data, 
  a + b xx + c yy + d xx^2 + e yy^2 + f xx yy, {a, b, c, d, e, f}, {xx, yy}]

You can make relative comparisons the models using AIC:

nlm1["AIC"]
(* -315.55322687806864 *)
nlm2["AIC"]
(* -342.34325852838634 *)
nlm3["AIC"]
(* -485.8345586581582 *)
nlm4["AIC"]
(* -569.8967227522896 *)

Here we see that model 4 with the smallest AIC gives the best fit of the 4 models. That doesn't mean you get a good fit. To get an idea of a more absolute fit one can look at the root mean square error:

nlm1["EstimatedVariance"]^0.5
(* 0.009818046569044022 *)
nlm2["EstimatedVariance"]^0.5
(* 0.007374767809867794 *)
nlm3["EstimatedVariance"]^0.5
(* 0.0017720941160463222 *)
nlm4["EstimatedVariance"]^0.5
(* 0.0007510869942900922 *)

Here model 4 performs best in that it has the smallest root mean square error but only you know if that value is small enough to meet your objective(s). For more details of model comparisons, Cross Validated is a more appropriate forum.

Update

Once you determine the approach for model fitting, you can use Mathematica to check on the model fit by looking at the residuals as one is assuming independent and normally distributed errors with a constant variance. Below is a comparison of Model 1 and Model 4 showing the predicted vs residual, histograms of residuals, and a quantile-quantile plot. Note that Model 4 has better behaved and smaller residuals than Model 1. (Why you want to use these visuals is better addressed in Cross Validated.)

(* Calclate residuals *)
residuals1 = z - nlm1["PredictedResponse"];
residuals4 = z - nlm4["PredictedResponse"];

(* Plot predicted response vs. residual *)
GraphicsRow[{ListPlot[
   Transpose[{nlm1["PredictedResponse"], residuals1}],
   PlotRange -> {{0, 1.4}, Full},
   PlotLabel -> Style["Model 1", Bold, Large], Frame -> True,
   FrameLabel -> {"Predicted", "Residual"}],
  ListPlot[Transpose[{nlm4["PredictedResponse"], residuals4}],
   PlotRange -> {{0, 1.4}, Full},
   PlotLabel -> Style["Model 4", Bold, Large], Frame -> True,
   FrameLabel -> {"Predicted", "Residual"}]},
 ImageSize -> Full]

(* Histograms of residuals *)
GraphicsRow[{
  Show[Histogram[residuals1, Automatic, "PDF", PlotRange -> Full, Frame -> True,
    FrameLabel -> {"Residual", "Density"}],
   Plot[PDF[NormalDistribution[Mean[residuals1], StandardDeviation[residuals1]], x], {x, -0.05, 0.02}, 
    PlotStyle -> Thick]],
  Show[Histogram[residuals4, Automatic, "PDF", PlotRange -> Full, Frame -> True,
    FrameLabel -> {"Residual", "Density"}],
   Plot[PDF[NormalDistribution[Mean[residuals4],  StandardDeviation[residuals4]], x], {x, -0.0025, 0.0025}, PlotStyle -> Thick]]},
 ImageSize -> Full]

(* Quantile-quantile plot *)
GraphicsRow[{QuantilePlot[residuals1, 
   NormalDistribution[Mean[residuals1], 
    StandardDeviation[residuals1]], 
   FrameLabel -> {"Normal theoretical quantiles", "Data quantiles"}],
  QuantilePlot[residuals4, 
   NormalDistribution[Mean[residuals4], StandardDeviation[residuals4]],
   FrameLabel -> {"Normal theoretical quantiles", "Data quantiles"}]},
  ImageSize -> Full]

Residual plots

Improve this answer
$\endgroup$
1
  • $\begingroup$ Very good and insightful answer. (+1) $\endgroup$
    – Anton Antonov
    Commented Jan 16, 2017 at 14:44
2
$\begingroup$

It helps to look at the data..your points fall very nearly on a plane, which can be found readily:

 plane = Last@Minimize[Total[(a X + b Y - Z + c)^2], {a, b, c}]

{a -> -0.00357726, b -> 0.00799625, c -> 0.267096}

poly = (Append[#, a #[[1]] + b #[[2]] + c /. plane] & /@ {{50, 
     19}, {76, 19}, {76, 230}, {50, 230}})
g1=Graphics3D[{ {Red, PointSize[.01], Point[Transpose[{X, Y, Z}]] }, 
  Opacity[.5], Polygon[poly] }, BoxRatios -> 1]

enter image description here

now it looks like the majority of your points are on a near perfect line, with some scatter.. so do a nonlinear fit to a line:

u = {ux, uy, uz};
sol = FindMinimum[
  Total[Norm[ # -  (# - {x0, y0, z0}).u u/Total[u^2]  - {x0, y0, 
        z0} ] & /@ Transpose[{X, Y, Z}]], {{x0, 75}, {y0, 
    50}, {z0, .37}, {ux, -1}, {uy, 0}, {uz, 0}}]

Show[g1, Graphics3D[
  Line[{{x0, y0, z0}, {x0, y0, z0} + 30 u} /. sol[[2]]], 
  BoxRatios -> 1]]

enter image description here

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ I agree plotting the data is essential and we don't know how the data was collected or why many of the sample points are on a line as you show. But what is the purpose of fitting that line when the objective is to predict Z from X and Y? In other words, just because that artifact appears in the data, why attempt to fit that line? What am I missing? It seems that that line of points is readily apparent without drawing a line through it. $\endgroup$
    – JimB
    Commented Dec 30, 2015 at 23:18
  • $\begingroup$ Assuming all those scattered points are significant, the plane is the answer to your question (Z=a X + b Y + c). I just thought it was interesting to show how to fit the line in 3d space. $\endgroup$
    – george2079
    Commented Dec 31, 2015 at 15:30
  • $\begingroup$ Thanks. That makes sense. (Sorry, I was just drawing a blank.) The x,y values are certainly not chosen to allow for a more uniform coverage. It looks like the "measuring device" was set to a particular x value and the the y was varying roughly uniformly and vice versa plus a few other points. And predictions should not taken too seriously in the region where there are no x and y values. Adding in confidence regions might highlight that. $\endgroup$
    – JimB
    Commented Dec 31, 2015 at 16:15

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.