1
$\begingroup$

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}}

Poly

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Axay Jul 4 at 10:36
  • $\begingroup$ Please ad links between this and the cross-post on Wolfram Community. $\endgroup$ – Daniel Lichtblau Jul 4 at 15:35
  • 2
    $\begingroup$ 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. $\endgroup$ – JimB Jul 4 at 23:48
4
$\begingroup$

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][
      {"BestFit", "AICc", "BIC", "AdjustedRSquared", "RSquared"}]], 
   {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]

enter image description here

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – JimB Jul 5 at 0:26
  • $\begingroup$ 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). $\endgroup$ – JimB Jul 5 at 0:43
  • $\begingroup$ @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. $\endgroup$ – kglr Jul 5 at 1:36
1
$\begingroup$

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] *) 
$\endgroup$
  • $\begingroup$ Or take the Norm of your result: Norm@((#[[-1]] - fit @@ (Most@#)) & /@ data) $\endgroup$ – Bob Hanlon Jul 4 at 13:02
1
$\begingroup$

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"}]

ListDensityPlot3D

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"}
     ]

DensityPlot3D of fit

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
     ]

DensityPlot3D of fit confined to same region as original data domain

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

$\endgroup$
1
$\begingroup$

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.

$\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.