How can I fit a function to my experimental data?

Hi, apologies if something has been asked similar before but I couldn't find anything that was that appropiate to my problem.

I have some experimental data that I am trying to fit to the function y = C1 + e^(dt) * (C2 cos(ωt)+ C3 sin(ωt)) where C1, C2, C3, ω and d are all constants which I need to validate my model.

I have tried to follow the lsqcurvefit example where I have done:

>> Data=...
[0,1.48737300000000;0.000977000000000000,1.48615200000000...........0.125000000000000,1.48325300000000];

t=Data(:,1);

>> y=Data(:,2);

>> % axis([0 2 -0.5 6])

>> % axis([1.47 1.5 0 0.125])

>> % hold on

>> plot(t,y,'ro')

>> title('Data points')

>> %hold off

>> F=@(x,xdata)x(1)+exp(x(2)*xdata)*(x(3)*cos(x(4)*xdata)+x(5)*sin(x(4)*xdata));

>> x0=[1 100 1 1 1];

>> [x,resnorm,~,exitflag,output] = lsqcurvefit(F,x0,t,y)

hold on

plot(t,F(x,t))

hold off

Unfortunately I am getting the following error:

??? Error using ==> mtimes
Inner matrix dimensions must agree.

Error in ==>
@(x,xdata)x(1)+exp(x(2)*xdata)*(x(3)*cos(x(4)*xdata)+x(5)*sin(x(4)*xdata))

Error in ==> lsqcurvefit at 209
            initVals.F =
            feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Caused by:
    Failure in initial user-supplied objective function
    evaluation. LSQCURVEFIT cannot continue.

Could anyone tell me where I am going wrong? And am I using the right function? I have seen fminsearch has been used on similar problems but couldn't find any examples to follow through. I should apologise in advance for my lack of experience in Matlab. If anyone could help I would be hugely appreciative!

Thank you

Answer 1

[Star Strider]

You need to vectorize your objective function. Experiment with this version:

F=@(x,xdata) x(1)+exp(x(2).*xdata).*(x(3).*cos(x(4).*xdata)+x(5).*sin(x(4).*xdata));

The dot operator modifies multiplication (.*), division (./), and exponentiation (.^) for element-by-element operations rather than matrix or vector operations.

[Joshua]

Thank you! I am definitely getting numbers out now but unfortunately it seems to be mainly straight lines/curves and I was hoping for a sort of sinusoid. Am I better off using a different function?

