3

I have tried to follow this tutorial to fit a curve to my dataset. The equation for the curve should be 

f(t) = log10((wpmcoeff./(t.^2)) + 
       ((1.038+3.*log(2.*pi.*1e6.*t)).*fpmcoeff./(t.^2))+(wfmcoeff./t) + 
       (ffmcoeff)+(rwfmcoeff.*t)).

I have created the following code: 

clock='atomicclockgpsworld.txt';
data=importdata(clock);

carrier=10e6;

sig=data(:,2);
t=data(:,1);
sigsq=log10(sig.^2);
fun = @(coeff)sseval(coeff,t,sigsq);
x0 = rand(5,1);
bestx = fminsearch(fun,x0);
wpmcoeff = bestx(1);
fpmcoeff = bestx(2);
wfmcoeff = bestx(3);
ffmcoeff = bestx(4);
rwfmcoeff = bestx(5);
yfit=log10((wpmcoeff./(t.^2))+((1.038+3.*log(2.*pi.*1e6.*t)).*fpmcoeff./(t.^2))+(wfmcoeff./t))+(ffmcoeff)+(rwfmcoeff.*t);
semilogx(t,sigsq,'x');
hold on
semilogx(t,yfit);
saveas(gcf,'fit','png');

and the corresponding function

function sse = sseval(coeff,t,sigsq)
    wpmcoeff = coeff(1);
    fpmcoeff = coeff(2);
    wfmcoeff = coeff(3);
    ffmcoeff = coeff(4);
    rwfmcoeff = coeff(5);
    sse = sum(sigsq - (log10((wpmcoeff./(t.^2))+((1.038+3.*log(2.*pi.*1e6.*t)).*fpmcoeff./(t.^2))+(wfmcoeff./t)+(ffmcoeff)+(rwfmcoeff.*t))));
end

But the fit produced is horrible (my y data should vary between approximately -20 to -22 but the fit produces a curve that reaches 1e59!). Can anyone suggest where I may be going wrong?

Current output vs data:

Link to the plot that is created.

Bethany Baxter
  • 31
  • 4

1 Answers1

3

In your function sseval the function is different then the one in your first script. In the function you take the log10 of the whole equation, while in the script log10 stops at (wfmcoeff./t):

first script:

log10((wpmcoeff./(t.^2))+((1.038+3.*log(2.*pi.*1e6.*t)).*fpmcoeff./(t.^2))+(wfmcoeff./t))+(ffmcoeff)+(rwfmcoeff.*t)

function:

log10((wpmcoeff./(t.^2))+((1.038+3.*log(2.*pi.*1e6.*t)).*fpmcoeff./(t.^2))+(wfmcoeff./t)+(ffmcoeff)+(rwfmcoeff.*t))

A second thing is that you didn't take the square of the difference. So in your function change the last line to

sse = sum((sigsq - (log10((wpmcoeff./(t.^2))+((1.038+3.*log(2.*pi.*1e6.*t)).*fpmcoeff./(t.^2))+(wfmcoeff./t))+(ffmcoeff)+(rwfmcoeff.*t))).^2);

Note: the fitting is sometimes good, sometimes it's really rubbish, so run the script a couple of times to get a good fitting. After a few times I got the following:

enter image description here

The parameters found are:

wpmcoeff: 10.898483535691309
wfmcoeff: 22.933722400252414
rwfmcoeff: 3.601059651996531e-05
fpmcoeff: -0.324473127267299
ffmcoeff: -21.497862719234053
ViG
  • 1,780
  • 1
  • 9
  • 13
  • Thank you so much @ViG! I really appreciate your help! My only issue now is that the whole function was actually supposed to be 'log'ed - I have done this in my code but can't reproduce a fit anywhere near as nice as the one you found.. mine looks like this after several attempts https://www.dropbox.com/s/pyms3zcizy3g937/fit.png?dl=0 – Bethany Baxter Feb 16 '18 at 15:53
  • 1
    The results I got when I log the whole function are even worse than yours. I don't know why it is that bad now.. – ViG Feb 19 '18 at 12:18