Add files via upload

fit_SEIQRDP.m

@@ -103,14 +103,25 @@
     kappaMin = guess(8:10)*0.95; % Constrain the guess if no R available
     lambdaMax = [1 1 100]; % bound the guess around the initial fit
     lambdaMin = [0 0 0]; % bound the guess around the initial fit   
+
+    if kappaMax(3)<1e-1
+        kappaMax(3) = 100;
+        kappaMin(3) = 0;
+     end
+
 else
-    kappaMax = guess(8:10)*3; % bound the guess around the initial fit
-    kappaMin = guess(8:10)/3; % bound the guess around the initial fit
-    lambdaMax = guess(5:7)*3; % bound the guess around the initial fit
-    lambdaMin = guess(5:7)/3; % bound the guess around the initial fit
+    kappaMax = guess(8:10)*2.0; % bound the guess around the initial fit
+    kappaMin = guess(8:10)/2.0; % bound the guess around the initial fit
+    lambdaMax = guess(5:7)*2.0; % bound the guess around the initial fit
+    lambdaMin = guess(5:7)/2.0; % bound the guess around the initial fit
+
+     if kappaMax(3)<1e-1
+        kappaMax(3) = 20;
+        kappaMin(3) = 0;
+     end
 
-    if lambdaMax(3)<1e-2
-        lambdaMax(3) = 100;
+    if lambdaMax(3)<1e-1
+        lambdaMax(3) = 20;
         lambdaMin(3) = 0;
     end
 end
@@ -273,14 +284,21 @@
                 % A death rate larger than 3 is abnormally high. It is not
                 % used for the fitting.                 
                 rate(abs(rate)>3)=nan;
+
+                % Remove death rate = 0 if the majority number of death is not
+                % zero
+                if numel(find(rate==0))/numel(rate) <0.5
+                    rate(abs(rate)==0)=nan;
+                end
+
                 [coeff1,r1] = lsqcurvefit(@(para,t) myFun1(para,t),...
                     guess(8:10),x(~isnan(rate)),rate(~isnan(rate)),[0 0 0],[1 1 100],opt);
                 [coeff2,r2] = lsqcurvefit(@(para,t) myFun2(para,t),...
                     guess(8:10),x(~isnan(rate)),rate(~isnan(rate)),[0 0 0],[1 1 100],opt);
                 [coeff3,r3] = lsqcurvefit(@(para,t) myFun3(para,t),...
                     guess(8:10),x(~isnan(rate)),rate(~isnan(rate)),[0 0 0],[1 1 100],opt);
 %                 
-%                  figure;plot(x,rate,x,myFun1(coeff1,x),'r',x,myFun2(coeff2,x),'g',x,myFun3(coeff3,x),'b')
+%                     figure;plot(x,rate,x,myFun1(coeff1,x),'r',x,myFun2(coeff2,x),'g',x,myFun3(coeff3,x),'b')
 
                 minR = min([r1,r2,r3]);
                 if r1==minR

getMultipleWaves.m

@@ -1,4 +1,4 @@
-function [Q,R,D,T] = getMultipleWaves(guess,Npop,time,Confirmed,Recovered,Deaths,tStart1,tStart2,tEnd)
+function [Q,R,D,T] = getMultipleWaves(guess,Npop,time,Confirmed,Recovered,Deaths,tStart1,tStart2,tEnd,Q0,E0,I0)
 % [Q,R,D,T] =
 % getMultipleWaves(guess,Npop,time,Confirmed,Recovered,Deaths,tStart1,tStart2,tEnd)
 % simulate the number of recovered, deaths and active cases for the
@@ -14,6 +14,9 @@
 %   tStart1: datetime [1x1]: Initial time for the first wave
 %   tStart2: datetime [1x1]: Initial time for the second wave
 %   tEnd: datetime [1x1]: Final time for the simulation
+%   Q0: datetime [1x1]: Initial number of quarantined cases
+%   E0: datetime [1x1]: Initial number of exposed cases
+%   I0: datetime [1x1]: Initial number of infectious cases
 % 
 % Outputs
 %   Q: double [1xN1]: Time histories of the quarantined/active cases 
@@ -41,12 +44,8 @@
 
 %% Simulate first wave
 % Initial conditions
-Q0 = Active(indT1(1));
-E0 = 0.25*Q0; % Initial number of exposed cases. Unknown but unlikely to be zero.
-I0 = 5*E0; % Initial number of infectious cases. Unknown but unlikely to be zero.
 R0 = Recovered(indT1(1)); 
 D0 = Deaths(indT1(1)); 
-
 [E1,I1,Q1,R1,D1,T1] = computeWave(Active(indT1),Recovered(indT1),...
     Deaths(indT1),E0,I0,Q0,R0,D0,time(indT1),tStart1,tStart2,guess);