# A Bimodal Lognormal Distribution Model for the Prediction of COVID-19 Deaths
## 3 Methods
### Numerical Implementation  
A nonlinear least-squares procedure is used to fit the available data and determine the model parameters. The fitting procedure is implemented into a MATLAB script called _COVID19_BiLog_World.m_, available in the Supplementary Materials.  
The OWID COVID-19 dataset is downloaded as a CSV file named _ovid-covid-data.csv_. For each location in the dataset, the daily number of deaths, \(d(t)\), and total number of deaths, \(D(t)\), are extracted. Here, the time variable, \(t\), is expressed in days. It runs from 1 (first day of dataset: 31 December 2019) to \(t_{\text{max}}=367\) (last day of simulations: 31 December 2020). The user can choose the days corresponding to the end of the fitting range, \(t_{\text{fit}}\), and the end of the extrapolation interval, \(t_{\text{forecast}}\).  
For each location, the procedure determines the days, \(t_{1}\) and \(t_{2}\), at which \(d(t)\) attains its first and second maxima. Then, the function representing the theoretical total number of deaths, \(F(t)\), is fitted against the reported total number of deaths, \(D(t)\), in the range between \(t=1\) and \(t=t_{\text{fit}}\). The procedure uses the MATLAB _lsqcurvefit_ function, which requires a minimum value, a guess value, and a maximum value for each parameter [19]. These are specified as detailed in Table 1.  
As a result, the following six parameters are determined for each country: \(A\), \(M_{1}=\exp\left(\mu_{1}-\sigma_{1}^{2}\right)\), \(\sigma_{t1}=\sqrt{\left(\exp\sigma_{1}^{2}-1\right)\exp\left(2\mu_{1}+\sigma _{1}^{2}\right)}\), \(M_{2}=\exp\left(\mu_{2}-\sigma_{2}^{2}\right)\), \(\sigma_{t2}=\sqrt{\left(\exp\sigma_{2}^{2}-1\right)\exp\left(2\mu_{2}+\sigma _{2}^{2}\right)}\), and \(\alpha\). Here, \(M_{1}\) and \(\sigma_{t1}\) respectively denote the mode and standard deviation of the first lognormal distribution. Likewise, \(M_{2}\) and \(\sigma_{t2}\) are the same quantities for the second distribution [18]. The obtained parameters are substituted into Equations (1) and (2). Thus, the theoretical estimates of the daily and total number of deaths are computed between \(t=1\) and \(t=t_{\text{max}}\). The procedure is repeated for each  
\begin{table}
\begin{tabular}{c c c c} \hline \hline
**Parameter** & **Minimum** & **Guess** & **Maximum** \\ \hline \(A\) (No.) & \(1.10\,D(t_{\text{fit}})\) & \(\max\left\{6\,D(t_{1}),\,2\,D(t_{\text{fit}})\right\}\) & \(\max\left\{12\,D(t_{1}),\,4\,D(t_{\text{fit}})\right\}\) \\ \(M_{1}\) (days) & \(1\) & \(t_{1}\) & \(t_{\text{max}}\) \\ \(\sigma_{t1}\) (days) & \(0\) & \(1\) & \(t_{\text{max}}\) \\ \(M_{2}\) (days) & \(1\) & \(t_{2}\) & \(t_{\text{max}}\) \\ \(\sigma_{t2}\) (days) & \(0\) & \(1\) & \(t_{\text{max}}\) \\ \(\alpha\) (–) & \(0\) & \(\frac{1}{1+d(t_{2})/d(t_{1})}\) & \(1\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: _lsqcurvefit_ range of fitting parameters.  
of the 211 locations. Hence, the predicted time trends for the six continents and the World are obtained by summing those computed for the belonging countries.  
The MATLAB script saves the results obtained for each location--in terms of fitted parameters, residual norm and exit flag returned by the _lsqcurvefit_ function, actual number of deaths (at \(t=t_{\text{fit}}\)), and forecast number of deaths (at \(t=t_{\text{forecast}}\))--into a CSV file called _COVID19_BiLog_World.csv_. Besides, for each (selected) country, continent, and the World, two plots are created and saved in PNG format: (i) the daily number of deaths and (ii) the total number of deaths, as functions of time.