7
$\begingroup$

Alright I have some data I want to fit using a function I know.

data={7783.52,7763.52,7804.15,7782.85,7606.1,7529.88,7541.92,7498.92,7544.81,7594.85,7469.69,7222.08,6958.88,6879.17,6889.06,6891.73,6805.94,6700.71,6478.79,6286.9,6025.6,5812.06,5700.92,5617.63,5562.48,5581.15,5482.96,5446.83,5248.6,5118.35,5030.38,4911.35,4717.25,4638.33,4450.75,4244.48,4109.38,4095.5,3989.44,4023.46,3960.42,3937.88,3890.94,3739.92,3556.4,3417.06,3282.96,3178.1,3129.85,2998.65,2943.69,2852.4,2824,2827.29,2862.67,2784.71,2658.15,2571.81,2490.88,2426.5,2405.88,2411.5,2476.65,2460.15,2391.88,2395.1,2335.58,2287.94,2287.44,2258.25,2213.81,2145.44,2083.92,2016.02,2007.75,2086.63,2128.38,2023.54,2008.46,2018.67,2066.35,2072.35,2080.73,2001.58,1898.58,1881.21,1930.35,1982.42,1933.96,1962.83,1970.31,1962.83,1933.96,1982.42,1930.35,1881.21,1898.58,2001.58,2080.73,2072.35,2066.35,2018.67,2008.46,2023.54,2128.38,2086.63,2007.75,2016.02,2083.92,2145.44,2213.81,2258.25,2287.44,2287.94,2335.58,2395.1,2391.88,2460.15,2476.65,2411.5,2405.88,2426.5,2490.88,2571.81,2658.15,2784.71,2862.67,2827.29,2824,2852.4,2943.69,2998.65,3129.85,3178.1,3282.96,3417.06,3556.4,3739.92,3890.94,3937.88,3960.42,4023.46,3989.44,4095.5,4109.38,4244.48,4450.75,4638.33,4717.25,4911.35,5030.38,5118.35,5248.6,5446.83,5482.96,5581.15,5562.48,5617.63,5700.92,5812.06,6025.6,6286.9,6478.79,6700.71,6805.94,6891.73,6889.06,6879.17,6958.88,7222.08,7469.69,7594.85,7544.81,7498.92,7541.92,7529.88,7606.1,7782.85,7804.15,7763.52,7783.52}

This is a dataset of Z values as function of degree, so I make a table for the fit.

datatofit=Transpose[{Table[n*Pi/180, {n, 0, 180}], data}]

Now I define the function I want to fit with.

orbital00[t_]:= 1+i*LegendreP[2,Cos[t]]
orbital22[t_]:= 1+j*LegendreP[2,Cos[t]]
orbital02[t_]:= 1+k*LegendreP[2,Cos[t]]
intensity[t_]:=a*((1-c^2)*orbital00[t]+c^2*orbital22[t]*WignerD[{2, 0, 0}, t]^2 +2*c*Sqrt[1-c^2]*orbital02[t]*WignerD[{2, 0, 0}, t])

Alright let's try to fit. And I'm helping with coefficients very close to what I am looking for.

NonlinearModelFit[datatofit,intensity[t], {{a, 3600}, {c, .45}, {i, -.3}, {j, 1}, {k, .4}}, t]

I get an error

The function value (some list) is not a list of real numbers with dimensions {181} at {a,c,i,j,k} = {-215905.,-53.542,151.348,302.937,136.279}

Now, I insists, these coefficients are actually very close to the ones I am looking for, as if the function intensity[t] is plotted while replacing variables by the values provided, it is extremely close.

When trying to fit the data with a function A+B*Cos[t], it works (but that's not the function I am looking for).

I don't understand the error I have, and why I have it. I understand the function is complex but I am giving good starting parameters.

Origin 8 manages to do it quite easily.

$\endgroup$
10
$\begingroup$

You need to use a constraint for c to avoid complex numbers.

fit = NonlinearModelFit[
  datatofit, {intensity[t, a, c, i, j, k], 
   c^2 < 1}, {{a, 3600}, {c, .45}, {i, -.3}, {j, 1}, {k, .4}}, t]

Show[ListPlot[datatofit], Plot[fit[t], {t, 0, 180}]]

enter image description here

Btw, as a matter of style, I would define all parameters as variables. That way you won't run into troubles if you assign a value to any of these parameters at some point.

$\endgroup$
  • $\begingroup$ Thanks, I never had to use a constrain on the variables. And your suggestion on the variables is a good idea that I will try to keep in mind. $\endgroup$ – A postdoc Apr 7 '17 at 0:00
8
$\begingroup$

There are some "identifiability" issues with fitting your model parameters with the data that you have. First, your data set is symmetric in that the first half of the data is identical to the second half so if any standard errors are to be believed, just half of the data needs to be used.

@Felix has given you the approach for getting NonlinearModelFit to converge. However, when looking at nlm["CorrelationMatrix"] one sees very large correlations (with some being suspiciously close to 1.0).

$$\left( \begin{array}{ccccc} 1. & -0.752256 & -0.103413 & 0.99936 & -0.927615 \\ -0.752256 & 1. & -0.441861 & -0.730288 & 0.695337 \\ -0.103413 & -0.441861 & 1. & -0.135651 & -0.0770068 \\ 0.99936 & -0.730288 & -0.135651 & 1. & -0.926308 \\ -0.927615 & 0.695337 & -0.0770068 & -0.926308 & 1. \\ \end{array} \right)$$

This is strongly suggestive that the model is overparameterized. Consider expanding and collecting terms in the function being fitted:

FullSimplify[Expand[intensity[t]]]

$$\frac{1}{256} a \left(3 \cos (2 t) \left(c^2 (-64 i+39 j+32)+64 \sqrt{1-c^2} c (k+2)+64 i\right)+2 c \left(8 \sqrt{1-c^2} (11 k+8)+c (-32 i+29 j-84)\right)+9 c \left(2 \cos (4 t) \left(8 \sqrt{1-c^2} k+3 c j+4 c\right)+3 c j \cos (6 t)\right)+64 (i+4)\right)$$

This is of the form

$$\text{a0}+\text{a1} \cos(2t)+\text{a2} \cos(4t)+\text{a3} \cos(6t)$$

We see that there are just 4 parameters needed to characterize the function rather than 5. We can now perform a fit on half of the data:

nlm = NonlinearModelFit[datatofit[[Range[1, 90], All]], a0 + a1 Cos[2 t] + a2 Cos[4 t] + a3 Cos[6 t], {a0, a1, a2, a3}, t];
nlm["BestFitParameters"]
(* {a0 -> 4094.8464407396978, a1 -> 2737.6662380018993, a2 -> 724.3603397546482, a3 -> 154.1540954251538} *)
nlm["CorrelationMatrix"] // MatrixForm 

$$\left( \begin{array}{cccc} 1. & -0.0157271 & 0.000698895 & -0.0157271 \\ -0.0157271 & 1. & -0.0222359 & 0.00074129 \\ 0.000698895 & -0.0222359 & 1. & -0.0222359 \\ -0.0157271 & 0.00074129 & -0.0222359 & 1. \\ \end{array} \right)$$

That correlation matrix looks close to ideal.

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