Fitting procedure for cochlear mechanics

Hello, I would like to fit some data to an analytic function. The data is attached. My code is:

clear all;
close all;

%%Fit after Lesser and Barkeley
%%Load data
filename = 'daten.xls';
A = xlsread(filename);

xdata=100*A(:,1)    %x-values in mm
ydata=A(:,2)        %y-values in nm

%%Start parameters
f0 = 500 ;      % Frequence of the measurement

F0 = 0.1;     % N
m1 = 0.5;       % kg/m^2
k1 = 10.^10;    % N/m^3
k2 = 3;        % 1/m
c1 = 30000;    % N*s/m^3
c2 = 1.5;      % 1/m

%%Parameterization

a1 = [F0, m1, k1, k2, c1, c2];           % Startparameter für die Fitfunktion
lb = [1e-2, 1e-3, 1e8, 1e-3, 1e3, 1e-3]; % Boarders
ub = [1e0, 1e2, 1e12, 1e3, 1e7, 1e3];          

fun = @(a,x) a(1)./(2*pi*f0*sqrt((a(5).*exp(-a(6)*x)).^2+(a(2).*2*pi*f0-((a(3).*exp(-a(4).*x))/(2*pi*f0))).^2));
options = optimset('TolX',1e-10,'TolFun',1e-12);
% ,'TypicalX', [1e-3, 0.1, 1e14, 0.1, 1000, 0.1]
a = lsqcurvefit(fun,a1,xdata,ydata,lb,ub,options);

%%Fitting curve
xfit = linspace(0.1, 14, 1000);
yfit = a(1)./(2*pi*f0*sqrt((a(5).*exp(-a(6)*xfit)).^2+(a(2).*2*pi*f0-((a(3).*exp(-a(4).*xfit))/(2*pi*f0))).^2));

%%Plotten der Daten
hold all;
plot(xdata,ydata,'k*', 'LineWidth',2);
plot(xfit,yfit,'r', 'LineWidth',2);
xlabel('Distance in mm');
ylabel('Distance in m');

Unfortunetely I am not able to create a proper fit of the data set. The thing is, my startig parameters are for a normal human model. My model is 4 times the size. So the starting parameters are not optimal choice but the only one choice. Does anybody has an idea, how to fit the parameters properly. It is common that the root is not negative. Maybe i can implement that option somehow in the code?

Best regards Henrik

Answer 1

[Star Strider]

I worked on this for a while, going to my favourite global optimisation routine, the genetic algorithm, the Global Optimization Toolbox ga function. It is ‘derivative-free’, so searches the parameter space for the best parameter set rather than using a gradient-descent algorithm that is exquisitely sensitive to the initial parameter estimates. If you have access to the Global Optimization Toolbox, and since you know your model and what the parameters represent, experiment with my code to get a better result than I got. (I did an unconstrained optimization, since using your supplied constraints did not produce what I considered to be an acceptable fit. I commented-out the constrained ga call here, and kept it in for you to use as you wish.)

The ‘fitfun’ function is the fitness function ga uses to optimise its parameter estimates. It will also work for other global optimisation routines.

The Code —

filename = 'Henrik Schädlich daten.xls';
A = xlsread(filename);
xdata=100*A(:,1);       % x-values in mm
ydata=A(:,2);           % y-values in nm

f0 = 500 ;      % Frequence of the measurement

F0 = 0.1;       % N
m1 = 0.5;       % kg/m^2
k1 = 10E+10;    % N/m^3
k2 = 3;         % 1/m
c1 = 30000;     % N*s/m^3
c2 = 1.5;       % 1/m

% %%Parameterization

a1 = [F0, m1, k1, k2, c1, c2];           % Startparameter für die Fitfunktion
lb = [1e-2, 1e-3, 1e8, 1e-3, 1e3, 1e-3]; % Boarders
ub = [1e0, 1e2, 1e12, 1e3, 1e7, 1e3];          

fun = @(a,x) a(1)./(2*pi*f0*sqrt((a(5).*exp(-a(6)*x)).^2+(a(2).*2*pi*f0-((a(3).*exp(-a(4).*x))/(2*pi*f0))).^2));
fitfun = @(a) norm(ydata - fun(a,xdata));

options = optimoptions('ga', 'InitialPopulationMatrix',a1, 'MaxGenerations', 1E+5, 'FunctionTolerance',1E-8);
% [a, fval, exitflag] = ga(fitfun, 6, [], [], [], [], lb, ub, [], options);
[a, fval, exitflag] = ga(fitfun, 6, [], [], [], [], [], [], [], options);

% %%Fitting curve
xfit = linspace(0.1, 14, 1000);
yfit = a(1)./(2*pi*f0*sqrt((a(5).*exp(-a(6)*xfit)).^2+(a(2).*2*pi*f0-((a(3).*exp(-a(4).*xfit))/(2*pi*f0))).^2));

figure(1)
plot(xfit,yfit,'r', 'LineWidth',2);
hold on
plot(xdata,ydata,'kp', 'LineWidth',2, 'MarkerFaceColor','k');
hold off
xlabel('Distance in mm');
ylabel('Distance in m');
set(gca, 'YLim', [0  3.5E-7])

fprintf('a = [%23.15E, %23.15E, %23.15E, %23.15E, %23.15E, %23.15E]; Residual Norm = %23.15E\n\n', a, fval)

This parameter set seems to me to produce the ‘best’ fit, although not as good as I would like. The ga function will come up with different parameter sets if you let it explore:

a = [  7.260159508994168E+00,  -3.718242469417371E+00,   1.000000001416294E+11,   1.672499867222527E+00,   9.739036285122584E+00,  -6.062517640300058E-01]; Residual Norm =   2.465399858581893E-07

Others that seem to work:

a = [  6.327863132225293E+00,  -3.473492727621029E+00,   1.000000001735362E+11,   1.827669273786148E+00,  -9.230281839238965E+01,  -4.002542268418363E-01]; Residual Norm =   2.842624031350991E-07
a = [  8.260691541604469E+00,   5.439710355062475E+00,   1.826455300281449E+01,  -1.393791538277103E+00,   3.006655643029021E+04,   1.364162877863433E-01]; Residual Norm =   2.838486667836503E-07
a = [  2.071960544990774E+01,  -1.259075316947824E+01,   1.000000001813016E+11,   1.463195931523877E+00,  -1.925367345857832E+01,   2.202329442425367E+01]; Residual Norm =   3.188882864786421E-07
a = [  8.182160873658439E-01,   7.151731072216396E-01,   1.002935536031425E+01,  -1.266360643079992E+00,  -1.846892606921009E+02,   7.447706710989812E-01]; Residual Norm =   3.344567656769925E-07

My apologies for the delay in responding. Good luck in your research!