# Not happy with my fit results. How to improve it?

I am currently working on a topic involving nuclear theory (the study of deformed nuclei). So for the past two days, I've been trying to fit some data; I developed a theoretical model that aims at describing the energy spectrum of a nucleus.

The analytical expressions for the energies have two free parameters, so I need to find those parameters such that the theoretical energies fit the experiment as much as possible. I am used to Wolfram Mathematica, so I decided to solve this problem with the help of it.

The structure of my problem is like this: I have four energy levels (pairs of x,y values, where x is the spin and y is the energy). The model has an analytical expression for the excitation energy, and it depends on two free parameters and some other values which I know beforehand.
Now, all I need is to find the two free parameters and then represent the results with a series of plots.

I am using the NonLinearModelFit function since the formulas are quite complex (and non-linear). With the obtained parameters, I go and give numerical values for the analytical expressions for the energy, the calculate the root mean square deviation (RMS) by using the well-known formula $$E_\text{RMS}=\sqrt{\frac{\sum_i(E_{exp}^i-E_{th}^i)^2}{N}}$$.

After some effort, I've managed to get around 0.5 MeV. However, in our problem, such a deviation is quite poor. We aim at obtaining an RMS between 0.2-0.3 (0.35 at worst). From what I can see, the issue is band-4, where half of the theoretical values are much smaller than the experimental ones.

I was also thinking about introducing some conditions in the fitting function (e.g. I already know that the second parameter s must be positive) but I don't know which kind of conditions would help me here. You can also see that the analytical expressions for energies have an If statement to them:

that is because when the NonLinearModelFit starts to compute the best-fit parameters, it goes to some values for which the energies are complex. So to avoid that, I just make the energy a big number whenever its Imaginary part is different from zero. In this way, I can avoid the situation where the program runs for "useless" regions of the free parameters.

Q: Is there any way of manipulating that Fit function to decrease that RMS? (no more than 0.3) Or maybe just consider another approach rather than the current one?

