0
$\begingroup$

I am trying to fit the following data to an expression of the form shown in the code below but it doesn't seem to fit that well and I don't know why it isn't as this kind of data is routinely fit to such expressions. I would greatly appreciate any help or pointer on what I am doing wrong. Thanks.

data = {{3., 92.1445}, {3.1, 5.8997}, {3.2, -44.8126}, {3.3, -72.4833}, {3.4, -85.5145}, {3.5, -89.4885}, {3.6, -88.0608}, {3.7, -83.587}, {3.8, -77.5604}, {3.9, -70.9073}, {4., -64.1871}, {4.1, -57.7239}, {4.2, -51.6913}, {4.3, -46.1686}, {4.4, -41.1779}, {4.5, -36.7066}, {4.6, -32.7236}, {4.7, -29.1886}, {4.8, -26.059}, {4.9, -23.2921}, {5., -20.8469}, {5.1, -18.6868}, {5.2, -16.7775}, {5.3, -15.0885}, {5.4, -13.5933}, {5.5, -12.2677}, {5.6, -11.0904}, {5.7, -10.0433}, {5.8, -9.11012}, {5.9, -8.27735}, {6., -7.53272}, {6.1, -6.86538}, {6.2, -6.26625}, {6.3, -5.72797}, {6.4, -5.24286}, {6.5, -4.80566}, {6.6, -4.41044}, {6.7, -4.05333}, {6.8, -3.72946}, {6.9, -3.43566}, {7., -3.1688}, {7.1, -2.92572}, {7.2, -2.70467}, {7.3, -2.50286}, {7.4, -2.31854}, {7.5, -2.14996}, {7.6, -1.99571}, {7.7, -1.85441}, {7.8, -1.72465}, {7.9, -1.60573}, {8., -1.49591}, {8.1, -1.39518}, {8.2, -1.30215}, {8.3, -1.21645}, {8.4, -1.13741}, {8.5, -1.06396}, {8.6, -0.996458}, {8.7, -0.933501}, {8.8, -0.875442}, {8.9, -0.821579}, {9., -0.771564}, {9.1, -0.724697}, {9.2, -0.681676}, {9.3, -0.641105}, {9.4, -0.603681}, {9.5, -0.568705}, {9.6, -0.535828}, {9.7, -0.505399}, {9.8, -0.476719}, {9.9, -0.450137}, {10., -0.425305}, {10.1, -0.401871}, {10.2, -0.380186}, {10.3, -0.35955}, {10.4, -0.340314}, {10.5, -0.322126}, {10.6, -0.305338}, {10.7, -0.289249}, {10.8, -0.274209}, {10.9, -0.260219}, {11., -0.246928}, {11.1, -0.234337}, {11.2, -0.222445}, {11.3, -0.211603}, {11.4, -0.20111}, {11.5, -0.190967}, {11.6, -0.181874}, {11.7, -0.17278}, {11.8, -0.164386}, {11.9, -0.156691}, {12., -0.149346}, {12.1, -0.142001}, {12.2, -0.135356}, {12.3, -0.12906}, {12.4, -0.123114}, {12.5, -0.117518}, {12.6, -0.111922}, {12.7, -0.106676}, {12.8, -0.102129}, {12.9, -0.0972324}, {13., -0.0930354}, {13.1, -0.0888383}, {13.2, -0.0846412}, {13.3, -0.0811436}, {13.4, -0.0772963}, {13.5, -0.0741485}, {13.6, -0.0706509}, {13.7, -0.0675031}, {13.8, -0.064705}, {13.9, -0.061907}, {14., -0.0591089}, {14.1, -0.0566606}, {14.2, -0.0542123}, {14.3, -0.051764}, {14.4, -0.0496655}, {14.5, -0.047567}, {14.6, -0.0454684}, {14.7, -0.0433699}, {14.8, -0.0416211}, {14.9, -0.0398723}, {15., -0.0381235}, {15.1, -0.0363747}, {15.2, -0.0349757}, {15.3, -0.0335767}, {15.4, -0.0321776}, {15.5, -0.0307786}, {15.6, -0.0293796}, {15.7, -0.0279806}, {15.8, -0.0269313}, {15.9, -0.025882}, {16., -0.024483}, {16.1, -0.0234337}, {16.2, -0.0223844}, {16.3, -0.0216849}, {16.4, -0.0206357}, {16.5, -0.0195864}, {16.6, -0.0188869}, {16.7, -0.0181874}, {16.8, -0.0171381}, {16.9, -0.0164386}, {17., -0.0157391}, {17.1, -0.0150396}, {17.2, -0.01434}, {17.3, -0.0136405}, {17.4, -0.012941}, {17.5, -0.0125913}, {17.6, -0.0118917}, {17.7, -0.0111922}, {17.8, -0.0108425}, {17.9, -0.010143}, {18., -0.0097932}, {18.1, -0.00909368}, {18.2, -0.00874393}, {18.3, -0.00839417}, {18.4, -0.00769465}, {18.5, -0.0073449}, {18.6, -0.00699514}, {18.7, -0.00664538}, {18.8, -0.00629563}, {18.9, -0.00594587}, {19., -0.00559611}, {19.1, -0.00524636}, {19.2, -0.0048966}, {19.3, -0.00454684}, {19.4, -0.00419708}, {19.5, -0.00384733}, {19.6, -0.00349757}, {19.7, -0.00349757}, {19.8, -0.00314781}, {19.9, -0.00279806}, {20., -0.0024483}, {20.1, -0.0024483}, {20.2, -0.00209854}, {20.3, -0.00174879}, {20.4, -0.00174879}, {20.5, -0.00139903}, {20.6, -0.00139903}, {20.7, -0.00104927}, {20.8, -0.000699514}, {20.9, -0.000699514}, {21., -0.000349757}, {21.1, -0.000349757}}

