# Dynamic number of parameter fits in Mathematica

I am trying to learn if in Mathematica one can perform fits of analytic functions that are expanded and written as power series. In particular, we can consider the following generic form:

$$f(B,x) = a_0 + \sum_{i=1}^{n}a_i [B^{1/\mu}(x-x_c)]^i + m_0B^{-k} \tag{1}$$

where $$\mu, k$$ and $$x_c$$ are known to be positive. To fit such a function to a given data set, one can of course assume up to which order the sum should to be taken, then redefine the fitting function by writing down the terms explicitly up to that order and re-doing the fit. But this would be quite inefficient if we had to go to orders of e.g. $$n=4$$ or $$n=8.$$

So I was wondering, if one could for a given order $$n$$, perform the fit w.r.t to its corresponding fit parameters without having to rewrite the terms of the sum explicitly by hand. For instance for $$n=3,$$ we would have $$8$$ fitting parameters ($$a_0,a_1,a_2,a_3,\mu, x_c, m_0, k$$), and for other $$n$$'s we would have a different number of these $$a$$ coefficients.

• Is there a way to define $$f$$ given in $$(1)$$ in Mathematica as a fitting function without explicitly writing the terms up to chosen order $$n?$$ Any hints or links to similarly solved problems would be very helpful, the goal is just to learn how such fittings could be performed, therefore, the provided example and dataset are just dummy examples to discuss with.