[John D'Errico]

You have not attached enough data to know what you are doing wrong. If the model is not fitting your data, then the answer is usually one of:

1. Your model is NOT a good choice for that data.
2. You have chosen poor starting values.

Note that my fminspleas will greatly help convergence on this model, as it reduces the number of starting values that you must provide. fminspleas is found for download from the file exchange.

[Star Strider]

I’m glad that helped. Without knowing more about your data, I can’t tell if your function is appropriate to it.

We need values for your Data array. I notice that there is very little change in both columns between rows for the values you posted. Either you have a huge array (and if that is the situation, I invite you to downsample it before you post it), or you may not have enough data to uniquely determine the frequencies denoted by x(4).

[Star Strider]

I‘ll comment here because I don’t feel like scrolling all the way to the end of this page.

This was a challenge because of the noise in your signal. That made it difficult to calculate the period of the sine wave. Using the second set of data you provided:

t = [  0.400391  0.402344  0.404297  0.406250 ...
y = [  1.473945  1.473030  1.469368  1.491188 ...

I got a reasonably good fit:

F(t) = 1.452 + exp(-1.924 * t) * (-0.032*cos(23.747 * t) - 0.068*sin(23.747 * t))

from this code:

% Estimate y-offset and period (t-units):
yofs = mean(y);
yz = y - yofs;
yprd = max(t( (yz .* circshift(yz, [0  1])) < 0))/2;

F3 = @(x,xdata) x(1)+exp(x(2)*xdata).*(x(3).*cos(x(4).*xdata)+x(5).*sin(x(4).*xdata));
opts = statset('MaxIter',10000);
B3 = nlinfit(t, y, F3, [yofs; -.01; 0.5; 2.5; 0.1], opts)

s = sprintf('F(\\beta,t) = %.3f + exp(%.3f * t) * (%.3f*cos(%.3f * t) %c %.3f*sin(%.3f * t))\n', B3(1:4), char(45-sign(B3(5))-1), abs(B3(3)), B3(4));

figure(1)
plot(t, y, '-b')
hold on
plot(t, F3(B3,t),'-r', 'LineWidth',2)
hold off
text(0.45,1.51,s,'FontSize',8)
grid

You mentioned fminsearch, so here’s a version using it:

% Using fminsearch: 
fcn = @(b) sum((F3(b,t) - y).^2);                           % Least-Squares cost function
opts4 = optimset('MaxFunEvals',50000, 'MaxIter',10000);
B4 = fminsearch(fcn, [yofs; -2; 0.5; 25.; 0.1])             % Minimise Least-Squares

[Joshua]

Youve been a great help! I really appreciate it!

[Star Strider]

My pleasure!

Answer 2

[Joshua]

I think that data ia a bad example. I have been using the following:

t=[0.400391
0.402344
0.404297
0.40625
0.408203
0.410156
0.412109
0.414062
0.416016
0.417969
0.419922
0.421875
0.423828
0.425781
0.427734
0.429687
0.431641
0.433594
0.435547
0.4375
0.439453
0.441406
0.443359
0.445312
0.447266
0.449219
0.451172
0.453125
0.455078
0.457031
0.458984
0.460937
0.462891
0.464844
0.466797
0.46875
0.470703
0.472656
0.474609
0.476562
0.478516
0.480469
0.482422
0.484375
0.486328
0.488281
0.490234
0.492187
0.494141
0.496094
0.498047
0.5
0.501953
0.503906
0.505859
0.507812
0.509766
0.511719
0.513672
0.515625
0.517578
0.519531
0.521484
0.523437
0.525391
0.527344
0.529297
0.53125
0.533203
0.535156
0.537109
0.539062
0.541016
0.542969
0.544922
0.546875
0.548828
0.550781
0.552734
0.554687
0.556641
0.558594
0.560547
0.5625
0.564453
0.566406
0.568359
0.570312
0.572266
0.574219
0.576172
0.578125
0.580078
0.582031
0.583984
0.585937
0.587891
0.589844
0.591797
0.59375
0.595703
0.597656
0.599609
0.601562
0.603516
0.605469
0.607422
0.609375
0.611328
0.613281
0.615234
0.617187
0.619141
0.621094
0.623047
0.625
0.626953
0.628906
0.630859
0.632812
0.634766
0.636719
0.638672
0.640625
0.642578
0.644531
0.646484
0.648437
0.650391
0.652344
0.654297
0.65625
0.658203
0.660156
0.662109
0.664062
0.666016
0.667969
0.669922
0.671875
0.673828
0.675781
0.677734
0.679687
0.681641
0.683594
0.685547
0.6875
0.689453
0.691406
0.693359
0.695312
0.697266
0.699219
0.701172
0.703125
0.705078
0.707031
0.708984
0.710937
0.712891
0.714844
0.716797
0.71875
0.720703
0.722656
0.724609
0.726562
0.728516
0.730469
0.732422
0.734375
0.736328
0.738281
0.740234
0.742187
0.744141
0.746094
0.748047
0.75
0.751953
0.753906
0.755859
0.757812
0.759766
0.761719
0.763672
0.765625
0.767578
0.769531
0.771484
0.773437
0.775391
0.777344
0.779297
0.78125
0.783203
0.785156
0.787109
0.789062
0.791016
0.792969
0.794922
0.796875
0.798828
0.800781
0.802734
0.804687
0.806641
0.808594
0.810547
0.8125
0.814453
0.816406
0.818359
0.820312
0.822266
0.824219
0.826172
0.828125
0.830078
0.832031
0.833984
0.835937
0.837891
0.839844
0.841797
0.84375
0.845703
0.847656
0.849609
0.851562
0.853516
0.855469
0.857422
0.859375
0.861328
0.863281
0.865234
0.867187
0.869141
0.871094
0.873047
0.875
0.876953
0.878906
0.880859
0.882812
0.884766
0.886719
0.888672
0.890625
0.892578
0.894531
0.896484
0.898437
0.900391
0.902344
0.904297
0.90625
0.908203
];
y=[1.473945
1.47303
1.469368
1.491188
1.482338
1.48249
1.5066
1.49073
1.483558
1.47837
1.466011
1.468299
1.465858
1.482185
1.496223
1.495155
1.49073
1.481117
1.482185
1.472877
1.468452
1.476387
1.478828
1.475776
1.483101
1.496071
1.482643
1.479438
1.481117
1.465248
1.464942
1.469825
1.477607
1.48127
1.479438
1.483864
1.49012
1.47013
1.474556
1.477455
1.46067
1.473487
1.47013
1.468147
1.478065
1.456092
1.464332
1.466774
1.446632
1.470283
1.462043
1.45594
1.469062
1.450904
1.4654
1.457008
1.456092
1.471199
1.460365
1.468299
1.457923
1.455329
1.457008
1.441291
1.441291
1.444038
1.429999
1.416876
1.437018
1.425421
1.414588
1.433661
1.415045
1.418708
1.424353
1.404669
1.421759
1.424811
1.414282
1.434577
1.430152
1.428016
1.442054
1.440528
1.451057
1.441138
1.434424
1.445258
1.435645
1.431067
1.441596
1.441596
1.439918
1.442664
1.440528
1.437476
1.436103
1.43183
1.434119
1.433814
1.429084
1.435035
1.437171
1.43473
1.437476
1.438087
1.436408
1.433967
1.42649
1.425116
1.426032
1.426032
1.425879
1.424964
1.426795
1.423438
1.420691
1.42359
1.425116
1.426795
1.433509
1.431525
1.437171
1.443885
1.437476
1.450141
1.448005
1.441749
1.462654
1.454719
1.457618
1.472572
1.458686
1.467384
1.468452
1.459754
1.469062
1.465858
1.465705
1.46891
1.464485
1.465705
1.473487
1.467079
1.462043
1.472114
1.465858
1.467384
1.471504
1.463416
1.47303
1.477607
1.478065
1.479896
1.47425
1.472114
1.475776
1.474556
1.476081
1.480354
1.479286
1.480049
1.478065
1.475166
1.475471
1.472419
1.470741
1.471351
1.474861
1.476387
1.472572
1.476081
1.476081
1.467384
1.473945
1.470436
1.460822
1.469825
1.461433
1.457008
1.469825
1.452277
1.453193
1.473945
1.453346
1.452888
1.459754
1.433509
1.444495
1.441749
1.431067
1.457008
1.442054
1.43946
1.456397
1.427405
1.434119
1.445563
1.427863
1.443885
1.440986
1.438849
1.451667
1.432746
1.446326
1.450294
1.432135
1.4477
1.442207
1.433661
1.444648
1.437476
1.437781
1.442664
1.434119
1.44358
1.449836
1.432593
1.435035
1.438239
1.425574
1.431067
1.438544
1.433967
1.439155
1.438849
1.434119
1.437629
1.44007
1.445258
1.448005
1.443122
1.443885
1.447242
1.446021
1.445869
1.447547
1.449226
1.451362
1.452277
1.45304
1.45182
1.446632
1.447547
1.444953
1.443122
1.44358
1.444495
1.446174
1.450294
1.447089
1.440986
1.447242
1.449683
1.446021
1.445716
1.446021
1.4448
1.44007
1.442207
1.447242
];

[Joshua]

Sorry about the length! I have tried playing around with the constants but don't think this is going to cut it. I am trying to follow an experiment that has been completed previously that used fminsearch but I can't find any examples that seem particularly similar to mine which is why I originally went down this method.

Answer 3

[John D'Errico]

This is a good example of why I wrote and provided my SLM toolbox.

- You have no good physical model for the process.
- The relationship is a complicated one.

The first point is clear. That you seem willing to try a different model tells me that you have no model that makes sense from physical reasoning.

There are valid reasons to need a nonlinear curve fit of course. But if all you need is a smooth curve that replicates the essential shape of your data, then a spline model is right. You can use it to create a plot, or to predict values of the relationship at any point. Of course, splines don't extrapolate well, but there is no sane way to extrapolate this relationship anyway. And if you want to see a functional form that you can write down to put in a paper, then again you are out of luck with a spline model.

With 10 equally spaced knots, we get a curve that replicates the essential shape, but does a fair amount of smoothing.

slm = slmengine(t,y,'plot','on','knots',10,'reg','c');

With 50 equally spaced knots, it has the flexibility to chase after some of the wiggles in the curve.

slm = slmengine(t,y,'plot','on','knots',50,'reg','c');

Use of cross validation here to automatically smooth the curve based on the amount of noise it sees seems the right choice, so with 25 knots, the resulting curve was not that much different from what I see with 50 knots.

As it turns out, I needed little of the many capabilities of SLM to fit your data, but a spline model seems best to me.

[Star Strider]

I would agree, except that Joshua mentioned in his original post that he has a model for his data but that he was having problems fitting it, most likely because of the noise. Empirically determining the offset, amplitude and period from the data and using them as initial parameter estimates for x(1), x(3), and x(4) allowed a simple nonlinear regression to fit his model to the other parameters as well.