# How to fit polynomial curve using Mathematica for multiple inputs?

I am new to mathematica. I have 3 inputs and 1 output. I want to find the 'N' th degree of polynomial which would approimately fit my dataset. I tried FindFit but it does not solve my problem. I also tried Neural Networks but want to do Polynomial Curve Fitting. I want to predict what degree of polynomial would fit to my data as the relation is Non-Linear(i.e. F(A)+F(B)+F(C) =Output)

Basically I want to find the 'n' degree polynomial which would approximate my data.

{{"A", "B", "C", "Output"}, {57.0088, 76.2105, 46.4823,
46.8503}, {56.2162, 75.5021, 44.8855, 47.1508}, {57.189, 74.434,
44.8579, 48.4662}, {57.551, 75.3068, 46.6976, 45.9044}, {55.9446,
75.4585, 47.2454, 46.3321}, {58.7154, 76.9871, 47.2429,
50.8479}, {59.7091, 79.1212, 48.5647, 50.7675}, {59.5948, 78.0647,
47.5596, 48.7605}, {59.9859, 77.9062, 50.9319, 48.0986}, {60.5902,
78.1553, 54.5241, 51.1845}, {59.6509, 77.6456, 53.0787,
52.9824}, {62.7154, 74.9587, 53.5889, 52.242}, {62.8849, 81.0653,
55.2839, 50.2263}, {63.2489, 78.4848, 54.355, 51.0166}, {64.6259,
78.9331, 54.9457, 53.4716}, {65.4176, 79.1202, 55.1678,
53.1001}, {64.962, 77.6563, 54.3707, 53.6528}, {63.8401, 79.6976,
51.0827, 54.067}, {64.4447, 82.2205, 52.9363, 55.5474}, {64.7958,
82.149, 55.1345, 57.0765}}


• Acually they are not a mesh, it is weight in (Kg), where the output is a physically measured one. and A,B,C are other weight functions affecting the output. So how to fit polynomial curve to the data? – axay Jul 4 '19 at 10:36
• Please ad links between this and the cross-post on Wolfram Community. – Daniel Lichtblau Jul 4 '19 at 15:35
• A quibble about terminology: A linear model is "linear" in the parameters rather than "linear" in the predictors. So a polynomial is a linear model and polynomials can be fit with LinearModelFit. – JimB Jul 4 '19 at 23:48

Modifying the approach from this answer:

