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
采纳的回答
Mathieu NOE
2023-7-19
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
3 个评论
Blanca Castells
2023-7-19
Alex Sha
2023-7-19
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
2023-7-19
as always, my pleasure !
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Least Squares 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