The model expression is:

vdc[RR_] := 
  Module[{R = RR, \[Alpha]0 = 16.36606, \[Alpha]1 = 
     0.70172, \[Beta]0 = 17.19338, \[Beta]1 = 0.09574},
   A[num1_Integer] := \[Alpha]0 num1^-\[Alpha]1;
   B[num2_Integer] := \[Beta]0 Exp[-\[Beta]1*num2];
   \[Rho] = 5.5 + 1.25 R0;
   R0 = 6.46;
   \[Chi][R_, 
     n_Integer] := (1 - Exp[-(A[n]*R/\[Rho]) - (B[n]*R^2/\[Rho]^2)])^n;
   v[R_] = Sum[\[Chi][R, n]*Subscript[CC, n]/R^n, {n, 6, 10, 2}]
   ];
vtotal[R_] := (a* Exp[-b * R]) - vdc[R]

and here is my trial code

dataplot = ListPlot[data, PlotRange -> {{2., 8.0}, {-400.0, 700}}, ImageSize -> Large];
fit = FindFit[data, 
  vtotal[R], {Subscript[CC, 6], Subscript[CC, 8], Subscript[CC, 10], 
   a, b}, R, Method -> NMinimize, AccuracyGoal -> Infinity];

Show[{dataplot, Plot[vtotal[R] /. fit, {R, 2., 8.0}]}, 
 AxesLabel -> {"R", 
   "Potential Energy,\!\(\*SuperscriptBox[\(cm\), \(-1\)]\)"}]
