In many countries across the world, there is an initial bottleneck in testing for COVID-19. This is due to a combination of two factors, first the availability of test kits, and second, of trained laboratory personnel to perform the tests. When the number of tests per million of population is limited, the exponential growth in infections at the beginning of an epidemic makes it impossible to estimate the prevalence of infections in any meaningful way. One aim of this paper is to establish and validate a Bayesian statistical inference procedure to make a probabilistic estimate of the infected population from the statistics of fatalities. This is done using data for Mumbai, which is one of the major COVID-19 epidemic hot spots in India.
The counts of fatalities due to COVID-19 also have errors. The most important is the fraction of mis-diagnosis. At the current time, it is unlikely that a patient who dies of severe acute respiratory illness (SARI) will be misdiagnosed. We know that the Municipal Corporation of Greater Mumbai (MCGM) tries to track, and report on, every such case, even if the diagnosis comes a few days after the fatality. Another source of systematic error in the count of fatalities may come from a part of the population with limited access to health care. This is not a trivial source of error, and corporations have to check their death records very thoroughly to make sure that no case is missed. In New York City, where the number of daily fatalities attributed to COVID-19 has been more than 100 since late March [1] the local government is wary of undercounts [2] . Similar undercounts have been reported from several European countries [3] , and a correction by 100-200% was found to be necessary. A retrospective correction in the record of fatalities is being carried out in Wuhan. Nevertheless, fatality counts may be in error by a factor 2-3, whereas testing could potentially have errors which are one or two orders of magnitude higher.
In view of this, it makes sense to estimate epidemic parameters from the daily count of fatalities, D(t), rather than less uncertain statistics such as the daily count of test positives T + (t). This is the primary aim of this paper. Daily fatality counts in Mumbai are currently less than 20. This is about one part in a million of the population, and an order of magnitude below peak rates in cities of comparable size elsewhere in the world. The simplifying assumption that the epidemic is still in the stage of exponential growth is then viable. The exponential growth rate can then be reliably extracted from D(t) without using any detailed model. Detailed analysis of the data also allows validation of this assumption, as we discuss in the results section of this paper.
Two more parameters, the infection fatality rate, IFR, and the interval from infection to death, T , relate the number of infections, I(t), and D(t). Since data on I(t) in Mumbai are currently unreliable, these two parameters are taken from the literature [4] to predict a credible interval for I(t). Alternative estimates of I(t) would, of course, constrain these parameters. First reports of such alternative estimates are now available from other locations [5] . We find that the ratio Π = T + (t)/I(t) is small, and has been almost independent of time in Mumbai. This ratio has also been reported in [5] . It is checked that the estimate of Π made here is compatible with that, provided that the different testing rates in different countries are taken into account.
These statistical inferences lead directly to estimates of the date at which the initial seeding of COVID-19 took place. The further assumption which goes into this estimate is that the growth rate has remained unchanged. Another More recent data has not been analyzed together with this, because retrospective corrections may be applied for those. Notice the log scale on the y-axis.
caveat should be kept in mind. Due to the statistical nature of infections, it is possible that multiple seedings occurred, and that the date obtained by the procedure followed here gives an average over such statistics. Stochastic models would refine our understanding of the early growth process. Genomic studies of multiple samples of the virus [6] in Mumbai would be independent data which could constrain such models. MCGM releases a daily count of fatalities due to COVID-19, and performs retrospective corrections due to late receipt of test results or other post-mortem analysis. Due to this, the count of fatalities may take a few days to settle down. Changes in municipal priorities may also affect resources assigned to various tasks, resulting in systematic changes in the statistical character of data streams. Any analysis of epidemic parameters must allow for checks of such confounding factors. Other data streams, such as the municipal records of deaths, are harder to access. However, such records will need to be subjected to statistical analysis later to ascertain whether there were possible cryptic fatalities due to COVID-19. The methods used here can also be utilized to analyze the results of such an appraisal.
Daily counts released by MCGM have been archived [7] . Here the summed and anonymous data on fatalities is analyzed for the two-week period between May 30 and April 14. The basic data on counts is displayed in Figure 1 . No clear trend can be distinguished in the first week, but a clear exponential growth, characteristic of the early stages of epidemics, can be seen in the second week.
During the early stages of an epidemic the number of fatalities is expected to grow exponentially, doubling in every τ days,
The parameter D 0 is an estimate of the number of fatalities at the time when the counting starts. Both τ and D 0 will be treated as parameters to be determined from the time series D(t) using a standard Bayesian inference process. Given the set of parameters π = {τ, D 0 }, and the set of data to be fitted {D 1 , D 2 , · · ·}, which are the fatalities on days 1, 2, etc., and the errors {σ 1 , σ 2 , · · ·}, the probability distribution function (PDF) of D i given π is defined as
where D(t) is the model of eq. (1), which depends on π. With the present disease surveillance regime, significant miscounts of D t would come only from small clusters of unnoticed infections, which would grow at roughly the same rate as the rest of the population. So the error would be a fraction of D t . We will take σ t = f D t , with two scenarios of f . In scenario L, the low error scenario, we take f = 0.1, i.e., assume 10% errors on D t . In scenario M, the medium error scenario, we double the errors and take f = 0.2. Bayes' theorem can be used to write the posterior PDF for π,
where P 0 (π) is called the prior PDF of the parameters The constant of proportionality, which is the prior PDF of the data, is trivial, since it does not depend on π, and serves only to normalize the posterior PDF. Maximizing the posterior PDF gives the best values of π. We will also quote the 95% credible intervals (95% CrI) obtained from the posterior distribution. The prior PDF is taken to be the product of two independent Gamma distributions, one for each of the parameters, p 0 (D 0 ) = Γ(D 0 ; k 0 , θ 0 ) and p 0 (τ ) = Γ(τ ; k, θ). We take the meta-parameter θ = θ 0 = 10, so that we don't suffer from overconfidence in priors [8] . The parameters k and k 0 are chosen to give the most probable prior values of τ = 10 days and D 0 = 5. It was checked that different choices of prior parameters gave results which were insensitive to the choices.
C. Mathematical formulation: estimation of number of infections I 0 , the number of infections on day zero, can be inferred from D 0 by a straightforward argument [9] . First, the infection fatality rate, IFR, gives a statistical estimate of the number of infections, I, which resulted in the observed fatalities on day zero. This is simply I = D 0 /IFR. The value of IFR could change from country to country because of differing age and gender structure of populations, or due to prevailing health conditions. However, a larger source of uncertainty is in the fraction of asymptomatic cases. As a result, IFR estimates vary from over 1% [10] to under 0.5% [11] . Based on an complete sampling of a well-defined population, it has been estimated that the age and gender averaged IFR is 0.657% [4] . This is the estimate that we use here. A recent preprint [12] updates the work in [9] by correcting this value using Indian census data to obtain a value which is approximately two thirds of this. Using this number would boost our median estimates by about 50%. Since the 95% CrI are much larger, we do not incorporate this correction.
From I it is possible to construct I 0 using the previous estimate of τ to write I 0 = 2 T /τ I. Also, T = t 1 + t 2 where t 1 is the incubation period, i.e., the interval between infection and the occurrence of symptoms, and t 2 is the interval between the occurrence of symptoms and death. Several independent estimates agree that the median value of t 1 is around 5 days [13, 14] . For t 2 we use the mean of 18.8 days and coefficient of variation of 0.25 [4] .
Putting these together we have
and the distributions of the three parameters are taken to be Γ(IFR; 6.3, 0.0011), Γ(t 1 ; 2.06, 2.93) and Γ(t 2 ; 16, 1.175). The amplification factor, within square brackets, has all the uncertainties associated with these inputs, unlike the parameters D 0 and τ . Nevertheless, I 0 is of high interest for public health reasons, so we will present estimates.
A. March 31 to April 4
The first analysis was performed with the short time series for fatalities between March 31 and April 4. The most likely model that emerges out of the data is 
The joint 95% and 99% CrI on the model parameters are shown in Figure 2 . Note that the posterior distribution of the model parameters is strongly non-Gaussian. From the posterior PDF for τ , marginalized over D 0 in scenario M, we find τ to be 8.75 days (median values) and [3.2, 30] days 95% CrI, giving a significant probability of very long doubling times. From the PDF for D 0 , marginalized over τ , we find a median of 0.98 and [0.6, 1.3] 95% CrI. One sees in Figure 2 that there is a broad tail to the posterior PDF which allows slightly larger values of D 0 to compensate for a much slower epidemic growth. Longer time series are needed to rule this out, as will be seen later. In Scenario L, after marginalizing over D 0 , τ is found to be 9.0 days (median) with [4.8, 27 ] days 95% CrI. For D 0 it was seen that D 0 = 1.02 and [0.8, 1.2] 95% CrI. The wide difference between the 95% CrI and the 99% CrI curves in scenario L is an indication that 10% Gaussian errors in scenario L are unable to give a good description of the jitter in the time series. In view of this, scenario L is not used in subsequent analyses.
The second part of the analysis used the data for the eight days between April 6 and 14, with day zero to be April 2. So D 0 and I 0 in this context are estimates of the inferred number of infected cases in Mumbai on April 2. Note that this time series is almost 50% longer than the first, and therefore regulates some of the peculiarities seen earlier. The most likely model parameters in Scenario M are τ = 5.1 days, D 0 = 2.5, which implies I 0 = 44000 [5500, 960000] 95% CrI. The 95% and 99% CrI on the joint distribution of model parameters are shown in Figure 3 . Although the contours are not as distorted as in Figure 2 , they are clearly not ellipses. So the posterior PDF is distinctly non-Gaussian. After marginalizing over D 0 , we find τ is 5.6 days (median) and [ Given the range of τ quoted above, in 3 days between March 31 and April 2, one expects D 0 on May 31 to increase by a factor [1.3, 1.7] 95% CrI, and become [1.0, 2.0] 95% CrI. While this interval overlaps with that quoted above, the probability that the two estimates are statistically indistinguishable is about 3.8%. The reason for this discrepancy can be traced to a change in the surveillance policy instituted on April 2 [15] . The increased vigilance after this date must be corrected for in order to compare the two data sets. One can account for this by changing the overall normalization of the earlier data set by a factor of 1.7. This then increases the probability of agreement between the two analyses by an order of magnitude.
If one restricts the time series between April 6 and 12, one sees "by eye" that there could be rapid exponential growth. For this specially selected data set, we find that the most probable τ is 2.8 days. There is only a 20% probability that τ is as large as 4 days or more. These alarming results are biased. It is important to remember that infection spreading is a stochastic process, and there are always short runs of low or high growth.
Third, the analysis of the full data set is presented. A reporting discrepancy between the series up to April 4 and the later part was discussed in the previous section. We applied the correction developed there, and rounded D(t) to the nearest integer. This corrected time series gave τ = 5.3 days, D 0 = 1.8, which implies I 0 = 27000 [3500, 520000] 95% CrI.
The 95% CrI of the joint distribution is shown in Figure 4 . The PDF for τ marginalized over D 0 gives τ = 5.5 days (median) with [4.6, 6.9] days 95% CrI, and that for D 0 marginalized over τ gives D 0 = 1.86 (median) with [1.45, 2.35] 95 . These are our primary results.
We have noted earlier that retrospective changes are expected often to the most recent fatality data. This is why the analysis of the data for April 15-19 is made separately. Taking day zero to be April 14, we find τ = 2.6 days, and D 0 = 1.59. The joint 95% CrI of the parameters is shown in Figure 5 . The doubling time is half of that in eq. (7), implying an uncontrolled growth rate for the epidemic. There is less than a 11% chance that this doubling rate is compatible with that in eq. (7). Even more disturbing is that D 0 is very low, with a negligible chance (less than one part in 10 −30 ) that it is compatible with the lower bound on the 95% CrI of fatalities allowed by eq. (7) on April 14. It is clear from this that the data is either not yet complete, or that we are confronted with a huge statistical deviation away from the previous part of the same time series.
A. The doubling time
The time series of data on fatalities in Mumbai was used to extract several characteristics of the epidemic in Mumbai. The most robust parameter that comes out of this is the doubling time, τ = 5.3 days (median, with [4.6, 6.9] 95% CrI). It is surprising that this rate of growth for COVID-19 is seen during the lock down. It is comparable to the doubling times reported in the early days of the epidemic in Wuhan, before a quarantine was imposed there [16] .
We also infer a parameter D 0 which is the number of estimated fatalities on day zero of the analysis, which we take to me March 30. We find D 0 = 1.8 (95% CrI, 1. 4-2.4) . This is in good agreement with the reported number of 2 on May 30. We note that this agreement indicates that the stochasticity in daily numbers is weak. One supporting evidence is that the cumulative count of fatalities has grown very smoothly over this period with the same doubling time as found here. This is a direct check that the assumption of exponential growth remains valid.
The doubling time can be matched to the basic reproduction rate, R 0 , in different epidemiological models. Here we match it to a simple SEIR model. This has been used widely to model the COVID-19 epidemic, since the disease is non-infective in the incubation state. In the early stages of an epidemic, when the fraction of susceptibles in the population is close to unity, one may write
where the distribution of t 2 used is given after eq. (4). Using the value of τ quoted in eq. (7), we find R 0 = 3.4 (median) and [2.4, 4.8] 95% CrI. This is consistent with estimates from a range of countries, but closer to the upper edge of the global spread in values. Note that the only inputs here are the doubling time τ and the case resolution time t 2 . There are public health implications which we discuss later. 
The process of extracting I 0 is indirect, since it depends on three different epidemiological parameters which require studies with patient databases more extensive than available in Mumbai now. We have taken them from [4] .The 95% CrI shown in Figure 6 should be understood as containing these uncertainties according to eq. (4).
MCGM's daily counts of test positives, T + (t), for Mumbai over the period of March 31 to April 14 were analyzed with exactly the same methods which were applied to D(t). It was found that the doubling time, τ , was 5.7 days (modal value). The modal number of test positives on March 30 was inferred to be T 0 = 33. Given the estimated I 0 on that day eq. (7), Π = T + (t)/I(t) can easily be estimated. It was found that on day zero Π = 0.13% (median). Since the doubling time of the numerator and denominator are the same, this ratio has remained constant through the epidemic.
Without a massive increase in the number of tests, there is no direct way to verify the estimates of infections presented here. Several alternatives have been proposed, ranging from testing viral load in sewage [17] to conducting sample surveys using different testing methods. An example of the latter is a recent serological sample study [5] in California, USA. This study reported Π to be 1.1-2%. When the test rate, i.e., the number of tests administered per million of population is small, Π is proportional to the test rate. On April 18, India's cumulative test rates were 2.69 × 10 −4 , USA had test rates of 1.1159 × 10 −2 [18] . If the Indian test rate is scaled up to the same value as the US, then, with the estimates of Π presented in the previous paragraph would have to be multiplied by a factor of 42. This scaled value of Π for Mumbai is compatible with that reported for California in [5] .
Given τ = 5.3 days, we can extrapolate from I 0 backwards to find when there was only a single case in Mumbai. This estimates the first seeding of COVID-19 in the city. We find that the median date for the initial seeding was around January 15, 2020. However, there is a 95% chance that the initial seeding occurred as early as December 20, 2019 or as late as February 1, 2020. Although the volume of traffic from China to Mumbai is not large, it is steady, and early seeding is not inconceivable, especially if the virus was spreading inside China already in mid-November. The first fatalities due to this cryptic spread of COVID-19 would not have happened before mid-February (median estimate). It is certainly possible to check municipal records of deaths in February to see whether any excesses over the same month in previous years are visible. If there is a statistically significant signal, it would support the idea of early cryptic seeding of COVID-19 in Mumbai.
It is interesting that the value of R 0 deduced from the growth rate during a lock down is among the higher end of values seen in other countries. The purpose of this drastic measure was to reduce the value of R 0 . This could be an indication that the epidemic growth is being driven by high density areas where physical distancing is impossible. It would be important to understand the structure in finer detail, ward-wise or at even finer scales, if possible. If there is indeed such geographical structure to the epidemic within Mumbai, then it is possible that post-lock down there could be a "second wave" surge of cases in areas where housing density is lower. Active steps must be taken to prevent this. The Maharashtra Government order of April 18, removing all relaxations of lock down in the Mumbai and Pune regions [19] may be seen in this context.
It should be noted that a significant fraction of the infected could be asymptomatic [20] [21] [22] . In this connection note that the symptom fatality ratio, SFR, defined as the probability of dying after the onset of symptoms is reported to be 1.4% (median) [0.9%, 2.1%] 95% CrI [16] . Comparing this to the IFR of 0.657% gives an estimate of the probability of being asymptomatic: P A = 1 − IFR/SFR. This is about 50%, and in the range quoted in literature. So it is possible that half the cases are not of significance for public health measures [22, 23] . Symptomatic cases are of significance. For these contact tracing, where possible, and quarantine, otherwise, remain the best tools available now for controlling spread by bringing down the growth rate.
However, the asymptomatic fraction is of importance for a long term solution to the COVID-19 epidemic. Barring a sudden discovery and large-scale deployment of a cure or vaccine, the epidemic can be only be stopped when the number of infected reaches a level where herd immunity begins to operate. Asymptomatic infections have to be counted in this. From this point of view one should see the upper ends of the bands in Figure 6 as the most optimistic scenarios for reaching this desirable end state to the epidemic.
Note that if the growth rate is not controlled quickly, the situation could become much worse. With the estimate of R 0 made here, it should be expected that herd immunity sets in when 71% (median and [58%, 79%] 95% CrI) of the population has been infected. Given IFR = 0.657%, this could imply around 100,000 fatalities in Greater Mumbai until herd immunity is reached. Note that conclusions like this may be model dependent. This extremely simple model is quoted here only for orientation.
Another way to look at these predictions is the following. The fatality rates till now have not been severe. However, at the current rate of doubling, it would take about 15 days for them to climb about 10-fold. This increase would give about 100 fatalities a day, and could climb above that. This is in the range of the peak number of daily COVID-19 fatalities in New York City [1].
If the fatality rate increases to such levels, the number of cases requiring critical care increases even more. Also critical care lasts longer for cases which recover. In the hypothetical cases discussed here, COVID-19 wards in the city would have to deal with several hundred admissions a day with care periods for each case ranging up to two months. Again, using the estimates given above, these rates would be reached in two to four weeks if the rate of growth of the epidemic does not come down.
I would like to thank Basudeb Dasgupta, Sandeep Krishna, S. Krishnaswamy, Rajdeep Sensarma, and R. Shankar for careful readings of the manuscripts and useful suggestions, and the ISRC (Indian Scientists' Response to COVID-19) mailing list for generating questions.