Optimization problem fitting arrays

Hi there!

I'm not very good at matlab and I'm having trouble with what I think is an easy optimization problem. I have an array of points (called alp) and I intend to fit it using an ecuation with 12 parameters. This ecuation leads to another array of points that should fit the original array (alp) by modifying the 12 parameters. Therefore I use fmincon trying to minimize the error calculated as the square sum of the original values subtracting the estimated values (SSE in statistics). I'm not getting any error but the error is not close to zero and the arrays do not fit.

Does anybody know what might be happening?

%load alp and T


initialGuess = [0,0,0,0,0,0,10,10,10,10,10,10,];
% Define the lower and upper bounds 
lowerBounds = [0,0,0,0,0,0,0,0,0,0,0,0];
upperBounds = [1,1,1,1,1,1,1000,1000,1000,1000,1000,1000,];

% Use fmincon to find the optimal parameters that minimize the error
options = optimoptions('fmincon', 'StepTolerance', 1e-14,'Display', 'iter', 'FiniteDifferenceStepSize', 1e-9);
paramsOpt = fmincon(@(params) objectiveFunction(params, T, alp), initialGuess, [], [], [], [], lowerBounds, upperBounds, [], options);

% Extract the optimal values for parameters
H1 = paramsOpt(1);
H2 = paramsOpt(2);
H3 = paramsOpt(3);
S1 = paramsOpt(4);
S2 = paramsOpt(5);
S3 = paramsOpt(6);
P1 = paramsOpt(7);
P2 = paramsOpt(8);
P3 = paramsOpt(9);
W1 = paramsOpt(10);
W2 = paramsOpt(11);
W3 = paramsOpt(12);



    Y1= H1*exp(-log(2)/S1^2*(log(1+2*S1*(T-P1)/W1.^2)));
    Y2= H2*exp(-log(2)/S2^2*(log(1+2*S2*(T-P2)/W2.^2)));
    Y3= H3*exp(-log(2)/S3^2*(log(1+2*S3*(T-P3)/W3.^2)));
    
   
    
    Y=Y1+Y2+Y3;
    
figure, clf
plot(T,alp)
hold on
plot(T,Y,'-o')

function SSE = objectiveFunction(params, T, alp)
   h1=params(1);
   h2=params(2);
   h3=params(3);
   s1=params(4);
   s2=params(5);
   s3=params(6); 
   p1=params(7);
   p2=params(8);
   p3=params(9);
   w1=params(10);
   w2=params(11);
   w3=params(12);
   
    y1 = h1*exp(-log(2)/s1^2*(log(1+2*s1*(T-p1)/w1.^2)));
    y2 = h2*exp(-log(2)/s3^2*(log(1+2*s2*(T-p2)/w2.^2)));
    y3 = h3*exp(-log(2)/s3^2*(log(1+2*s3*(T-p3)/w3.^2)));

    V = y1+y2+y3;
    %SSR = sum((V-mean(alp)).^2);
    SSE = sum((alp - V).^2);
    
    
    %Rsq=(SSR/(SSE+SSR))^(1/2);
end

Answer 1

[Mathieu NOE]

hello

this is the poor man's solution with fminsearch (as I don't have the optimization toolbox but you can easily use fmincon instead of fminsearch

I had to modify your fit model

your version :

y1 = h1*exp(-log(2)/s1^2*(log(1+2*s1*(T-p1)/w1.^2)));

my version

y1 = h1*exp(-log(2)/s1^2*(log(1+2*s1*((T-p1)/w1).^2)));

so you get a peak centered around x= p1 ( your original model does not show any peak)

load('alp.mat')
load('T.mat')

x = T;
y = alp;

% curve fit using fminsearch
% We would like to fit the function
f = @(a,b,c,d,e,f,g,h,k,l,m,n,x) a*exp(-log(2)/b^2*(log(1+2*b*((T-c)/d).^2))) + e*exp(-log(2)/f^2*(log(1+2*f*((T-g)/h).^2))) + k*exp(-log(2)/l^2*(log(1+2*l*((T-m)/n).^2)));
obj_fun = @(par) norm(f(par(1),par(2),par(3),par(4),par(5),par(6),par(7),par(8),par(9),par(10),par(11),par(12),x)-y);
C1_guess = [0.005 0.5 275 50 0.004 0.5 315 30 0.0015 0.5 500 300];
sol = fminsearch(obj_fun, C1_guess); % 
h1 = sol(1); % H1
s1 = sol(2); % S1
p1 = sol(3); % P1
w1 = sol(4); % W1
h2 = sol(5); % H2
s2 = sol(6); % S2
p2 = sol(7); % P2
w2 = sol(8); % W2
h3 = sol(9); % H3
s3 = sol(10); % S3
p3 = sol(11); % P3
w3 = sol(12); % W3

y_fit = f(h1,s1,p1,w1,h2,s2,p2,w2,h3,s3,p3,w3,x);
Rsquared = my_Rsquared_coeff(y,y_fit); % correlation coefficient

plot(x, y_fit, '-',x,y, 'r .-', 'MarkerSize', 15)

title(['Fit equation to data : R² =  ' num2str(Rsquared) ], 'FontSize', 20)
xlabel('x data', 'FontSize', 20)
ylabel('y data', 'FontSize', 20)
legend('fit ','data');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Rsquared = my_Rsquared_coeff(data,data_fit)
    % R2 correlation coefficient computation
    
    % The total sum of squares
    sum_of_squares = sum((data-mean(data)).^2);
    
    % The sum of squares of residuals, also called the residual sum of squares:
    sum_of_squares_of_residuals = sum((data-data_fit).^2);
    
    % definition of the coefficient of correlation is
    Rsquared = 1 - sum_of_squares_of_residuals/sum_of_squares;
    
    
end

[Blanca Castells]

Than you very much Mathieu!! Now it works perfectly!!

[Alex Sha]

As Mathieu pointed out, "(T-p1)/w1.^2" should be "((T-p1)/w1).^2", same as "(T-p2)/w2.^2" should be "((T-p2)/w2).^2" and "(T-p3)/w3.^2" should be "((T-p3)/w3).^2", the result will be:

Sum Squared Error (SSE): 1.68645694617025E-6
Root of Mean Square Error (RMSE): 8.19692135912861E-5
Correlation Coef. (R): 0.998583864294734
R-Square: 0.997169734029805

Parameter	Best Estimate      
---------	-------------      
h1       	0.00149375797077984
h2       	0.00494557714226805
h3       	0.00415875950017576
s1       	0.0013350334579832 
s2       	0.353521666078417  
s3       	0.610593297493903  
p1       	453.849315355541   
p2       	273.19404888193    
p3       	316.833648995382   
w1       	5827.85677112038   
w2       	54.4256102795608   
w3       	28.3759648901107 

[Mathieu NOE]

@Blanca Castells

as always, my pleasure !