# How can i find the optimal parameters that minimize SSE?

12 views (last 30 days)
Riccardo Bartoloni on 30 Apr 2019
Answered: Alan Weiss on 2 May 2019
Hi, i am working for my university on this code and i have to minimize the last line, the variable "SQ_TOT". The solution of this exercise is the vector "Param" (that give me minimum SSE --> SQ_TOT= 3.1082e+05) but i don't know how to find this vector. How can i find the optimal parameters that minimize SQ_TOT? Thank you all!
clc
clear
t=0.33333:0.33333: 76.33333;
Param=[0.5254010 11.4402640 29861.1830810 2.0010250 -0.051606000 -0.057373000 0.036517000 -0.099394000 1.647586000 81.195773000 1.422753000 -0.078617000 0.036363000 0.106596000 0.128553000]
%Param=ones(1,15);
DUMqA=[Param(5:8)];
expbeta=exp(DUMqA*DUM);
t1=t.^Param(4);
t2=t.^Param(4);
t2(length(t))=[];
t2=[0 t2];
tc=t1-t2;
lt=length(t1);
cohort=zeros(1,lt);
for i=1:length(t1);
summ=[zeros(1,i-1) expbeta(i:length(t1)).*tc(1:length(t)+1-i)];
cohort=[cohort;summ];
end;
cohort(1,:)=[];
B_A=[];
HEY_A=[];
F_A=[];
for k=1:lt;
for i=1:lt;
summ=sum(cohort(k,1:(lt+1-i)));
B_A=[summ B_A];
end;
HEY_A=[HEY_A;B_A];
B_A=[];
end;
F_A=[];
for k=1:lt;
for i=1:lt;
summ=(1-(Param(3)/(Param(3)+HEY_A(k,i)))^Param(2))*(1-Param(1));
B_A=[B_A summ];
end;
F_A=[F_A; B_A];
B_A=[];
end;
DUMqR=[Param(12:15)];
expbeta_R=exp(DUMqR*DUM);
t1R=t.^Param(11);
t2R=t.^Param(11);
t2R(length(t))=[];
t2R=[0 t2R];
tc_R=t1R-t2R;
lt=length(t1);
cohort_R=[];
for i=1:length(t1);
summ=[zeros(1,i-1) expbeta_R(i:length(t1)).*tc_R(1:length(t)+1-i)];
cohort_R=[cohort_R;summ];
end;
cohort_R(1,:)=[];
B_R=[];
HEY_R=[];
for k=1:lt-1;
for i=1:lt;
summ=sum(cohort_R(k,1:(lt+1-i)));
B_R=[summ B_R];
end;
HEY_R=[HEY_R;B_R];
B_R=[];
end;
F_R=[];
for k=1:lt-1;
for i=1:lt;
summ=(1-(Param(10)/(Param(10)+HEY_R(k,i)))^Param(9));
B_R=[B_R summ];
end;
F_R=[F_R; B_R];
B_R=[];
end;
HouseCha1=housed;
housed(229)=[];
HouseCha2=[0 housed];
ChaninPop=HouseCha1-HouseCha2;
F_AZ=[zeros(lt,1) F_A];
EXP_A=[];
EXP_AA=[];
EC_R=zeros(1,lt);
EC=[];
F_RZ=[zeros(228,1) F_R];
EXP_R=[];
EXP_RR=[];
EC_RVERO=[];
EC_A=[zeros(lt,1)];
EXP_AA=[];
EXP_A=[];
EXP_RR=[];
EC=[];
EC_R=[zeros(lt,1)];
EC_AB=[0];
EC_TRAN=[];
EC_RTRAN=[];
for i=1:lt;
for k=1:lt-1;
EXP_AATRAN=(F_A(k,i)-F_AZ(k,i))*(ChaninPop(k)+EC_R(k));
EXP_A=[EXP_A EXP_AATRAN];
ECTRAN=sum(EXP_A(:));
EXP_RRTRAN=(F_R(k,i)-F_RZ(k,i))*EC_A(k);
EXP_R=[EXP_R EXP_RRTRAN];
EC_RVEROTRAN=sum(EXP_R(:));
end;
EXP_AA=[EXP_AA (EXP_A)'];
EXP_A=[];
EC_TRAN=[EC_TRAN ECTRAN];
EC_A=[EC_TRAN zeros(1,lt)];
EC=[EC ECTRAN];
EXP_RR=[EXP_RR (EXP_R)'];
EXP_R=[];
EC_RTRAN=[EC_RTRAN EC_RVEROTRAN];
EC_R=[0 EC_RTRAN zeros(1,lt)];
EC_RVERO=[EC_RVERO EC_RVEROTRAN];
end
EC3m_R=[0 0];
for i=3:lt;
summ=sum(EC_RVERO(i-2:i));
EC3m_R=[EC3m_R summ];
end;
EC3m=[EC(1) sum(EC(1:2))];
for i=3:lt;
summ=sum(EC(i-2:i));
EC3m=[EC3m summ];
end;
DIVEC3=DIV_.*EC3m;
DIV_A2=(DIV_A).^2;
DIV=LOSS~=0;
DIV_L=DIV.*EC3m_R;
DIV_LOSS=(DIV_L-LOSS)
DIV_LOSS2=(DIV_LOSS).^2
SQ_LOSS=sum(DIV_LOSS2(1,:));
EXP_E=EC-EC_RVERO;
EXP_END=[];
for i=1:lt;
summ=sum(EXP_E(1:i));
EXP_END=[EXP_END summ];
end;
DIV1=END~=0;
DIV2=DIV1-DIV_;
DIV_E=DIV2.*EXP_END;
END_DIV=DIV2.*END;
DIV_E(1)=EXP_END(1)
SQ_EN=[(EXP_END(1))^2];
pos=find(DIV2==1);
for i=4:3:pos(length(pos));
DIV_END=((DIV_E(i)-DIV_E(i-3))-(END_DIV(i)-END_DIV(i-3)))^2;
SQ_EN=[SQ_EN DIV_END];
end;
SQ_END=sum(SQ_EN(1,:));
SQ_END=SQ_END+(END(length(END))-EXP_END(length(EXP_END)))^2;

Alan Weiss on 1 May 2019
To solve an optimization problem you have to force your problem into the form required by optimization solvers. Sorry, that's just the way it is. You have to decide first which variables you can optimize over. These are usually called "decision variables", and are typically called x = [x(1),x(2),...,x(n)], where n is the number of decision variables.
Then you write your objective function, the thing to minimize, as a function of x. See, for example, Optimizing Nonlinear Functions.
Then you guess a starting value of x, say x0, and call an optimization solver, such as
sol = fminsearch(@objective,x0)
If you have constraints on your decision variables, such as they must all be positive, then you need an Optimization Toolbox™ license, and use an appropriate solver.
Good luck,
Alan Weiss
MATLAB mathematical toolbox documentation

Riccardo Bartoloni on 2 May 2019
Have i to do something like that? But matlab say "Not enough input arguments. Error in Myobjectivefunction (line 6)
DUMqA=[Param(5:8)];"
function [SQ_TOT] = Myobjectivefunction(Param)
t=0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333:0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333:76.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333;
%Param=[0.5254010 11.4402640 29861.1830810 2.0010250 -0.051606000 -0.057373000 0.036517000 -0.099394000 1.647586000 81.195773000 1.422753000 -0.078617000 0.036363000 0.106596000 0.128553000];
%Param=Param(1:15);
DUMqA=[Param(5:8)];
expbeta=exp(DUMqA*DUM);
t1=t.^Param(4);
t2=t.^Param(4);
t2(length(t))=[];
t2=[0 t2];
tc=t1-t2;
lt=length(t);
cohort=zeros(1,lt);
for i=1:length(t1);
summ=[zeros(1,i-1) expbeta(i:lt).*tc(1:lt+1-i)];
cohort=[cohort;summ];
end;
cohort(1,:)=[];
HEY_A=cumsum(cohort,2);
F_A=(1-(Param(3)./(Param(3)+HEY_A)).^Param(2)).*(1-Param(1));
%%%%Retention Model%%%%%
DUMqR=[Param(12:15)];
expbeta_R=exp(DUMqR*DUM);
t1R=t.^Param(11);
t2R=t.^Param(11);
t2R(length(t))=[];
t2R=[0 t2R];
tc_R=t1R-t2R;
cohort_R=[];
%cohort_R= matrice riga 1000
for i=1:length(t1);
summ=[zeros(1,i-1) expbeta_R(i:length(t1)).*tc_R(1:length(t)+1-i)];
cohort_R=[cohort_R;summ];
end;
cohort_R(1,:)=[];
HEY_R=cumsum(cohort_R,2);
F_R=(1-(Param(10)./(Param(10)+HEY_R)).^Param(9));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
HouseCha1=housed;
%HouseCha1= riga 19
housed(229)=[];
HouseCha2=[0 housed];
ChaninPop=HouseCha1-HouseCha2;
F_AZ=[zeros(lt,1) F_A];
EXP_A=[];
EXP_AA=[];
EC_RVERO=zeros(1,lt);
EC=[];
%EC= riga 976 oppure colonna IA e riga 1712
F_RZ=[zeros(228,1) F_R];
F_AZ(:,end)=[];
F_RZ(:,end)=[];
ChaninPopR=ChaninPop(1:end-1);
EC_RVERO=[0];
for i=1:lt
EXP_AA=(F_A-F_AZ).*(ChaninPop+EC_RVERO)';
EC=sum(EXP_AA);
EC(lt)=[];
EXP_RR=(F_R-F_RZ).*(EC)';
EC_RVERO=sum(EXP_RR);
EC_RVERO=[0 EC_RVERO];
EC_RVERO(end)=[];
end
EC=[EC sum(EXP_AA(:,end))];
EC_RVERO(1)=[];
EC_RVERO=[EC_RVERO sum(EXP_RR(:,end))];
EC3m_R=[0 0];
for i=3:lt;
summ=sum(EC_RVERO(i-2:i));
EC3m_R=[EC3m_R summ];
end;
EC3m=[EC(1) sum(EC(1:2))];
for i=3:lt;
summ=sum(EC(i-2:i));
EC3m=[EC3m summ];
end;
DIVEC3=DIV_.*EC3m;
DIV_A2=(DIV_A).^2;
DIV=LOSS~=0;
DIV_L=DIV.*EC3m_R;
DIV_LOSS=(DIV_L-LOSS);
DIV_LOSS2=(DIV_LOSS).^2
SQ_LOSS=sum(DIV_LOSS2(1,:));
EXP_E=EC-EC_RVERO;
EXP_END=cumsum(EXP_E);
DIV1=END~=0;
DIV2=DIV1-DIV_;
DIV_E=DIV2.*EXP_END;
END_DIV=DIV2.*END;
DIV_E(1)=EXP_END(1)
SQ_EN=[(EXP_END(1))^2];
pos=find(DIV2==1);
for i=4:3:pos(length(pos));
DIV_END=((DIV_E(i)-DIV_E(i-3))-(END_DIV(i)-END_DIV(i-3)))^2;
SQ_EN=[SQ_EN DIV_END];
end;
SQ_END=sum(SQ_EN(1,:));
SQ_END=SQ_END+(END(length(END))-EXP_END(length(EXP_END)))^2;
end
Riccardo Bartoloni on 2 May 2019
Sorry now it work thank you very much to you all!!
Riccardo Bartoloni on 2 May 2019
My last problem is to find the minimum. "fminsearch" arrive at a minimum of 500000 but the real minimum is more or less 180000. How can i find this minimum? This is my code. I put the data in only one file as you said.
function [SQ_TOT] = Myobjectivefunction(Param)
t=0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333:0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333:76.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333;
%Param=[0.5254010 11.4402640 29861.1830810 2.0010250 -0.051606000 -0.057373000 0.036517000 -0.099394000 1.647586000 81.195773000 1.422753000 -0.078617000 0.036363000 0.106596000 0.128553000];
%Param=Param(1:15);
DUMqA=[Param(5:8)];
expbeta=exp(DUMqA*DUM);
t1=t.^Param(4);
t2=t.^Param(4);
t2(length(t))=[];
t2=[0 t2];
tc=t1-t2;
lt=length(t);
cohort=zeros(1,lt);
for i=1:length(t1);
summ=[zeros(1,i-1) expbeta(i:lt).*tc(1:lt+1-i)];
cohort=[cohort;summ];
end;
cohort(1,:)=[];
HEY_A=cumsum(cohort,2);
F_A=(1-(Param(3)./(Param(3)+HEY_A)).^Param(2)).*(1-Param(1));
%%%%Retention Model%%%%%
DUMqR=[Param(12:15)];
expbeta_R=exp(DUMqR*DUM);
t1R=t.^Param(11);
t2R=t.^Param(11);
t2R(length(t))=[];
t2R=[0 t2R];
tc_R=t1R-t2R;
cohort_R=[];
%cohort_R= matrice riga 1000
for i=1:length(t1);
summ=[zeros(1,i-1) expbeta_R(i:length(t1)).*tc_R(1:length(t)+1-i)];
cohort_R=[cohort_R;summ];
end;
cohort_R(1,:)=[];
HEY_R=cumsum(cohort_R,2);
F_R=(1-(Param(10)./(Param(10)+HEY_R)).^Param(9));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
HouseCha1=housed;
%HouseCha1= riga 19
housed(229)=[];
HouseCha2=[0 housed];
ChaninPop=HouseCha1-HouseCha2;
F_AZ=[zeros(lt,1) F_A];
EXP_A=[];
EXP_AA=[];
EC_RVERO=zeros(1,lt);
EC=[];
%EC= riga 976 oppure colonna IA e riga 1712
F_RZ=[zeros(228,1) F_R];
F_AZ(:,end)=[];
F_RZ(:,end)=[];
ChaninPopR=ChaninPop(1:end-1);
EC_RVERO=[0];
for i=1:lt
EXP_AA=(F_A-F_AZ).*(ChaninPop+EC_RVERO)';
EC=sum(EXP_AA);
EC(lt)=[];
EXP_RR=(F_R-F_RZ).*(EC)';
EC_RVERO=sum(EXP_RR);
EC_RVERO=[0 EC_RVERO];
EC_RVERO(end)=[];
end
EC=[EC sum(EXP_AA(:,end))];
EC_RVERO(1)=[];
EC_RVERO=[EC_RVERO sum(EXP_RR(:,end))];
EC3m_R=[0 0];
for i=3:lt;
summ=sum(EC_RVERO(i-2:i));
EC3m_R=[EC3m_R summ];
end;
EC3m=[EC(1) sum(EC(1:2))];
for i=3:lt;
summ=sum(EC(i-2:i));
EC3m=[EC3m summ];
end;
DIVEC3=DIV_.*EC3m;
DIV_A2=(DIV_A).^2;
DIV=LOSS~=0;
DIV_L=DIV.*EC3m_R;
DIV_LOSS=(DIV_L-LOSS);
DIV_LOSS2=(DIV_LOSS).^2
SQ_LOSS=sum(DIV_LOSS2(1,:));
EXP_E=EC-EC_RVERO;
EXP_END=cumsum(EXP_E);
DIV1=END~=0;
DIV2=DIV1-DIV_;
DIV_E=DIV2.*EXP_END;
END_DIV=DIV2.*END;
DIV_E(1)=EXP_END(1)
SQ_EN=[(EXP_END(1))^2];
pos=find(DIV2==1);
for i=4:3:pos(length(pos));
DIV_END=((DIV_E(i)-DIV_E(i-3))-(END_DIV(i)-END_DIV(i-3)))^2;
SQ_EN=[SQ_EN DIV_END];
end;
SQ_END=sum(SQ_EN(1,:));
SQ_END=SQ_END+(END(length(END))-EXP_END(length(EXP_END)))^2;
end

Alan Weiss on 2 May 2019
I think that you did not understand the comment about not loading the data within your objective function. You should load it once before you call the objective, and just pass in the values, like this:
fun = @(Param)Myobjectivefunction(Param,DUM);
best = fminsearch(fun,x0) % Of course, you have to define x0 first
% Separately, you change your function like so:
function SQ_TOT = Myobjectivefunction(Param,DUM)
% The rest of your function goes here, without the load and DUM = statements
end
To search for a smaller minimum, you can try starting you search from a variety of initial points x0. See Nonlinear Data-Fitting. Also, fminsearch is not really good when you have so many parameters. Maybe you can use lsqnonlin?
Alan Weiss
MATLAB mathematical toolbox documentation