Re: Data fitting with fmins

[Andreas Stahel]

Tweety <ftjp <at> gmx.de> writes:

> 
> Dear group,
> 
> I am trying to fit measured data with a function using this example:
>
http://fweb.wallawalla.edu/class-wiki/index.php/How_to_use_Octave_to_Fit_an_Arbitrary_Function_with_fmins
> 
> Problem is, the fitting parameters all remain 1 except for the last one 
> so my fitting function values all become INF.
> In my understanding, a more complex function may require more time to 
> converge but may not cause a fundamental issue? I also tried a "tanh" 
> function, but could not get the files to work.
> 
> What am I missing? Do I need to run some kind of update (similar to 
> "texhash"?) so that octave knows the "model.m" file exists? It is 
> located in the same folder as the data and the fitting file. Thanks for 
> your remarks!
> 
> Jan
> 
> My model.m looks like this:
> 
> function sum_square_errors = model(p,x,y)
> y_trial = p(1) + p(2)*exp(x.*p(3)) + p(4)*exp(x.*p(5));
> difference = y-y_trial;
> sum_square_errors = sum(difference.^2);
> 
> The file which I am trying to fit my data with looks like this:
> 
> clear; clc;
> E = dlmread('data.txt');
> x = E(:,1);
> y_data = E(:,2);
> p0  = [1 1 1 1 1];
> p = fmins('model', p0, [], [], x, y_data);
> y_fit = p(1) + p(2)*e.^(x*p(3)) + p(4)*e.^(x*p(5));
> figure(1)
> plot(x,y_data,'.', x,y_fit)
> 
> The data.txt file looks like this (typical saturation curve):
> 1329    0.1480062623
> 3981    0.4818839118
> 6711    0.691154131
> 9363    0.805889982
> 12093    0.8781710629
> 14745    0.9218823009
> 17475    0.9490571055
> 20127    0.9649592347
> 22857    0.975116832
> 25509    0.9815304812
> 28161    0.9857597316
> 30891    0.9887945353
> 33543    0.9910305769
> 36273    0.9927147068
> 38925    0.9939278847
> 41655    0.9951150653
> 44307    0.9961515922
> 47037    0.9970010438
> 49689    0.9977598722
> 52341    0.9983067517
> 55071    0.9987808017
> 57723    0.9990164135
> 60453    0.9992141382
> 63105    0.9993880715
> 65835    0.9995269449
> 68487    0.9996117692
> 71217    0.9996160582
> 73869    0.9995567346
> 76599    0.9994854522
> 
> Octave 3.8.1 on LinuxMint
> 
Dear Jan

You approch is fine, but the data does not fit well to a double exponential.
May I suggest to use the funtion leasqr() from the optimization package.
Find illustative example for this type of problems in my lectures notes 
in Section 2.2.13 on nonlinear regression, startin on page 167 in
https://web.ti.bfh.ch/~sha1/Labs/PWF/Documentation/OctaveAtBFH.pdf

The code below might illustrate the problem.

Hope this helps

Andreas
============
clear; 
%clc;
E = dlmread('data.txt');
x = E(:,1);
y_data = E(:,2);

if 0
  p0  = [1 1 1 1 1];
  p = fmins('model', p0, [], [], x, y_data)
  y_fit = p(1) + p(2)*e.^(x*p(3)) + p(4)*e.^(x*p(5));
  figure(1)
  plot(x,y_data,'.', x,y_fit)
endif



function res = f1(x,p)
  res = p(1) + p(2)*exp(p(3)*x);
endfunction
[fr,p,kvg,iter,corp,covp] = leasqr(x,y_data,[1,-1,-1e-4],'f1');
parameters = [p';sqrt(diag(covp))']
  
y_fit = f1(x,p);

figure(1)
plot(x,y_data,'.',x,y_fit)
xlabel('x'); ylabel('values')

figure(2)
plot(x,y_data-y_fit)
xlabel('x'); ylabel('difference')

%% now with two exponential, which will fail miserably
function res = f2(x,p)
  res = p(1) + p(2)*exp(p(3)*x) + p(4)*exp(p(4)*x);
endfunction

fr,p,kvg,iter,corp,covp] = leasqr(x,y_data,[1,-1,-1e-4,0,0],'f2');
parameters = [p';sqrt(diag(covp))']
  
y_fit = f2(x,p);
  
figure(3)
plot(x,y_data,'.',x,y_fit)
xlabel('x'); ylabel('values')

figure(4)
plot(x,y_data-y_fit)
xlabel('x'); ylabel('difference')