how to estimate paramters.

조회 수: 13 (최근 30일)
Chinwendu Madubueze
Chinwendu Madubueze 2020년 5월 28일
답변: Chinwendu Madubueze 2020년 6월 2일
I am try to estimate parameters using the following matlab code. it keep giving me this message below.
Error using ==> lsqncommon at 101
LSQCURVEFIT cannot continue because user supplied objective function failed with the following error:
Undefined function or variable "x".
Error in ==> lsqcurvefit at 182
[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...
Error in ==> lassafitw at 75
[x,Rsdnrm] = lsqcurvefit(@Model_Inc,x0,xdata,ydata,LB,UB,optimset);
function lassafitw(do_estimation)
warning off;
P_Data(:,1)=[1:80];
P_Data(:,2)=[ 0 16 117 252 305 386 454 525 574 646 717 783 848 910 1026 1261 1361 1424 1620 1673 1813 1940 2021 2304 2437 2950 3252 3410 3706 ...
3896 5235 5338 4759 4862 5368 5586 6073 6190 6599 7109 7312 7897 8356 9004 9446 9780 10124 10934 11103 11466 11751 11841 11974 12022 12244 ...
12313 12440 12593 127666 12735 12827 12884 12962 13012 13093 13150 13241 13402 13518 13586 13638 13701 13785 13894 13928 13945 13992 ...
14001 14052 14122]./3800;
piH = 20;
muH = 0.02;
d1=0.2;
piR=5;
eR=0.01;
rhoR=0.01;
muR=0.02;
d2=0.15;
rho3=0.05;
rho6=0.6;
rho2=0.6;
rho1=0.02;
rho7=0.2;
rho5=0.03;
omega=0.6;
rho4=0.5;
gamma=0.03;
beta1=0.2;
beta2=0.1;
beta3=0.08;
beta4=0.2;
beta5=0.02;
Sh0=100000; Eh0=2000; Eth0=0; Enh0=0; Ih0=0; Iqh0=0;Dh0=0; Sr0=500;Er0=0;Ir0=0;
INITIAL=[Sh0, Eh0, Eth0, Enh0, Ih0,Iqh0,Dh0, Sr0,Er0,Ir0];
Istart=1; %month to start the model simulation
Iend=Istart + 81;
OPTIONS=odeset('AbsTol',0.0001,'RelTol',0.0001,'MaxStep',1/12);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimate parameters
%by minimizing the sum of squares
%when fitting modeled to real prevalence data
do_estimation=1;
if(do_estimation)
xdata=P_Data(:,1)';
ydata=P_Data(:,2)';
x0(1,1)= 20; %Lambda
x0(1,2)=0.02; %mu
x0(1,3)= 0.2; %beta
x0(1,4)= 5; %0.1; %sigma
x0(1,5)=0.01; %0.005; %gamma
x0(1,6)=0.01; %0.005; %p
x0(1,7)= 0.02; %0.00001; %theta
x0(1,8)= 0.15; %0.0008; %delta1
x0(1,9)= 0.05; %0.004; %delta2
x0(1,10)= 0.6; %c
x0(1,11)= 0.6; %omega
x0(1,12)= 0.02; %alpha1
x0(1,13)= 0.2; %alpha2
x0(1,14)= 0.03; %alpha3
x0(1,15)= 0.6; %alpha4
x0(1,16)= 0.5; %alpha5
x0(1,17)= 0.03;
x0(1,18)= 0.2;
x0(1,19)= 0.1;
x0(1,20)= 0.08;
x0(1,21)= 0.2;
x0(1,22)= 0.02;
LB=[10 0 0.0001 5 0.0 0.000 0.0060 0.0060 0.0060 1 0.008 0.003 0 0 0 0 0 0 0 0 0 0 ];
UB=[50 1 1.5 20 0.9 0.9 0.9 0.9 0.9 10 1 1 1 1 1 1 1 1 1 1 1 5 5 ];
[x,Rsdnrm] = lsqcurvefit(@Model_Inc,x0,xdata,ydata,LB,UB,optimset);
'estimated parameters'
piH=x(1);
muH=x(2);
d1=x(3);
piR=x(4);
eR=x(5);
rhoR=x(6);
muR=x(7);
d2=x(8);
rho3=x(9);
rho6=x(10);
rho2=x(11);
rho1=x(12);
rho7=x(13);
rho5=x(14);
omega=x(15);
rho4=x(16);
gamma=x(17);
beta1=x(18);
beta2=x(19);
beta3=x(20);
beta4=x(21);
beta5=x(22);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t y] = ode45(@lassafitw,[0:1/12:(Iend-Istart)], INITIAL);
Sh=y(:,1);
Eh=y(:,2);
Eth=y(:,3);
Enh=y(:,4);
Ih=y(:,5);
Iqh=y(:,6);
Dh=y(:,7);
Sr=y(:,8);
Er=y(:,9);
Ir=y(:,10);
incidence= beta1*y(:,5) + beta2*y(:,6);
%incidence= y(:,4)*gamma;
close all;
figure(1);
hold on
h_l=plot(Istart+t,incidence,'r-');
set(h_l,'linewidth',2);
h_l=plot(P_Data(:,1),P_Data(:,2),'bo','Markersize',6);
set(h_l,'linewidth',2);
axis([Istart Iend 0 9]);
%ylim([0 45.00]) % The line added to change the y limit
%title('Ebola Prevalence(%)');
xlabel('Days','fontsize',15)
ylabel('Ebola Incidence','fontsize', 15)
hold off
function [ydot]=lassafitw(t,y)
Sh=y(1);
Eh=y(2);
Eth=y(3);
Enh=y(4);
Ih=y(5);
Iqh=y(6);
Dh=y(7);
Sr=y(8);
Er=y(9);
Ir=y(10);
ydot(1)=piH - beta4*y(9)*y(1)/(y(8)+y(9)+y(10)) -(beta1*y(5)+beta2*y(6)+beta3*y(7))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7))-muH*y(1)+rho3*y(3)+rho6*y(7);
ydot(2)=beta4*y(9)*y(1)/(y(8)+y(9)+y(10))+(beta1*y(5)+beta2*y(6)+beta3*y(7))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)) -((1-omega)*rho2+ omega*rho1+muH)*y(2);
ydot(3)=(1-omega)*rho2*y(2)-(rho3+rho4+muH)*y(3);
ydot(4)=omega*rho1*y(2)-(rho5+muH)*y(4);
ydot(5)=rho5*y(4)-(rho7+d1+muH)*y(5);
ydot(6)=rho7*y(5)+rho4*y(3)-(rho6+d2+muH)*y(6);
ydot(7)=d1*y(5)+d2*y(6) - gamma*x(7);
ydot(8)=piR - beta5*y(9)*y(8)/NR - (muR+rhoR)*y(8);
ydot(9)=beta5*y(9)*y(8)/NR -(muR+rhoR+eR)*y(9);
ydot(10)=eR*y(9)-(muR+rhoR)*y(10);
ydot = ydot';
end % function
function Inc=Model_Inc(x0,xdata)
% Inc=0; %initialization of this not to have an empty array
piH=x0(1);
muH=x0(2);
d1=x0(3);
piR=x0(4);
eR=x0(5);
rhoR=x0(6);
muR=x0(7);
d2=x0(8);
rho3=x0(9);
rho6=x0(10);
rho2=x0(11);
rho1=x0(12);
rho7=x0(13);
rho5=x0(14);
omega=x0(15);
rho4=x0(16);
gamma=x0(17);
beta1=x0(18);
beta2=x0(19);
beta3=x0(20);
beta4=x0(21);
beta5=x0(22);
%%% Initial values
Sh0=100000; Eh0=2000; Eth0=0; Enh0=0; Ih0=0; Iqh0=0;Dh0=0; Sr0=500;Er0=0;Ir0=0;
INITIAL=[Sh0, Eh0, Eth0, Enh0, Ih0,Iqh0,Dh0,Sr0,Er0,Ir0];
OPTIONS=odeset('AbsTol',0.0001,'RelTol',0.0001,'MaxStep',1/12);
tdur=50;
[t y] = ode45(@lassafitw, [0:1/12:Iend-Istart], INITIAL);
Sh=y(:,1);
Eh=y(:,2);
Eth=y(:,3);
Enh=y(:,4);
Ih=y(:,5);
Iqh=y(:,6);
Dh=y(:,7);
Sr=y(:,8);
Er=y(:,9);
Ir=y(:,10);
incidence= beta1*y(:,5) + beta2*y(:,6)
%incidence= y(:,4)*gamma;
for i=1:length(xdata)
ind=find(xdata(i)-Istart==t);
Sh=y(ind,1);
Eh=y(ind,2);
Eth=y(ind,3);
Enh=y(ind,4);
Ih=y(ind,5);
Iqh=y(ind,6);
Dh=y(ind,7);
Sr=y(ind,8);
Er=y(ind,9);
Ir=y(ind,10);
Inc(i)= beta1*Ih + beta2*Iqh; %incidence
%Inc(i)= Is*gamma; %incidence
end
end % end of Model_prev
Q1= gamma + mu;
Q2= mu + theta + delta1;
Q3 = delta2 + sigma + mu;
mstar= Lambda*gamma1 / omega* mu;
R_M = beta*c*gamma*(rho + eta*sigma)*(1-mstar)*(p*theta + Q2*(1-p))/(rho*Q1*Q2*Q3)
end %end of ebola_fit
  댓글 수: 4
