MATLAB Answers

0

How can i find the optimal parameters that minimize SSE?

Asked by Riccardo Bartoloni on 30 Apr 2019
Latest activity Answered by 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);
load('DUM.mat');%DUM is 4x229
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;
load('housed.mat');%housed is 1x229
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;
load('DishADD.mat'); %ADD is 1x229
DIV_=ADD~=0;
DIVEC3=DIV_.*EC3m;
DIV_A=(DIVEC3-ADD);
DIV_A2=(DIV_A).^2;
SQ_ADD=sum(DIV_A2(1,:));
load('LOSS.mat');%LOSS is 1x229
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;
load('END.mat'); %END is 1x229
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;
SQ_TOT= SQ_LOSS+SQ_END+SQ_ADD

  0 Comments

Sign in to comment.

2 Answers

Answer by 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

  6 Comments

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];
load('DUM.mat'); %4x229
%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));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Addition%%%%
load('housed.mat'); %1x229
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;
load('DishADD.mat'); %1x229
DIV_=ADD~=0;
DIVEC3=DIV_.*EC3m;
DIV_A=(DIVEC3-ADD);
DIV_A2=(DIV_A).^2;
SQ_ADD=sum(DIV_A2(1,:));
load('LOSS.mat'); %1x229
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);
load('END.mat');
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;
SQ_TOT=SQ_LOSS+SQ_END+SQ_ADD
end
Sorry now it work thank you very much to you all!!
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];
load('DATADISH12.mat');
DUM=DATADISH12(1:4,:);
%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));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Addition%%%%
housed=DATADISH12(5,:);
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;
ADD=DATADISH12(6,:);
DIV_=ADD~=0;
DIVEC3=DIV_.*EC3m;
DIV_A=(DIVEC3-ADD);
DIV_A2=(DIV_A).^2;
SQ_ADD=sum(DIV_A2(1,:));
LOSS=DATADISH12(7,:);
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);
END=DATADISH12(8,:);
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;
SQ_TOT=SQ_LOSS+SQ_END+SQ_ADD
end

Sign in to comment.


Answer by 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:
load('DATADISH12.mat');
DUM = DATADISH12(1:4,:);
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

  0 Comments

Sign in to comment.