vars = {x, y, z};
maxdegree = 3;
cols = Join @@ (MonomialList[(Plus @@ vars)^#] /. _Integer x_ :> x & /@ Range[0, maxdegree])


{1, x, y, z, x^2, x y, x z, y^2, y z, z^2, x^3, x^2 y, x^2 z, x y^2, x y z, x z^2, y^3, y^2 z, y z^2, z^3}

Specify the maximum number of terms (nparams) including the intercept (1) and construct all possible models with up to nparams terms

nparams = 6;
models = Subsets[cols, {1, nparams}];
Length @ models


60459

Use LinearModelFit for each model with the option IncludeConstantBasis -> False:

fits = Table[Join[{j}, {Length @ j},
LinearModelFit[data, j, vars, IncludeConstantBasis -> False][
{j, models}];


Take the top 20 (say) by "AICc" and display:

topTwentyByAICc = SortBy[fits, #[[4]] &][[;; 20]];
Style[# /. x_Real :> Round[x, .00001]] &@
Grid[{{"Model", "Length", "BestFit", "AICc", "BIC",
"AdjustedRSquared", "RSquared"}, ## & @@ topTwentyByAICc},
Dividers -> All]


Note: As JimB noted in comments, when the model does not include a constant bases (1 is not the basis) R^2 and AdjustedR^2 are meaningless.

• In the SortBy command you don't want the minus sign if front of #[[4]]. (AICc has the log likelihood multiplied already by -1 so low values indicate better models.) You've produced the bottom 20 models rather than the top 20 models. – JimB Jul 5 '19 at 0:26
• And while I'm being picky, $R^2$ and adjusted $R^2$ are nonsensical (or at best misleading) when the intercept is forced to be zero (IncludeConstantBasis -> False). – JimB Jul 5 '19 at 0:43
• @JimB, great points , thank you! I fixed the table and added a note re R^2 when the model does not contain an intercept term. – kglr Jul 5 '19 at 1:36

Polynomial Fit (data as defined in @Thies Heidecke 's answer)

fit[x_, y_, z_] :=Fit[data,
{1, x, y, z, x^2, x y, y^2, x z, y z, z^2, x^3, x^2 y,x y^2, y^3, x^2 z, x y z, y^2 z, x z^2, y z^2, z^3} (* polynomial basis *)
, {x, y, z}] //Evaluate


fits the data very well

Map[#[[-1]] - Apply[fit, Most[#]] &, data]
(* O[10^-10] *)

• Or take the Norm of your result: Norm@((#[[-1]] - fit @@ (Most@#)) & /@ data) – Bob Hanlon Jul 4 '19 at 13:02

This is not a final analysis but just to get you started, let's look at the data

data = {{57.0088, 76.2105, 46.4823, 46.8503}, {56.2162, 75.5021,
44.8855, 47.1508}, {57.189, 74.434, 44.8579, 48.4662}, {57.551,
75.3068, 46.6976, 45.9044}, {55.9446, 75.4585, 47.2454,
46.3321}, {58.7154, 76.9871, 47.2429, 50.8479}, {59.7091, 79.1212,
48.5647, 50.7675}, {59.5948, 78.0647, 47.5596, 48.7605}, {59.9859,
77.9062, 50.9319, 48.0986}, {60.5902, 78.1553, 54.5241,
51.1845}, {59.6509, 77.6456, 53.0787, 52.9824}, {62.7154, 74.9587,
53.5889, 52.242}, {62.8849, 81.0653, 55.2839, 50.2263}, {63.2489,
78.4848, 54.355, 51.0166}, {64.6259, 78.9331, 54.9457,
53.4716}, {65.4176, 79.1202, 55.1678, 53.1001}, {64.962, 77.6563,
54.3707, 53.6528}, {63.8401, 79.6976, 51.0827, 54.067}, {64.4447,
82.2205, 52.9363, 55.5474}, {64.7958, 82.149, 55.1345, 57.0765}}


via ListDensityPlot3D

p1 = ListDensityPlot3D[data, AxesLabel -> {"A", "B", "C"}]


Ok, this looks like it might be approximated sufficiently by a linear gradient. We can get a fit via LinearModelFit or NonlinearModelFit (also Fit and FindFit but they are not as versatile). Let's use NonlinearModelFit in case we want to extend the model to something more complicated later:

fit = NonlinearModelFit[
data,
w[1] a + w[2] b + w[3] c + w[4],
{w[1], w[2], w[3], w[4]},
{a, b, c}
]


FittedModel[-14.9611+0.730883 a+0.294494 b-0.0325676 c]

and plot it in the same way

p2 = DensityPlot3D[
fit[a, b, c],
{a, 55, 66}, {b, 74, 83}, {c, 44, 56},
PlotRange -> AbsoluteOptions[p1, PlotRange][[1, 2]],
AxesLabel -> {"A", "B", "C"}
]


We can use the ConvexHullMesh of the data points to create a RegionMember function to plot our fit in the same region as we got earlier with ListDensityPlot3D:

rm = RegionMember[ConvexHullMesh[data[[All, ;; 3]]]]
rmn[x_?NumericQ, y_?NumericQ, z_?NumericQ, f_] := rm[{x, y, z}]
p3 = DensityPlot3D[
fit[a, b, c], {a, 55, 66}, {b, 74, 83}, {c, 44, 56},
PlotRange -> AbsoluteOptions[p1, PlotRange][[1, 2]],
AxesLabel -> {"A", "B", "C"},
RegionFunction -> rmn
]


From here we could go into different directions (more complicated model, verifying the quality of the fit, etc.), but this might be a start!

Following this question, we can find fit with minimal Akaike Information Criterion ( AIC ) by calculating a series of fits with 1, 2, ...maxn fits, and selecting the fit with the smallest "AIC" as defined in the NonlinearModelFit documentation.

kvar[n_]:={x^n,y^n,z^n};
kvar[0]=1;
kpar[n_]:={a[n],b[n],c[n]};
kpar[0]=d[0];
gmodel[n_Integer]:=kpar[0]+Sum[kvar[i].kpar[i],{i,1,n}];
gpars[n_Integer]:=Flatten@Array[kpar,n+1,{0,n}]
fitg[data_,maxn_Integer]:=MinimalBy[Table[{#,#["AIC"]}&@NonlinearModelFit[data,gmodel[n],gpars[n],{x,y,z}],{n,maxn}],Last][[1,1]]


The data

dat={{57.0088,76.2105,46.4823,46.8503},{56.2162,75.5021,44.8855,47.1508},{57.189,74.434,44.8579,48.4662},{57.551,75.3068,46.6976,45.9044},{55.9446,75.4585,47.2454,46.3321},{58.7154,76.9871,47.2429,50.8479},{59.7091,79.1212,48.5647,50.7675},{59.5948,78.0647,47.5596,48.7605},{59.9859,77.9062,50.9319,48.0986},{60.5902,78.1553,54.5241,51.1845},{59.6509,77.6456,53.0787,52.9824},{62.7154,74.9587,53.5889,52.242},{62.8849,81.0653,55.2839,50.2263},{63.2489,78.4848,54.355,51.0166},{64.6259,78.9331,54.9457,53.4716},{65.4176,79.1202,55.1678,53.1001},{64.962,77.6563,54.3707,53.6528},{63.8401,79.6976,51.0827,54.067},{64.4447,82.2205,52.9363,55.5474},{64.7958,82.149,55.1345,57.0765}};


Noy we try fits up to order 10.

Quiet@fitg[dat,10]
(* FittedModel[-14.9611+0.730883 x+0.294494 y-0.0325676 z] *)


And we find the best fit is linear.