Sorry for the long text, but I wanted to give you some context. The attached document contains all my code, with some explanations (a dropbox link for the file here. The analytic expressions for the energies and everything else are irrelevant. The only essential part is the Fit function, which is marked by the green background.

Some details about the document: It has the four excited bands and a joint data set with an index attached to them so that each of the four analytic formulas energyTSD1, energyTSD2, energyTSD3, energyTSD4 are applied on the corresponding data points.

• Note that RMS can be obtained from using the "EstimatedVariance" option of the result of NonlinearModelFit (although the numerator is the sample size minus the number of parameters rather than just the number of observations). The two best ways to have a better fit is to have a better model or better data. Short of that, there's nothing one can do as the NonlinearModelFit does what it is supposed to do with the model and data. Restricting the parameters will only give you a worse fit.
– JimB
Mar 25, 2020 at 5:19

While getting good fits with theoretical models is most desirable/satisfying, if good predictions is the overall objective, why not just fit the data with a simple model?

Also the objective in your linked notebook seems to be the estimation of two parameters ($$a$$ and $$b$$) which are assumed to be common to 4 datasets. That assumption might not likely be supported by the data and getting precise estimates to those parameters for a model that doesn't predict well for any of the datasets can be a recipe for disaster.

The figures in your notebook show poor fits with the theoretical models:

Fitting a simple quadratic reproduces the data much, much better:

(* Obtain separate datasets *)
dataX
(* {{1, 8.5, 0.1966}, {1, 10.5, 0.4597}, {1, 12.5, 0.7746}, {1, 14.5,
1.1609}, {1, 16.5, 1.6112}, {1, 18.5, 2.1265}, {1, 20.5,
2.7051}, {1, 22.5, 3.3441}, {1, 24.5, 4.0411}, {1, 26.5,
4.7937}, {1, 28.5, 5.5992}, {1, 30.5, 6.457}, {1, 32.5, 7.3667}, {1,
34.5, 8.3293}, {1, 36.5, 9.3458}, {1, 38.5, 10.4169}, {1, 40.5,
11.5431}, {1, 42.5, 12.7224}, {1, 44.5, 13.9491}, {1, 46.5,
15.2181}, {1, 48.5, 16.5221}, {2, 13.5, 1.3394}, {2, 15.5,
1.7467}, {2, 17.5, 2.2184}, {2, 19.5, 2.7527}, {2, 21.5,
3.3484}, {2, 23.5, 4.003}, {2, 25.5, 4.7143}, {2, 27.5, 5.4805}, {2,
29.5, 6.3004}, {2, 31.5, 7.1733}, {2, 33.5, 8.0998}, {2, 35.5,
9.08}, {2, 37.5, 10.1147}, {2, 39.5, 11.2036}, {2, 41.5,
12.3466}, {2, 43.5, 13.5441}, {2, 45.5, 14.7911}, {3, 16.5,
2.1237}, {3, 18.5, 2.6293}, {3, 20.5, 3.1973}, {3, 22.5,
3.8243}, {3, 24.5, 4.5094}, {3, 26.5, 5.2506}, {3, 28.5,
6.0465}, {3, 30.5, 6.8963}, {3, 32.5, 7.7988}, {3, 34.5,
8.7546}, {3, 36.5, 9.7638}, {3, 38.5, 10.8268}, {3, 40.5,
11.9392}, {3, 42.5, 13.0861}, {4, 23.5, 4.58}, {4, 25.5,
5.2251}, {4, 27.5, 5.9273}, {4, 29.5, 6.6819}, {4, 31.5,
7.4919}, {4, 33.5, 8.3573}, {4, 35.5, 9.2778}, {4, 37.5,
10.2535}, {4, 39.5, 11.2851}, {4, 41.5, 12.3701}} *)
data1 = Select[dataX, #[[1]] == 1 &][[All, {2, 3}]];
data2 = Select[dataX, #[[1]] == 2 &][[All, {2, 3}]];
data3 = Select[dataX, #[[1]] == 3 &][[All, {2, 3}]];
data4 = Select[dataX, #[[1]] == 4 &][[All, {2, 3}]];


Now fit the separate datasets with a quadratic polynomial:

nlm1 = NonlinearModelFit[data1, a0 + a1 x + a2 x^2, {a0, a1, a2}, x];
nlm2 = NonlinearModelFit[data2, a0 + a1 x + a2 x^2, {a0, a1, a2}, x];
nlm3 = NonlinearModelFit[data3, a0 + a1 x + a2 x^2, {a0, a1, a2}, x];
nlm4 = NonlinearModelFit[data4, a0 + a1 x + a2 x^2, {a0, a1, a2}, x];

(* Estimates of root mean square error *)
#["EstimatedVariance"]^0.5 & /@ {nlm1, nlm2, nlm3, nlm4}
(* {0.03234, 0.0113766, 0.010037, 0.000625011} *)

(* Show results *)
GraphicsGrid[{{Show[ListPlot[data1],
Plot[a0 + a1 x + a2 x^2 /. nlm1["BestFitParameters"], {x, 0, 50}], PlotLabel -> "1"],
Show[ListPlot[data2],
Plot[a0 + a1 x + a2 x^2 /. nlm2["BestFitParameters"], {x, 0, 50}], PlotLabel -> "2"]},
{Show[ListPlot[data3],
Plot[a0 + a1 x + a2 x^2 /. nlm3["BestFitParameters"], {x, 0, 50}], PlotLabel -> "3"],
Show[ListPlot[data4],
Plot[a0 + a1 x + a2 x^2 /. nlm4["BestFitParameters"], {x, 0, 50}], PlotLabel -> "4"]}}]


• Than you for the answer. The issue is that I can't do separate plots since the entire data should be reproduced by only two parameters, which must be the same across all three energy bands (I'm referring to their analytic formulas). As you said, the behavior of the theoretical curves is indeed dictated by a quadratic function, but those parameters are the crucial part in my problem. Mar 26, 2020 at 4:54
• If you fit the 4 models separately, you'll find that 2 and 4 fit very well but with very different estimates for $b$. Again, the lack of fit is not due to the fitting procedure (or any fitting procedure). The models and the data just don't work well together.
– JimB
Mar 26, 2020 at 5:46