the cyclist
the cyclist 2020년 6월 2일
I formatted the code for you.
Chinwendu Madubueze
Chinwendu Madubueze 2020년 6월 2일
okay. Thank you very much. I will run it again.

댓글을 달려면 로그인하십시오.

채택된 답변

Chinwendu Madubueze
Chinwendu Madubueze 2020년 6월 2일
I need help concerning parameter estimation and curve fitting using data. I am trying to estimate 22 paramters of the model using the the following matlab code. I keep geting this error in my work. Please I need help to correct my code. Thank you.
function lassafitw1(do_estimation)
warning off;
P_Data(:,1)=[1:80];
P_Data(:,2)=[ 0 16 117 252 305 386 454 525 574 646 717 783 848 910 1026 1261 1361 1424 1620 1673 1813 1940 2021 2304 2437 2950 3252 3410 3706 ...
3896 5235 5338 4759 4862 5368 5586 6073 6190 6599 7109 7312 7897 8356 9004 9446 9780 10124 10934 11103 11466 11751 11841 11974 12022 12244 ...
12313 12440 12593 127666 12735 12827 12884 12962 13012 13093 13150 13241 13402 13518 13586 13638 13701 13785 13894 13928 13945 13992 ...
14001 14052 14122]./3800;
piH = 0.20;
muH = 0.02;
d1=0.2;
piR=0.05;
eR=0.01;
rhoR=0.01;
muR=0.02;
d2=0.15;
rho3=0.05;
rho6=0.6;
rho2=0.6;
rho1=0.02;
rho7=0.2;
rho5=0.03;
omega=0.6;
rho4=0.5;
gamma=0.03;
beta1=0.2;
beta2=0.1;
beta3=0.08;
beta4=0.2;
beta5=0.02;
Sh0=100000; Eh0=2000; Eth0=0; Enh0=0; Ih0=0; Iqh0=0;Dh0=0; Sr0=500;Er0=0;Ir0=0;
INITIAL=[Sh0, Eh0, Eth0, Enh0, Ih0,Iqh0,Dh0, Sr0,Er0,Ir0];
Istart=1; %month to start the model simulation
Iend=Istart + 81;
OPTIONS=odeset('AbsTol',0.001,'RelTol',0.001,'MaxStep',1/12);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Estimate parameters
%by minimizing the sum of squares
%when fitting modeled to real prevalence data
do_estimation=1;
if(do_estimation)
xdata=P_Data(:,1)';
ydata=P_Data(:,2)';
x0(1,1)= 0.20; %Lambda
x0(1,2)=0.02; %mu
x0(1,3)= 0.2; %beta
x0(1,4)= 0.05; %0.1; %sigma
x0(1,5)=0.01; %0.005; %gamma
x0(1,6)=0.01; %0.005; %p
x0(1,7)= 0.02; %0.00001; %theta
x0(1,8)= 0.15; %0.0008; %delta1
x0(1,9)= 0.05; %0.004; %delta2
x0(1,10)= 0.6; %c
x0(1,11)= 0.6; %omega
x0(1,12)= 0.02; %alpha1
x0(1,13)= 0.2; %alpha2
x0(1,14)= 0.03; %alpha3
x0(1,15)= 0.6; %alpha4
x0(1,16)= 0.5; %alpha5
x0(1,17)= 0.03;
x0(1,18)= 0.2;
x0(1,19)= 0.1;
x0(1,20)= 0.08;
x0(1,21)= 0.2;
x0(1,22)= 0.02;
LB=[10 0 0.0001 5 0.0 0.000 0.0060 0.0060 0.0060 1 0.008 0.003 0 0 0 0 0 0 0 0 0 0 ];
UB=[50 1 1.5 20 0.9 0.9 0.9 0.9 0.9 10 1 1 1 1 1 1 1 1 1 1 1 5 5 ];
[x,Rsdnrm]=lsqcurvefit(@Model_Inc,x0,xdata,ydata,LB,UB,optimset);
'estimated parameters'
piH=x(1)
muH=x(2)
d1=x(3)
piR=x(4)
eR=x(5)
rhoR=x(6)
muR=x(7)
d2=x(8)
rho3=x(9)
rho6=x(10)
rho2=x(11)
rho1=x(12)
rho7=x(13)
rho5=x(14)
omega=x(15)
rho4=x(16)
gamma=x(17)
beta1=x(18)
beta2=x(19)
beta3=x(20)
beta4=x(21)
beta5=x(22)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t y] = ode45(@lassafitw1,[0:1/12:(Iend-Istart)], INITIAL);
Sh=y(:,1);
Eh=y(:,2);
Eth=y(:,3);
Enh=y(:,4);
Ih=y(:,5);
Iqh=y(:,6);
Dh=y(:,7);
Sr=y(:,8);
Er=y(:,9);
Ir=y(:,10);
incidence= beta1*y(:,5) + beta2*y(:,6);
%incidence= y(:,4)*gamma;
close all;
figure(1);
hold on
h_l=plot(Istart+t,incidence,'r-');
set(h_l,'linewidth',2);
h_l=plot(P_Data(:,1),P_Data(:,2),'bo','Markersize',6);
set(h_l,'linewidth',2);
axis([Istart Iend 0 9]);
%ylim([0 45.00]) % The line added to change the y limit
%title('Ebola Prevalence(%)');
xlabel('Days','fontsize',15)
ylabel('Ebola Incidence','fontsize', 15)
hold off
function [ydot]=lassafitw1(t,y)
Sh=y(1);
Eh=y(2);
Eth=y(3);
Enh=y(4);
Ih=y(5);
Iqh=y(6);
Dh=y(7);
Sr=y(8);
Er=y(9);
Ir=y(10);
ydot(1)=piH - beta4*y(9)*y(1)/(y(8)+y(9)+y(10)) -(beta1*y(5)+beta2*y(6)+beta3*y(7))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7))-muH*y(1)+rho3*y(3)+rho6*y(7);
ydot(2)=beta4*y(9)*y(1)/(y(8)+y(9)+y(10))+(beta1*y(5)+beta2*y(6)+beta3*y(7))*y(1)/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)) -((1-omega)*rho2+ omega*rho1+muH)*y(2);
ydot(3)=(1-omega)*rho2*y(2)-(rho3+rho4+muH)*y(3);
ydot(4)=omega*rho1*y(2)-(rho5+muH)*y(4);
ydot(5)=rho5*y(4)-(rho7+d1+muH)*y(5);
ydot(6)=rho7*y(5)+rho4*y(3)-(rho6+d2+muH)*y(6);
ydot(7)=d1*y(5)+d2*y(6) - gamma*x(7);
ydot(8)=piR - beta5*y(9)*y(8)/(y(8)+y(9)+y(10)) - (muR+rhoR)*y(8);
ydot(9)=beta5*y(9)*y(8)/(y(8)+y(9)+y(10)) -(muR+rhoR+eR)*y(9);
ydot(10)=eR*y(9)-(muR+rhoR)*y(10);
ydot = ydot';
end % function
function Inc=Model_Inc(x0,xdata)
% Inc=0; %initialization of this not to have an empty array
piH=x0(1);
muH=x0(2);
d1=x0(3);
piR=x0(4);
eR=x0(5);
rhoR=x0(6);
muR=x0(7);
d2=x0(8);
rho3=x0(9);
rho6=x0(10);
rho2=x0(11);
rho1=x0(12);
rho7=x0(13);
rho5=x0(14);
omega=x0(15);
rho4=x0(16);
gamma=x0(17);
beta1=x0(18);
beta2=x0(19);
beta3=x0(20);
beta4=x0(21);
beta5=x0(22);
%%% Initial values
Sh0=100000; Eh0=2000; Eth0=0; Enh0=0; Ih0=0; Iqh0=0;Dh0=0; Sr0=500;Er0=0;Ir0=0;
INITIAL=[Sh0, Eh0, Eth0, Enh0, Ih0,Iqh0,Dh0,Sr0,Er0,Ir0];
% OPTIONS=odeset('AbsTol',0.0001,'RelTol',0.0001,'MaxStep',1/12);
% tdur=50;
[t y] = ode45(@lassafitw1, [0:1/12:Iend-Istart], INITIAL);
Sh=y(:,1);
Eh=y(:,2);
Eth=y(:,3);
Enh=y(:,4);
Ih=y(:,5);
Iqh=y(:,6);
Dh=y(:,7);
Sr=y(:,8);
Er=y(:,9);
Ir=y(:,10);
incidence= beta1*y(:,5) + beta2*y(:,6)
%incidence= y(:,4)*gamma;
for i=1:length(xdata)
ind=find(xdata(i)-Istart==t);
Sh=y(ind,1);
Eh=y(ind,2);
Eth=y(ind,3);
Enh=y(ind,4);
Ih=y(ind,5);
Iqh=y(ind,6);
Dh=y(ind,7);
Sr=y(ind,8);
Er=y(ind,9);
Ir=y(ind,10);
Inc(i)= beta1*Ih + beta2*Iqh; %incidence
%Inc(i)= Is*gamma; %incidence
end
end % end of Model_prev
% Q1= gamma + mu;
% Q2= mu + theta + delta1;
% Q3 = delta2 + sigma + mu;
% mstar= Lambda*gamma1 / omega* mu;
% R_M = beta*c*gamma*(rho + eta*sigma)*(1-mstar)*(p*theta + Q2*(1-p))/(rho*Q1*Q2*Q3)
end %end of ebola_fit

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Descriptive Statistics에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by