```
$\endgroup$
6
  • $\begingroup$ Have a look at NonlinearModelFit. $\endgroup$ Commented Oct 14, 2020 at 13:26
  • $\begingroup$ How important is to use "a model expression" for your fit? If it is not, you can try FindFormula. $\endgroup$ Commented Oct 14, 2020 at 13:44
  • $\begingroup$ Anton Antonov the model if very important in this problem. $\endgroup$
    – fred85
    Commented Oct 14, 2020 at 17:45
  • 1
    $\begingroup$ Your data are not asymptotically consistent with you model: the model goes as $1/R^6$ for large $R$ (verify with Asymptotic[vtotal[R], R -> ∞]) but your data do not: try ListPlot[{#[[1]], #[[1]]^6*#[[2]]} & /@ data] (does not approach an asymptote for large $R$). Maybe your long-range dispersion forces are not calculated correctly. $\endgroup$
    – Roman
    Commented Oct 14, 2020 at 18:00
  • $\begingroup$ Hi Roman, the model is a double many body expansion (DMBE) and converges a 1/R$^6$ as you correctly predict. This data is a slice of a 3D potential fit to an Ab Initio data calculated at CCSD(T)-F12b/CBS. I think the data is accurate enough and I have edited it to exclude the highly repulsive region and one can as well exclude points beyond R=21 $\endgroup$
    – fred85
    Commented Oct 14, 2020 at 18:50

2 Answers 2

3
$\begingroup$

If in your original dataset you had retained the large values associated with $R<3$, then @DanielHuber 's approach is a good one (which is one way to account for varying precision of the observations).

But with your modified dataset with $R\geq 3$, there is less need for weighting. However, because the data and model structure induces high correlations among the parameter estimators and because of the wildly different scales of the parameter values, having good starting values for the parameters is essential.

nlm = NonlinearModelFit[data, vtotal[R], 
  {{a, 4.75 10^8}, {b, 4.8}, {Subscript[CC, 6], 900000}, {Subscript[CC, 8], -10^8},
   {Subscript[CC, 10], 2.6 10^9}}, R];
nlm["CorrelationMatrix"] // MatrixForm

Correlation matrix

nlm["ParameterTable"]

Parameter table

Show[ListPlot[data, PlotRange -> {{3, 8}, {-100, 100}}], 
 Plot[nlm[R], {R, 3, 8}]]

Data and fit

One can see the high correlations large standard errors for most of the parameter estimators.

Using NonlinearModelFit allows one to obtain diagnostic statistics which allows checks on the appropriateness of the fit. You can't get that with FindFit. (Although FindFit now has the NormFunction which is great for robust regression models.)

And to confess: I didn't use magic or hallucinogenic drugs to come up with the starting values. The model being fit is essentially a linear model with the exception of the Exp[-b R] term. So I just tried a bunch of values for b, used LinearModelFit on the remaining parameters, and looked at the best fit. Then I plugged those values into NonlinearModelFit.

$\endgroup$
5
$\begingroup$

Look at your data:

enter image description here

If you try a least square fit to this data, MMA mainly tries to fit the leftmost end because that is where the large errors are. Therefore you must in some way tell MMA, that the errors at the left end are less important. That is you need a weighted fit. In FindFit you may specify a norm function, that the actual error to be minimized calculates. So, we give some function (you may try different functions), that weigh the data at at the left end less:

fit = FindFit[data, 
   vtotal[R], {Subscript[CC, 6], Subscript[CC, 8], Subscript[CC, 10], 
    a, b}, R, NormFunction -> (Norm[#^-0.4 #] &)];
Plot[vtotal[R] /. fit, {R, 2., 8.0}, Epilog -> Point[data], 
 PlotRange -> All]

enter image description here

$\endgroup$
10
  • $\begingroup$ You hit the nail on the head about the large response values determining the solution with least squares. But even with MaxIterations -> 5000 I don't get convergence with your code. Also, what does #^-0.4 # mean? Does it mean #^0.6 ? $\endgroup$
    – JimB
    Commented Oct 14, 2020 at 16:42
  • $\begingroup$ The region of interest is in the range R >=3 thus I have edited the data to remove the high energy data points (where R < 3) $\endgroup$
    – fred85
    Commented Oct 14, 2020 at 17:55
  • $\begingroup$ @JimB I wrote Norm[ #^-0.4 # ] for didactical reasons to separate the weight from the error. $\endgroup$ Commented Oct 14, 2020 at 18:14
  • $\begingroup$ @DanielHuber. Sorry that 4-syllable word is too much for me. Is there something in the documentation that tells me what Norm[ #^-0.4 # ] means if it doesn't mean Norm[#^0.6] ? $\endgroup$
    – JimB
    Commented Oct 14, 2020 at 20:20
  • $\begingroup$ @JimB do not worry, it is the same. $\endgroup$ Commented Oct 15, 2020 at 8:21

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.