Forecasting the number of deaths - for all T

Modules

Dead

Dead[DataT_,Ind_:0]:=Module[{},
danIni=First[DataT][[{1,2}]]; danEnd=Last[DataT][[{1,2}]];
delta=10;step=10;
Lpar=Gpar=0;
ticksx=Table[{i,Data[[i,{1,2}]]},{i,delta,Length[Data],step}];
If[Ind≠0,
novoz=ListPlot[DataT[[All,4]],Joined->True, Filling->{1->{1/2,LightRed}} ,Ticks->{ticksx,Automatic},GridLines->{tabx,Automatic},PlotLabel->"New cases from "<>ToString[danIni]<>" to "<>ToString[danEnd]];
umrli=ListPlot[DataT[[All,5]],Joined->True, Filling->{1->{1/2,LightGray}} ,Ticks->{ticksx,Automatic},GridLines->{tabx,Automatic},PlotLabel->"Deceased from "<>ToString[danIni]<>" to "<>ToString[danEnd]];
Print[GraphicsGrid[{{novoz,umrli}},ImageSize->600]]
 ];
delU=15;  (* Data delay for deceased and hospitalized *)
pod=Table[{DataT[[i,3]],DataT[[i,4]]},{i,Length[DataT]}] ;  (* Daily data of new cases *)
dead=Table[{DataT[[i,3]],DataT[[i,5]]},{i,delU,Length[DataT]}] ;  (* Daily data of deceased with delay *)
axesstyle={Thickness[.004],Thickness[.004]};
slpod=ListPlot[pod,PlotStyle->{Gray,PointSize[.015]},AxesStyle->axesstyle];
sldead=ListPlot[dead,PlotStyle->{Gray,PointSize[.015]},AxesStyle->axesstyle];
				(* Cumulative data of new cases *)
POD=pod;
Do[
POD[[i]]={pod[[i,1]],pod[[i,2]]+POD[[i-1,2]]}
,{i,2,Length[pod]}];
sltocN=ListPlot[POD,PlotRange->All,PlotStyle->{Gray,PointSize[.015]}];
			(* Logistic model-function of new cases:  DIRECT  *)
Print[Style["Application of the Logistic model function:",16,Red]];
f[t_,a_,b_,c_]=a/(1+b Exp[-c t]);
FCompile=Compile[{{a,_Real},{b,_Real},{c,_Real}},
Evaluate[ Sum[(f[POD[[i,1]],a,b,c] -POD[[i,2]])^2,{i,Length[POD]}]]];	
F[a_,b_,c_]:=Apply[FCompile,{a,b,c}];	
dom={{1,2Last[POD][[2]]},{0,5Last[POD][[2]]/Mean[Table[POD[[s,2]],{s,5}]]},{0,10}};
Tmin=Timing[min=DIRECT[F,dom,DMin->10.^(-3),MaxIter->100, Eps->10^(-4)]];
If[Ind≠0,Print[Style["Parameter estimation of the logistic model function:", Blue,16] ,"\nDIRECT:" ,min]];
a0=min[[1,1]];b0=min[[1,2]]; c0=min[[1,3]];
					(* Logistic model-function:  NonlinearModelFit *)
nlm=NonlinearModelFit[POD,f[t,a,b,c],{{a,a0},{b,b0},{c,c0}},t];aN1=a/.nlm[[1,2]]; bN1=b/.nlm[[1,2]]; cN1=c/.nlm[[1,2]];
I1={Log[bN1]/cN1,aN1/2};  
tB=Log[bN1/(2+Sqrt[3])]/cN1;
tC=Log[bN1/(2-Sqrt[3])]/cN1;
 B={tB,f[tB,aN1,bN1,cN1]};  
CC={tC,f[tC,aN1,bN1,cN1]};
Print["NonlinearModelFit Parameters: (SuperStar[a]) = ",aN1,"; (SuperStar[b]) = ",bN1,"; (SuperStar[c]) = ",cN1];
If[Ind≠0,
Print["I = ",Round[I1,-1],";  B = ",Round[B,-1],";  C = ",Round[CC,-1],";   Saturation level: ",Round[aN1,-1]];Print[Show[sltocN,ListPlot[{I1,B,CC},PlotStyle->{Blue,PointSize[.03],Opacity[.8]}],Plot[{aN1,nlm[t]},{t,0,2I1[[1]]},PlotStyle->Gray],PlotRange->{0,aN1+1},AxesStyle->axesstyle,Ticks->{ticksx,Automatic},GridLines->{tabx,Automatic},PlotLabel->"Logistic model function of new cases",ImageSize->Medium]];
	];
(* The function ϕ of deceased according to logistic model-function and prediction; 
Initial appr. by using DIRECT, and then NonlinearModelFit *)
phi[t_]:=aN1/(1+bN1 Exp[-cN1 t]); (* Logistic model-function of new cases *)
						(* DIRECT *)
FCompile=Compile[{{u,_Real},{al,_Real},{tau,_Real}},
Evaluate[ Sum[(u+al phi'[dead[[i,1]]-tau] -dead[[i,2]])^2,{i,Length[dead]}]]];	
F[u_,al_,tau_]:=Apply[FCompile,{u,al,tau}];	
dom={{0,1000},{0,1},{0,20}};
Tmin=Timing[min=DIRECT[F,dom,DMin->10.^(-4),MaxIter->100, Eps->10^(-4)]];
If[Ind≠0,Print[Style["\nParameter estimation of the model function ϕ:",16,Blue],"\nDIRECT:", min]];
u1=min[[1,1]];al1=min[[1,2]];  tau1=min[[1,3]];
								(* NonlinearModelFit *)
nlm=NonlinearModelFit[dead,u+al phi'[t-tau],{{u,u1},{al,al1},{tau,tau1}},t];
u1=u/.nlm[[1,2]];
al1=al/.nlm[[1,2]];alL=Round[100 al1,-.01];
tau1=tau/.nlm[[1,2]];
Print["NonlinearModelFit parameters: (SuperStar[u]) = ",u1,"; (SuperStar[α]) = ",al1,"; (SuperStar[τ]) = ",tau1];
If[Ind≠0,Print[Plot[u1+al1 phi'[t-tau1],{t,0,Length[Data]+10},Ticks->{ticksx,Automatic},GridLines->Automatic,PlotLabel->"Model function ϕ",ImageSize->300]]];
Rtau1=Round[tau1];tauL=Rtau1;
If[u1>0 && al1>0 &&Rtau1≥ 1,
If[Rtau1<22,
PP=Table[Round[ u1+al1 phi'[s-tau1],-1],{s,Length[pod]+1,Length[pod]+Rtau1}];
If[Ind≠0,slphid=Plot[u1+al1 phi'[t-tau1],{t,0,Last[dead][[1]]+tau1},PlotStyle->{Red}];
tickstau=Table[{tabx[[i]],DataT[[10i,{1,2}]]},{i,(Floor[Length[pod]])/10}];
(* tabx=Table[i,{i,First[dead][[1]]-5, T,10}];*)
taby=Table[i,{i,20,Round[ Max[dead[[All,2]]]],20}];
					(* BANDS *)
{bands90[ t_],bands95[ t_],bands99[ t_]}=Table[nlm["SinglePredictionBands",ConfidenceLevel->cl],{cl,{.90,.95,.99}}];Print[Show[ Plot[{nlm[t],bands90[t],bands95[t],bands99[t]},{ t,First[dead][[1]]-5, T},Filling->{2->{1},3->{2},4->{3}},PlotLabel->"Confidence bands",PlotLegends->{"model function ϕ", "90% confidence bands","95% confidence bands","99% confidence bands"}],ListPlot[dead],
GridLines->{tabx,taby},Ticks->{tickstau,Automatic}]];
	];
SL=Sum[PP[[s]],{s,Rtau1}]; (* Sum of (logistic) deceased *)
Print["According to logistic model function:","\nProportion of deaths among cases confirmed earlier: ",Round[100 al1,-.1],"%"];Print["Predicted number of deaths during ",Rtau1," days from ",danEnd," is ",SL];
If[T+Rtau1<Length[Data],
SAL=Sum[Data[[s,5]],{s,T+1,T+Rtau1}]; 
Lape=Round[Abs[100.(SL-SAL)/SL],-.1];
Print["Actual deaths = ",SAL,";  Absolute Percentage Error = ",Lape,"%"];
],
Lpar=13;Print["GOP (5) for logistic model function does not give a sensible solution (predicting for more than 3 weeks makes no sense): (OverHat[τ]) = ",Rtau1," >  21"]
],
Lpar=13;Print["GOP (5) for logistic model function does not give a sensible solution: \none or several parameters are not positive"]
];
			(* Gompertz model-function of new infections *)
Print[Style["\nApplication of Gompertz model function",Red,18]];
f[t_,a_,b_,c_]:=Exp[a-b Exp[-c t]];
							(* DIRECT *)
FCompile=Compile[{{a,_Real},{b,_Real},{c,_Real}},
Evaluate[ Sum[(f[POD[[i,1]],a,b,c] -POD[[i,2]])^2,{i,Length[POD]}]]];	
F[a_,b_,c_]:=Apply[FCompile,{a,b,c}];
a0=Log[1.5Last[POD][[2]]]//N;	
dom={{0,a0},{0,a0-Log[Mean[Table[POD[[s,2]],{s,5}]]]},{0,10}};
Tmin=Timing[min=DIRECT[F,dom,DMin->10.^(-3),MaxIter->100, Eps->10^(-4)]];
If[Ind≠0,Print[Style["Parameter estimation of the Gompertz model function: ",Blue,16],"\nDIRECT: ",min]];
						(* NonlinearModelFit *)
a0=min[[1,1]];b0=min[[1,2]]; c0=min[[1,3]];
nlm=NonlinearModelFit[POD,f[t,a,b,c],{{a,a0},{b,b0},{c,c0}},t];aN1=a/.nlm[[1,2]]; bN1=b/.nlm[[1,2]]; cN1=c/.nlm[[1,2]];
If[aN1>0 && bN1>0 && cN1>0,
I1={Log[bN1]/cN1 ,Exp[aN1-1]};    
tB=Log[bN1(3-Sqrt[5])/2]/cN1;
tC=Log[bN1(3+Sqrt[5])/2]/cN1;
 B={tB,f[tB,aN1,bN1,cN1]}; 
CC={tC,f[tC,aN1,bN1,cN1]};
slA=Plot[Exp[aN1],{t,0,tC+10}]; 
Print["NonlinearModelFit Parameters:  (SuperStar[a]) = ",aN1,"; (SuperStar[b]) = ",bN1,"; (SuperStar[c]) = ",cN1];
If[Ind≠0,
Print["I = ",Round[I1,-1],";  B = ",Round[B-1],";  C = ",Round[CC,-1],";   Saturation level: ",Round[Exp[aN1],-1]];Print[Show[sltocN,slA,ListPlot[{I1,B,CC},PlotStyle->{Blue,PointSize[.03],Opacity[.8]}],Plot[{aN1,nlm[t]},{t,0,tC+10},PlotStyle->Gray],PlotRange->{0,Exp[aN1]},AxesStyle->axesstyle,Ticks->{ticksx,Automatic},GridLines->{tabx,Automatic},PlotLabel->"Gompertz model function",ImageSize->Medium]];
	];
(* The function ϕ of deceased according to Gompertz model-function and prediction; 
Initial appr. first using DIRECT, and then using NonlinearModelFit*)
phi[t_]:=Exp[aN1-bN1 Exp[-cN1 t]];       (* Gompertz model-function of new cases *)
								(* DIRECT *)
FCompile=Compile[{{u,_Real},{al,_Real},{tau,_Real}},
Evaluate[ Sum[(u+al phi'[dead[[i,1]]-tau] -dead[[i,2]])^2,{i,Length[dead]}]]];	
F[u_,al_,tau_]:=Apply[FCompile,{u,al,tau}];	
dom={{0,1000},{0,2},{0,50}};
Tmin=Timing[min=DIRECT[F,dom,DMin->10.^(-4),MaxIter->15, Eps->10^(-4)]];
If[Ind≠0,Print[Style["\nParameter estimation of the model function ϕ:",Blue,16],"\nDIRECT: ",min]];
u1=min[[1,1]]; al1=min[[1,2]]; tau1=min[[1,3]];
								(* NonlinearModelFit *)
nlm=NonlinearModelFit[dead,u+al phi'[t-tau],{{u,u1},{al,al1},{tau,tau1}},t];
u1=u/.nlm[[1,2]];
al1=al/.nlm[[1,2]];alG=Round[100 al1,-.01];
tau1=tau/.nlm[[1,2]];
Gtau1=Round[tau1];tauG=Gtau1;
Print["NonlinearModelFit parameters: (SuperStar[u]) = ",u1,"; (SuperStar[α]) = ",al1,"; (SuperStar[τ]) = ",tau1];
If[u1>0 && al1>0 &&Gtau1≥1,
If[Gtau1<22,
PP=Table[Round[u1+ al1 phi'[s-tau1],-1],{s,Length[pod]+1,Length[pod]+Gtau1}];
If[Ind≠0,slphid=Plot[u1+al1 phi'[t-tau1],{t,0,Last[dead][[1]]+tau1},PlotStyle->{Red}];
{bands90[ t_],bands95[ t_],bands99[ t_]}=Table[nlm["SinglePredictionBands",ConfidenceLevel->cl],{cl,{.90,.95,.99}}];Print[Show[ Plot[{nlm[t],bands90[t],bands95[t],bands99[t]},{ t,First[dead][[1]]-5, T},Filling->{2->{1},3->{2},4->{3}},PlotLabel->"Confidence bands",PlotLegends->{"model function ϕ", "90% confidence bands","95% confidence bands","99% confidence bands"}],ListPlot[dead],ListPlot[dead], GridLines->{tabx,taby},Ticks->{tickstau,Automatic}]];
	];
SG=Total[PP];
Print["According to Gompertz model function:"];
Print["Proportion of deaths among cases confirmed earlier: ",Round[100 al1,-.1],"%"];Print["Predicted number of deaths during ",Gtau1," days from ",danEnd," is ",SG];
If[T+Gtau1<Length[Data],
SAG=Sum[Data[[s,5]],{s,T+1,T+Gtau1}]; 
Gape=Round[Abs[100.(SG-SAG)/SG],-.1];
Print["Actual deaths = ",SAG,";  Absolute Percentage Error = ",Gape,"%"];
    ],
Gpar=15;Print["GOP (5) for Gompertz model function does not give a sensible solution (predicting for more than 3 weeks makes no sense): (OverHat[τ]) = ",Gtau1," >  21"]
],
Gpar=15;Print["GOP (5) for Gompertz model function does not give a sensible solution:\none or several parameters are not positive "]
],
Gpar=15;
Print["The optimal parameters of the Gompertz model function were not found\nThe calculation process for Gompertz model function cannot be continue"]
];
]

DIRECT Modul

Options[DIRECT]={Eps->0,MaxIter->20,DMin->10.^(-2)};
DIRECT[f_,dom_,OptionsPattern[]]:=Module[{n,ff,c,d,dmax,fc,dmaxrazni,pot0,pot1,pot2,temp,ind,novic,novid,novifc,red,TMin,iter=0,dtest,anim={}},
dtest=N[OptionValue[DMin]/Max[(dom[[All,2]]-dom[[All,1]])]];
DIRECT::maxiter="The maximal number of iterations (`1`) has been reached.";
n=Length[dom];
(* Scaling the domain to the unit hypercub *)
ff=Function[N[f@@({##}(dom[[All,2]]-dom[[All,1]])+dom[[All,1]]),50]];
(* Initial center and diameters *)
c={Table[1/2,{i,n}]};
d={Table[1,{i,n}]};
dmax={Max[d[[1]]]};
dmaxrazni=Union[dmax];
fc={ff@@c[[1]]};
TMin={c[[1]],fc[[1]]};
(* MAIN LOOP *) 
Monitor[
While[(dmaxrazni[[1]]>dtest)&&(iter<OptionValue[MaxIter]),
iter++;
(* STEP 1 :searching for potentially optimal hyper-rectangles *)
dmaxrazni=Union[dmax];
(* All data grouped according to d *)
pot0=Table[Select[Range[Length[fc]],dmax[[#]]==dmaxrazni[[i]]&],{i,Length[dmaxrazni]}];
(* Data with minimal f grouped according to d *)
pot1=Table[pot0[[i,Position[fc[[pot0[[i]]]],Min[fc[[pot0[[i]]]]]][[All,1]]]],{i,Length[dmaxrazni]}];
(* Checking other criteria *)
pot2={pot1[[-1]]};
Do[temp=Min[Table[(fc[[pot1[[i,1]]]]-fc[[pot2[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot2[[j,1]]]]),{j,Length[pot2]}]];
If[temp>OptionValue[Eps]&&(i==1||Max[Table[(fc[[pot1[[i,1]]]]-fc[[pot1[[j,1]]]])/(dmax[[pot1[[i,1]]]]-dmax[[pot1[[j,1]]]]),{j,i-1}]]<temp)
,AppendTo[pot2,pot1[[i]]]];
,{i,Length[dmaxrazni]-1,1,-1}];
pot2=Flatten[pot2,1];
(* STEP 2: Dividing potentially optimal hyper-rectangles *)
Do[
(* Dividing directions *)
ind=Flatten[Position[d[[pot]],Max[d[[pot]]]]];
novic=Table[c[[pot]],{2*Length[ind]}];
novid=Table[d[[pot]],{2*Length[ind]}];
d[[pot,ind]]/=3;
dmax[[pot]]/=3;
Do[
novic[[2i-1,ind[[i]]]]-=d[[pot,ind[[i]]]];
novic[[2i,ind[[i]]]]+=d[[pot,ind[[i]]]];
,{i,Length[ind]}];
novifc=Table[ff@@novic[[i]],{i,Length[novic]}];
red=Ordering[Table[Min[{novifc[[2i-1]],novifc[[2i]]}],{i,Length[ind]}]];
Do[novid[[2red[[j]]-1,ind[[red[[i]]]]]]=novid[[2red[[j]],ind[[red[[i]]]]]]/=3,{i,1,Length[red]},{j,i,Length[red]}];
c=Join[c,novic];
d=Join[d,novid];
dmax=Join[dmax,Map[Max,novid]];
fc=Join[fc,novifc];
If[novifc[[2red[[1]]-1]]<TMin[[2]],TMin=N[{novic[[2red[[1]]-1]],novifc[[2red[[1]]-1]]},10]];
If[novifc[[2red[[1]]]]<TMin[[2]],TMin=N[{novic[[2red[[1]]]],novifc[[2red[[1]]]]},10]];
,{pot,pot2}];
(*Print[iter,": ",{TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]],TMin[[2]]}];*)
If[iter==OptionValue[MaxIter],Message[DIRECT::maxiter,OptionValue[MaxIter]];];
(* Figure for the 2D case *)
(*If[n==2,
AppendTo[anim,
Graphics[{{Green,Rectangle[]},Table[{Blue,Rectangle[c[[i]]-d[[i]]/2+d[[i]]/10,c[[i]]+d[[i]]/2-d[[i]]/10]},{i,Length[d]}],Point[c]}]
];
];*)
];
,{iter,Length[pot2]}];
(* Back to the original domain *)
TMin[[1]]=TMin[[1]](dom[[All,2]]-dom[[All,1]])+dom[[All,1]];
(* Printing figures for the 2D case *)
(*If[n==2,Print[ListAnimate[anim,AnimationRunning->False]]];*)
Return[TMin];
]

Data input

SetDirectory[NotebookDirectory[]]; 
Data=ReadList["Croatia.txt"];
(*Print["|Data| = ",Length[Data]]*)
(*Print["First datum: ",First[Data]]*)
(*Print["Last datum: ",Last[Data]]*)

(* To change the initial datum set k>1 (default is k=1) *)
k=1;
Data=Table[Data[[i]],{i,k,Length[Data]}];
Print["|Data| = ",Length[Data]]
Print["First datum: ",First[Data]]
Print["Last datum: ",Last[Data]]

(*      Graphical representation of all new cases             *)
dayIni=First[Data][[{1,2}]]; dayEnd=Last[Data][[{1,2}]];
ticks=Table[i,{i,0,Length[Data],10}];
tabx=Table[10i,{i,150}]; fig1=ListPlot[Data[[All,4]],Joined->True, Filling->{1->{1/2,LightRed}} ,Ticks->{ticks,Automatic},GridLines->{ticks,Automatic},PlotLabel->"New cases from "<>ToString[dayIni]<>" to "<>ToString[dayEnd],ImageSize->Medium];
(* 7-day moving average of new cases *)
DataAv=Table[N[Sum[Data[[All,4]][[n]],{n,n-3,n+3}]/7],{n,0,Length[Data]-3}];
fig2=ListPlot[DataAv,Joined->True,Ticks->{ticks,Automatic},GridLines->{ticks,Automatic},PlotLabel->"7-day moving average of new cases from "<>ToString[dayIni+{4,0}]<>" to "<>ToString[dayEnd-{2,0}],ImageSize->Medium];
Print[GraphicsGrid[{{fig1,fig2}},ImageSize->800]]

Exe with details

(* INPUT:25 ≤  T ≤ Length[Data], the point which divides training and testing data in Data *)
T=55;
If[25≤ T≤  Length[Data],
DataT=Table[Data[[s]],{s,T}];
Dead[DataT,1],
Print[Style["T must be between 25 and ",Red,12],Length[Data]]];

Exe without details

(* INPUT:25 ≤  T ≤ Length[Data], the point which divides training and testing data in Data *)
T=65;
If[25≤ T≤  Length[Data],
DataT=Table[Data[[s]],{s,T}];
Dead[DataT,0],
Print[Style["T must be between 25 and ",Red,12],Length[Data]]];

All training data sets

(* Input: jIni ≥ 25 and jEnd ≤ Length[Data] *)
jIni=35; jEnd=112;DEAD={};
If[jIni≥  25,
If[jEnd≤   Length[Data]&& jIni<jEnd,
Do[
Print[Style["\nT = ",Blue,18],Style[j,Blue,18]];
T=j;
DataT=Table[Data[[s]],{s,j}];
Dead[DataT,0];
If[Lpar== 13,
tsal=tsl=Lape=tauL=alL=-10,
If[j+Rtau1≥ Length[Data],
tsal=tsl=Lape=-10,
tsal=Total[SAL]; tsl=Total[SL]
  ]];
If[Gpar== 15,
tsag=tsg=Gape=tauG=alG=-10,
If[j+Gtau1≥ Length[Data],
tsag=tsg=Gape=-10,
tsag=Total[SAG]; tsg=Total[SG]
  ]];
DEAD=Append[DEAD,{tsal,tsl,tsag,tsg,Lape,Gape,tauL,tauG,alL,alG}];
,{j,jIni,jEnd}];
						(* Graphic representations *)
tab=Table[i,{i,100}];
ticksX=Table[{s,jIni+s-1},{s,1,jEnd-jIni+1,5}];Print[ListPlot[{DEAD[[All,1]],DEAD[[All,2]]},PlotStyle->{{Green,PointSize[.02],Opacity[.5]},{Blue,PointSize[.02],Opacity[.5]}},PlotRange->{0,1.1Max[DEAD[[All,1]]]},Ticks->{ticksX,Automatic},GridLines->{tab,Automatic},AxesLabel->{"T","Deceased"},PlotLabel->"Deceased (according to Logistic)"<>ToString[ {jIni,jEnd}],PlotLegends->{Actual,Logistic},ImageSize->300]];
Print[ListPlot[{DEAD[[All,3]],DEAD[[All,4]]},PlotStyle->{{Green,PointSize[.02],Opacity[.5]},{Red,PointSize[.02],Opacity[.5]}},Ticks->{ticksX,Automatic},PlotRange->{0,1.1Max[DEAD[[All,3]]]},GridLines->{tab,Automatic},AxesLabel->{"T","Deceased"},PlotLabel->"Deceased (according to Gompertz)"<>ToString[ {jIni,jEnd}],PlotLegends->{Actual,Gompertz},ImageSize->300]];
Print[Show[ListPlot[{DEAD[[All,5]],DEAD[[All,6]]},PlotStyle->{{Blue,PointSize[.02],Opacity[.5]},{Red,PointSize[.02],Opacity[.5]}},Ticks->{ticksX,Automatic},PlotRange->{0,1.5Max[DEAD[[All,5]]]},GridLines->{tab,Automatic},AxesLabel->{"T","%"},PlotLabel->"Error "<>ToString[ {jIni,jEnd}],PlotLegends->{Logistic,Gompertz},ImageSize->300],ListPlot[DEAD[[All,5]],PlotStyle->{Blue}],ListPlot[DEAD[[All,6]],PlotStyle->{Red}]]];
Print[ListPlot[{DEAD[[All,7]],DEAD[[All,8]]},PlotStyle->{{Blue,PointSize[.02],Opacity[.5]},{Red,PointSize[.02],Opacity[.5]}},Ticks->{ticksX,Automatic},PlotRange->{0,1.1Max[DEAD[[All,7]]]},GridLines->{tab,Automatic},AxesLabel->{"T","Days"},PlotLabel->"Delay "<>ToString[ {jIni,jEnd}],PlotLegends->{Logistic,Gompertz},ImageSize->300]];
Print[ListPlot[{DEAD[[All,9]],DEAD[[All,10]]},PlotStyle->{{Blue,PointSize[.02],Opacity[.5]},{Red,PointSize[.02],Opacity[.5]}},Ticks->{ticksX,Automatic},PlotRange->{0,1.1Max[DEAD[[All,9]]]},GridLines->{tab,Automatic},AxesLabel->{"T","%"},PlotLabel->"Proportion of deceased "<>ToString[ {jIni,jEnd}],PlotLegends->{Logistic,Gompertz},ImageSize->300]],
Print[Style["Upper bound of jEnd must be ≤ ",Red,12],Length[Data],Style[" and jIni < jEnd",Red,12]]],
Print[Style["Lower bound of jIni must be ≥ 25",Red,12]]];
Print["Gaps in displayed figures occur for those T's for which the GOP (9) does not give a correct solution."]