As sample data, I have created the following set ($$x$$ values, yvals are corresponding $$f$$ values and $$B=200,$$ for other $$B$$'s, the rest of the dataset can be taken from here):

xvals = {0.300000, 0.312245, 0.324490, 0.336735, 0.348980, 0.361224,
0.373469, 0.385714, 0.397959, 0.410204, 0.422449, 0.434694,
0.446939, 0.459184, 0.471429, 0.483673, 0.495918, 0.508163,
0.520408, 0.532653, 0.544898, 0.557143, 0.569388, 0.581633,
0.593878};
yvals = {0.005568, 0.005930, 0.006655, 0.008103, 0.010992, 0.016728,
0.028021, 0.049881, 0.090847, 0.163149, 0.278213, 0.434257,
0.605564, 0.755021, 0.861433, 0.926721, 0.963188, 0.982494,
0.992427, 0.997461, 0.999994, 1.001263, 1.001898, 1.002215,
1.002374};

• A plot of your data show an "S" shaped curve. Polynomials don't tend to provide good fits to such shapes. In any event you can include a variable number of parameters in NonlinearModelFit using something like this: Join[{{\[Mu], 1.04}, {xc, 0.2}, {m0, 44}, {k, 44}}, Table[{a[i], 0}, {i, 0, n}]] where you've defined n as the order of the polynomial. – JimB Jun 14 '19 at 18:39
• @JimB hi, thx a lot for the helpful comment. Indeed ultimately fitting the whole dataset is not the intention, the form of $f$ as shown in the Eq is supposed to be valid near the point $x_c$ (mid-way through the curve) i think your approach might do the trick and i ve already learnt something new. Would you be so kind to write your approach in a more expanded form as an answer so I can accept and also it might be useful for future readers. Thanks again! – user929304 Jun 14 '19 at 21:09
• The way you have written the equations you really only have n+2 fitting parameters. Using a new variable c you can always write ai * B^1/\[Mu] as ci since B and \[Mu] are constants. Similarly a0 + m0*B^-k can be written as c0. The parameters B, \[Mu] and k are not needed. – Jack LaVigne Jun 14 '19 at 22:32
• @JackLaVigne very good point. But in this particular case, for instance $\mu$ itself is a critical exponent to be estimated from the fit, so could we even in such a case go through with the redefinition of constant? – user929304 Jun 15 '19 at 1:29
• You've now added in a data set which might make $\mu$ estimable. You should then remove the "Accept" from my answer to encourage others to help as I'm off playing for much of the weekend. – JimB Jun 15 '19 at 16:29

Following @JackLaVigne 's insightful comment, you've recognized that to be able to estimate the parameters of interest multiple values of $$B$$ are required and then provided that data. But that is not the end of your troubles.

Here is one way to estimate the coefficients in the model that you propose:

Construct the overall dataset:

xvals100 = {0.300000, 0.312245, 0.324490, 0.336735, 0.348980, 0.361224, 0.373469, 0.385714,
0.397959, 0.410204, 0.422449, 0.434694, 0.446939, 0.459184, 0.471429, 0.483673, 0.495918,
0.508163, 0.520408, 0.532653, 0.544898, 0.557143, 0.569388, 0.581633, 0.593878};
yvals100 = {0.026580, 0.032527, 0.040959, 0.052835, 0.069416, 0.092282, 0.123279, 0.164343,
0.217116, 0.282346, 0.359204, 0.444817, 0.534413, 0.622262, 0.703059, 0.773115, 0.830819,
0.876372, 0.911141, 0.937004, 0.955875, 0.969452, 0.979121, 0.985957, 0.990765};
xvals150 = {0.300000, 0.312245, 0.324490, 0.336735, 0.348980, 0.361224, 0.373469, 0.385714,
.397959, 0.410204, 0.422449, 0.434694, 0.446939, 0.459184, 0.471429, 0.483673, 0.495918,
0.508163, 0.520408, 0.532653, 0.544898, 0.557143, 0.569388, 0.581633, 0.593878};
yvals150 = {0.009346, 0.012428, 0.017519, 0.025888, 0.039520, 0.061408, 0.095748, 0.147719,
0.222233, 0.321193, 0.440052, 0.566725, 0.685633, 0.784671, 0.859267, 0.911305, 0.945695,
.967617, 0.981271, 0.989653, 0.994753, 0.997840, 0.999701, 1.000822, 1.001496};
xvals200 = {0.300000, 0.312245, 0.324490, 0.336735, 0.348980, 0.361224, 0.373469, 0.385714,
0.397959, 0.410204, 0.422449, 0.434694, 0.446939, 0.459184, 0.471429, 0.483673, 0.495918,
0.508163, 0.520408, 0.532653, 0.544898, 0.557143, 0.569388, 0.581633, 0.593878};
yvals200 = {0.005568, 0.005930, 0.006655, 0.008103, 0.010992, 0.016728, 0.028021, 0.049881,
0.090847, 0.163149, 0.278213, 0.434257, 0.605564, 0.755021, 0.861433, 0.926721, 0.963188,
0.982494, 0.992427, 0.997461, 0.999994, 1.001263, 1.001898, 1.002215, 1.002374};
data100 = Transpose[{xvals100, ConstantArray[100, Length[xvals100]], yvals100}];
data150 = Transpose[{xvals150, ConstantArray[150, Length[xvals150]], yvals150}];
data200 = Transpose[{xvals200, ConstantArray[200, Length[xvals200]], yvals200}];
data = Join[data100, data150, data200];


Define regression function and estimate coefficients:

f[x_, a_, b_, μ_, xc_, m0_, k_] :=
Sum[a[[i + 1]] (b^(1/μ) (x - xc))^i, {i, 0, Length[a] - 1}] + m0 b^(-k)

n = 6;
nlm = NonlinearModelFit[
data, {f[x, Table[a[i], {i, 0, n}], b, μ, xc, m0, k], xc > 0 && μ > 0},
Join[{{μ, 1}, {xc, Min[data[[All, 1]]]}, {m0, 1}, {k, 1}},
Table[{a[i], 1}, {i, 0, n}]], {x, b}];

Show[ListPlot[{data100[[All, {1, 3}]], data150[[All, {1, 3}]],
data200[[All, {1, 3}]]},
PlotLegends -> {"B = 100", "B = 150", "B = 200"}],
Plot[{nlm[x, 100], nlm[x, 150], nlm[x, 200]},
{x, Min[data[[All, 1]]], Max[data[[All, 1]]]}]]


You can see that the fits when fitting all data with the same $$a_i$$ coefficients are not very good and yet the individual fits with a single value of $$B$$ are not too horrible. This suggests that a common set of parameters for different values of $$B$$ is not supported by the data.

If your curves are S-shaped with asymptotes, you probably want to try something other than polynomials. And no matter what kind of curve you have you should be careful (or even avoid) high order polynomials because of numerical precision issues that could arise.

The approach that @OkkesDelgurci will provide you a much better fit than your restricted model. The potential downside is that you don't get the parameterization that you desire. But as you can see from the poor fits above the model that generates the data is more complicated than the model you propose. Individual fits for each value of $$B$$ are better but then you can't estimate $$B$$, $$\mu$$, $$m_0$$, and $$k$$. But getting values does not do one good if the fits are poor and apparently inconsistent among datasets. No point in getting a value of $$B$$ under a very wrong model.

My advice is to consult with a statistician and bring not just the data but the rationale for the need to estimate your "critical" parameter.

• The issue is one of estimability and not whether $\mu$ is a critical parameter. If I can elaborate on what @JackLaVigne stated: You can estimate $a_i B^{i/\mu }$ but there are an infinite number of pairs of $(a_i, \mu)$ that result in the same estimate of $a_i B^{i/\mu }$ from which you can't chose the "right one". The data and model you have just can't separate out unique estimates of $\mu$ and $a_i$. To belabor the point if you just knew that $x+y=5$, you can't determine $x$ or $y$ but you can determine zillions of pairs $(x,y)$ that satisfy $x+y=5$. If you also knew $x-y=7$, then... – JimB Jun 15 '19 at 5:09
• Wow neat problem, now i see what you mean! Would we be able to progress towards the estimation of $\mu$ if we had multiple datasets of x and y, each for a different constant $B?$ (e.g. in the sample data above we know $B=200$ and therefore not itself a fit parameter). To better illustrate the kind of fit im trying to perform on the data, you can e.g see the Appendix B of this (arxiv.org/abs/1106.5514) paper. – user929304 Jun 15 '19 at 9:34
• Thank you so much JimB for the elaborate and super clear approach, learned so much already, and for your patience with me. Small question about the defs, it appears b is given twice in the f[] and also in nlm, is this intentional? I couldn't agree more with all your points. Admittedly, I realise now I wasn't clear enough about why this model and how to fit, since I was trying to ask the gist of the thing without giving unncecessary details. In fact, the form of f as given, is a power series which stems from the expansion of the underlying scaling law at the critical point (...) – user929304 Jun 16 '19 at 11:08
• (...) where it is expected to behave analytically as we have finite systems ($B$'s are those chosen sizes). This is similar to (bond/site) percolation transitions, where in inf systems it is a step function, but in finite ones analytic, and therefore we can fit this form of f to it (where f can be the percolation strength and x the occupation probability). So in our case here, we want to limit the fit between y-axis values roughly from 0.2 to 0.8, and not fit the entire curve, since indeed for the entire S shape we have other options as suggested in the other answer, or a tanh (...) – user929304 Jun 16 '19 at 11:09
• (...) but that would miss the critical behaviour. I hope this gives it a bit of background, of interest is to extract the common intersection point which yields the threshold xc and the critical exponent mu. And in literature people often argue, that by finding m0,k and subtracting it from the data, then a common crossing point can be obtained, as finite size errors are removed. I wonder, is there a built-in way to estimate the chi-squared per dof of the fit? – user929304 Jun 16 '19 at 11:10

"All models are wrong, but some are useful" George Box

When I saw your data, I thought Hill model fit pretty good see here.

ClearAll["Global*"]
xvals200 = {0.300000, 0.312245, 0.324490, 0.336735, 0.348980,
0.361224, 0.373469, 0.385714, 0.397959, 0.410204, 0.422449,
0.434694, 0.446939, 0.459184, 0.471429, 0.483673, 0.495918,
0.508163, 0.520408, 0.532653, 0.544898, 0.557143, 0.569388,
0.581633, 0.593878};
yvals200 = {0.005568, 0.005930, 0.006655, 0.008103, 0.010992,
0.016728, 0.028021, 0.049881, 0.090847, 0.163149, 0.278213,
0.434257, 0.605564, 0.755021, 0.861433, 0.926721, 0.963188,
0.982494, 0.992427, 0.997461, 0.999994, 1.001263, 1.001898,
1.002215, 1.002374};

data200 = Transpose[{xvals200, yvals200}];

g[x_] :=  (x^p)/(xs^p + x^p)

obj = Sum[(g@data200[[k, 1]] - data200[[k, 2]])^2, {k,
Length@data200}];

fit = NMinimize[{obj, {p > 0, xs > 0}}, {p, xs},
Method -> "DifferentialEvolution", MaxIterations -> 1000]


{0.00111351, {p -> 24.8378, xs -> 0.438781}}

Thread[{p, xs} = {p, xs} /. Last@fit];

Show[Plot[g[x], {x, 0.3, 0.6}], ListPlot[data200, PlotStyle -> Red]]
`

• Nice spot, i did not know about this model! Thx for performing the fit. In fact the model i ve given for $f$ is only supposed to describe the data near the critical point, and indeed not the entire S curve, which has this shape since in infinite limit it would be step function but finite systems smear out this behaviour. I ve written more details in comments above with JimB. – user929304 Jun 16 '19 at 11:14