The Monte Carlo simulations [8,9,10] have been designed to possibly take into account the measurement error in each daily number of the cumulative positive cases of COVID-19 in Italy. This error should describe the uncertainty in the process of measuring the daily number of positive cases due to fluctuations in the measurement procedures (such as a different number of performed daily nasopharyngeal swabs of one day with respect to another day); of course, this error does not describe the difference between the actual total positive cases and the diagnosed ones which can be very large [4]. However, the diagnosed cases are hypothesized to be a representative sample of the actual population (i.e., of the total number of positive cases, which is unknown). To get an estimate of the uncertainty in each daily number, we applied the following heuristic approach. We have assumed a measurement uncertainty in the total number of positive cases equal to 10% of each daily number (Gaussian distributed).

The second step was to generate a random matrix \(\left( {m\times n} \right) \), where n (columns) is the number of observed days and m (rows) is the number of random outcomes, which we have chosen to be 150. Each number in the matrix is part of a Gaussian distribution with mean equal to 1 and sigma equal to 0.1 (i.e., 10% of 1), both row-wise and column-wise. In such a way, each day will be characterized by 150 simulated outcomes that allow to apply a statistical approach. The 150 outcomes represent a reasonably large number of simulated deviations from the official data.

We then multiplied the nominal value of the diagnosed positive cases of the jth day for the 150 numbers of the jth column of the random matrix mentioned above. In such a way, each day will be associated with a series of 150 numbers (Gaussian distributed with a 10% standard deviation), which simulate the statistical nature of a single datum. The index j will run from February 15, 2020, to March 26, 2020. Finally, we integrated the daily data to obtain 150 series of the cumulative diagnosed positive cases that allowed to perform the statistical analysis.

In summary, starting from the n nominal values of the daily data, we generated n Gaussian distributions with 150 outcomes, with mean equal to the n nominal values and with 10% standard deviation. Then, for each of the 150 simulations, these n values (corresponding to the cumulative positive cases of n days) were fitted with a four-parameter function of the type of the Gauss error function, and we then determined the date of the flex with such fitted function for each simulation. Using the fitted function, we also determined the date at which the number of daily positive cases will be at a certain threshold that, for example, we have chosen to be 100. Finally, we calculated the standard deviation of these 150 simulations. In Fig. 5 and Fig. 6, we report the values (red dots) and the mean (horizontal solid line) of the Monte Carlo simulations, respectively, for the date of the flex and for the date of a substantial reduction in the number of daily positive cases.

Fig. 5 Monte Carlo simulations: each red dot corresponds to the day of occurrence of the flex (after which there is a reduction in the number of daily cases, i.e., there is a deceleration in the number of daily cases), obtained with each of the 150 Monte Carlo simulations. The vertical axis reports the number of days from February 15, 2020 Full size image

Fig. 6 Monte Carlo simulations: each red dot corresponds to the day in which a substantial reduction in the number of daily cases (about 100) occurs. Each red dot is obtained with each of the 150 Monte Carlo simulations. The vertical axis reports the number of days from February 15, 2020 Full size image

Using \(n = 41\hbox { days}\) (i.e., the number of daily diagnosed positive cases up to March 26, 2020), the mean of the 150 Monte Carlo simulations gives the expected dates of March 25, 2020, and April 22, 2020, for the flex and the day of a substantial reduction in the number of daily cases (about 100), respectively. We then obtained a standard deviation (1-sigma) of 1 day for the date of the flex and of 2.3 days for the date in which a substantial reduction in the number of daily cases would be about 100. This result corresponds to a probability of 68.2% that the date of the flex will be at a certain date plus or minus 1 day and that the date of a substantial reduction in the number of cases will be at a certain date plus or minus 2.3 days. A 2-sigma standard deviation will give a more robust probability of 95.4% of the day of the flex and of the day of a substantial reduction in the number of cases. The 2-sigma values correspond to plus or minus 2 days for the day of the flex and plus or minus 4.6 days for the day of a substantial reduction in the number of cases. A similar uncertainty was obtained for the day of the flex and the day of a substantial reduction in the number of fatalities.