Coronavirus disease 2019 (COVID- 19) , which is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) (1), was detected in the United States (US) in January 2020 (2) . In February, COVID-19-caused deaths were detected (3) . Thereafter, surveillance testing expanded nationwide (4) . These and other efforts revealed community spread across the US and exponential growth of new COVID-19 cases throughout most of March with a doubling time of 2 to 3 d (5) , similar to that of the initial outbreak in China (6) . This situation led to government mandates closing schools, prohibiting public gatherings, and restricting commercial activities, as well as broad adoption of social-distancing practices, such as working-from-home, curtailing of travel, and face mask-wearing (7) . Although the US became a hot spot of the COVID-19 pandemic, detection of new cases peaked in late-April and steadily declined until mid-June (4), suggesting that public-health mandates and social-distancing practices were effective at slowing COVID-19 transmission. Attempts to quantify the impacts of these measures suggest substantial benefits in terms of reduction of disease burden (8, 9) .
At present, the number of new daily cases in the US is again increasing (4) . It is imperative that we effectively monitor ongoing COVID-19 transmission, so that dangerous upticks in cases can be responded to as quickly as possible.
To contribute to situational awareness of COVID-19 transmission dynamics, we developed two distinct mathematical models for the regional COVID-19 epidemic in each of the 15 most populous US metropolitan statistical areas (MSAs) (10) . The two models are 1) a fitting function, which can reproduce the shape of an epidemic curve having two timescales (e.g., fast growth and slow decay), and 2) a compartmental model. The latter model is composed of ordinary differential equations (ODEs) characterizing the population dynamics of susceptible, presymptomatic, symptomatic, asymptomatic, and recovered populations and various subpopulations, including mixing and protected-by-social-distancing populations, individuals who are quarantined or self-isolated, individuals incubating virus without symptoms at multiple stages of disease progression, and individuals with mild and severe forms of disease. The model also tracks hospitalizations and deaths.
In an ongoing online-learning effort, we are calibrating our models for regional COVID-19 epidemics on a daily basis for consistency with historical case reports. We have also been applying Bayesian methods to quantify uncertainties in predicted detection of new cases. In the face of variability in case detection, this approach allows for early objective identification of new epidemic trends beyond what can be reasonably explained by a well-calibrated model that is consistent with historical case reporting, enabling an evidence-based response by policy makers.
The COVID-19 surveillance data used to parameterize models-reports of new confirmed cases, the earliest widely available indicators of local trends-are obtained daily (at variable so as to allow consideration of multiple periods of social distancing with distinct setpoints for the protected population size. In the model, there is always an initial period of social distancing. The number of additional social-distancing periods is indicated by . Here, we only consider two cases: = 0 and = 1.
For = 0 (i.e., only one social-distancing period), the compartmental model and auxiliary measurement model have a total of 20 parameters. We take six of these parameters to have adjustable values (Table 2 ) and 14 to have fixed values (Table 3 ). In the Appendix, we describe each parameter and explain the fixed parameter settings, which are based on information in Refs. (12) - (20) and assumptions. The adjustable model parameters are as follows: ! , the start time of the local epidemic; > ! , the time at which social distancing begins; ! , which establishes a setpoint for the quasi-stationary fraction of the total population practicing social distancing; ! , which establishes a timescale for adoption of social-distancing practices; and , which characterizes the rate of disease transmission in the absence of social distancing. The measurement-model parameter, " , represents the time-averaged fraction of new cases detected, which is related to the intensity of local surveillance efforts. As for the curve-fitting model, inference of adjustable parameter values is based on a negative binomial likelihood function (Equation (31) in the Appendix). The dispersal parameter of the likelihood is taken to be adjustable; its value is jointly inferred with those of ! , , ! , ! , , and " . All inferred parameter values are inference-time-dependent and region-specific. Inferences are conditioned on the model and fixed parameter estimates.
For > 0 (i.e., distinct social-distancing periods following an initial or first-phase socialdistancing period), the compartmental model has three additional adjustable parameters for each additional period of social distancing considered beyond the initial period. For one additional period of social distancing ( = 1), the additional adjustable parameters are as follows: # > , the onset time of second-phase social-distancing; # , the second-phase quasi-stationary setpoint parameter; and # , which determines the timescale for transition from first-to second-phase social-distancing behavior. If adherence to effective social-distancing practices begins to relax at time = # , then # < ! .
The curve-fitting and compartmental models are both deterministic, despite COVID-19 transmission being a stochastic process. We interpret each model to predict the expected number of new confirmed COVID-19 cases reported daily by the relevant public health authorities. In other words, we assume that the number of new cases reported over a 1-d period is a random variable and its expected value follows a deterministic trajectory. We further assume that day-today fluctuations in the random variable are independent and characterized by a negative binomial distribution, which we will denote as NB( , ). We use NB( , ) to statistically model noise in reporting (and case detection, which influences reporting) because its support-the non-negative integers-is natural for populations and, in addition, its shape is flexible enough to recapitulate a wide range of unimodal empirical distributions. With these assumptions, we obtain a likelihood function (Equation (31) in the Appendix) taking the form of a product of probability mass functions of NB( , ). Formulation of a likelihood is a prerequisite for standard Bayesian inference 2 .
We use Bayesian inference to learn adjustable model parameter values consistent with the time-series of daily MSA-specific reports of new COVID-19 cases available up to the point of inference. Inferences with both models are performed daily for each MSA of interest. In each inference, we assume a uniform prior and use an adaptive Markov chain Monte Carlo (MCMC) algorithm (21) to generate samples of the posterior distribution for the adjustable parameters. A full description of our inference procedure is provided in the Appendix.
The maximum a posteriori (MAP) estimate of a parameter value is the value of the parameter corresponding to the mode of its marginal posterior, where probability mass is highest.
In addition to inferring parameter values, we quantify uncertainty in predicted trajectories of daily case reports. A predictive inference of the expected number of new cases detected on a given day is derived from a model by parameterizing it using a randomly chosen parameter posterior sample generated in MCMC sampling. The number of these cases that are detected is then predicted by adding a noise term, drawn from a negative binomial distribution NB( , ), where is set at the randomly sampled value-recall that the value of is jointly inferred with adjustable model parameter values-and is set such that the mean of NB( , ) corresponds to the predicted number of new active cases, i.e., the value of is set using Equation (30) in the Appendix. A full description of our predictive inference procedure is provided in the Appendix.
Predictions of the curve-fitting model are obtained by evaluating the sum in Equation (2) in the Appendix. Predictions of the compartmental model are obtained by using LSODA (22) to numerically integrate Equations (5)- (21) and (27) in the Appendix; the initial condition is defined by the inferred value of ! (Table 2 ) and the fixed settings for ! and ! ( Table 3 ).
The 95% credible interval for the predicted number of new case reports on a given day is the central part of the marginal predictive posterior capturing 95% of the probability mass. This region is bounded above by the 97.5 percentile and below by the 2.5 percentile
The objective of our ongoing study is to detect significant new trends in new COVID-19 cases as early as possible by 1) systematically and regularly updating mathematical models capturing historical trends in regional COVID-19 epidemics through Bayesian inference and 2) making forecasts with rigorously quantified uncertainties through Bayesian predictive inference.
An important aspect of how we are analyzing COVID-19 data is our focus on the populations of US cities and their surrounding areas (i.e., MSAs) vs. the regional populations within other political boundaries, such as those of the US States. The boundaries of MSAs are defined on the basis of social and economic interactions (10) , which suggests that the population of an MSA is likely to be more uniformly affected by the COVID-19 pandemic than, for example, the population of a State. In accordance with this expectation, daily reports of new COVID-19 cases for New York City (Figure 3 , panel A) are more temporally correlated than for the US State of New York (Figure 3, For each of the 15 most populous MSAs in the US, we parameterized a curve-fitting model and a compartmental model using MSA-specific surveillance data, namely, aggregated Countylevel reports indicating the number of new confirmed COVID-19 cases within a given MSA each day. Bayesian parameterization and forecasting with uncertainty quantification (UQ)-predictive inference-are now being performed daily for each of the 15 MSAs, with both models.
Results of Bayesian predictive inference are exemplified in Figure 4 . Predictions are obtained in the form of a predictive posterior distribution and are based on all case reports available at the time of inference. In Figure 4 , the entire shaded region indicates the 95% credible interval of the predictive posterior as a function of time. In other words, this band indicates where 95% of predictions of daily case reports fall at the times indicated. Predictions vary because of the uncertainties in adjustable model parameter estimates, which are characterized quantitatively through Bayesian inference. The colors within the shaded band indicate other credible intervals (10%, 20%, etc.) and also the median of all predictions as a function of time.
Predictive inferences for all 15 MSAs of interest are shown in Figures 5 and 6. The predictions of Figure 5 are conditioned on the curve-fitting model, and the predictions of Figure  6 are conditioned on the compartmental model. These results demonstrate that the curve-fitting and compartmental models are each capable of reproducing empirical epidemic curves for multiple MSAs, which vary in shape. However, the curve-fitting model is less flexible, as can be seen by comparing its predictions to those of the compartmental model for the Atlanta MSA, where there is high variability in the daily number of new cases detected. Although there is no clear downward trend in the data, the curve-fitting model nevertheless predicts a peak in late-April/early-May and a downward trend ever since. This prediction is obtained because the model, by design, is only capable of generating single-peaked epidemic curves that rise and then fall. Hereafter, we will focus on the compartmental model.
Recall that predictive inferences are made daily. In Figure 7 , we show predictive inferences for New York City and Phoenix made over a series of progressively later dates. These results illustrate that accurate short-term predictions are possible but continual updating of parameter estimates is required to maintain accuracy. In Videos 1 and 2, daily predictions for New York City and Phoenix, respectively, are shown as animations.
In practice, we find that the adjustable parameters of the compartmental model have identifiable values 3 , meaning that their marginal posteriors are unimodal. This finding is illustrated in Figure 8 , which displays a 7 × 7 matrix of 1-and 2-dimensional projections of the 7-dimensional MCMC posterior samples underlying predictive inferences for New York City. Plots of marginal posteriors are shown on the diagonal extending from top left to bottom right. The other 2-dimensional plots reveal correlations between pairs of parameter estimates (if any). The significance of identifiability is that, despite uncertainties in parameter estimates, we can expect predictive inferences of daily new-case reports to cluster around a central trajectory.
How does learning the region-specific adjustable parameter values of the compartmental model and a subsequent predictive inference (i.e., a forecast with UQ) improve situational awareness? In the vast majority of cases, when we forecast with UQ, the empirical new-case count for the day immediately following our inference (+1), and very often for each of several additional days, falls within the 95% credible interval of the predictive posterior. When the reported number of new cases falls outside the 95% credible interval and above the 97.5% percentile, we interpret this event, which we will call an upward-trending rare event, to have a probability of 0.0275 or less assuming the model is both explanatory (i.e., satisfactorily consistent with historical data) and predictive of the near future. If the model is predictive of the near future, the probability of two consecutive rare events is far smaller, less than 0.001. Thus, consecutive upward-trending rare events, which we will refer to as an upward-trending anomaly, can be reasonably taken as a sign that the model is not in fact predictive. Indeed, an anomaly suggests that the rate of COVID-19 transmission has increased beyond what can be explained by the model.
For New York City, anomalies are not seen, as illustrated in Figure 9 , panel A. However, for Phoenix, there are several recent anomalies, which preceded rapid and sustained growth in the number of new cases reported per day in June (Figure 9 , panel B).
In an attempt to explain the anomalies, we allowed the compartmental model to account for a distinct second social-distancing period, i.e., we increased the setting for from 0 to 1. With this change, the number of adjustable parameters increases from 7 to 10, because in the model, each distinct social-distancing period is characterized by 3 parameters, which must be inferred. Recall that these parameters are an onset time, rate parameter, and setpoint parameter. The value of the setpoint parameter determines the quasi-stationary fraction of the total regional population that is adhering to effective social-distancing practices. In the model, ! is the setpoint parameter of the initial social-distancing period starting at time = > ! , and # is the setpoint parameter of a distinct subsequent social-distancing period starting at time = # > . As can be seen by comparing the plots in Figure 10 , panels A and B, the compartmental model with two socialdistancing periods ( Figure 10 , panel B) better explains the Phoenix surveillance data available than the compartmental model with just one social-distancing period ( Figure 10 , panel A) 4 . Furthermore, the MAP estimate for # (~0.38) is less than that for ! (~0.49) (cf. panels C and D, Figure 10 ), and the marginal posteriors for these parameters are largely non-overlapping ( Figure 10 , panel D). These findings suggest that the recent increase in COVID-19 cases in Phoenix can be explained by relaxation in social distancing, which is quantified by our estimates for ! and # . The MAP estimate of the start time of the second period of social distancing corresponds to 24-May-2020 5 . Intriguingly, 8 of the 9 anomalies noted in Figure 10 , panel B occurred after this period, with the first of these occurring on 02-June-2020.
We considered the hypothesis that a one-time event (viz., a mass gathering) generating 1000's of new infections might trigger a new upward trend in COVID-19 transmission. Simulations for New York City and Phoenix do not support this hypothesis (Appendix Figure 2 ).
Besides Phoenix, four other MSAs have recent trends explainable by relaxation of social distancing (Appendix Figure 3 and Appendix Table 1 ). Upward-trending anomalies were detected for these MSAs (Appendix Figure 4 , panels A-D), but not for three of four other MSAs having epidemic curves consistent with sustained social distancing (Appendix Figure 4 , panels E-H). Daily predictions for the MSAs considered in Appendix Figure 4 are shown in Videos 3-10.
Daily online learning of model parameter values from real-time surveillance data is feasible for mathematical models for COVID-19 transmission. Furthermore, predictive inference of the daily number of new cases reported is feasible for the regional COVID-19 epidemics occurring in multiple US metropolitan areas. Our model-based analyses are ongoing and daily forecasts are being disseminated (23).
We have suggested how our predictive inferences can be used to identify harbingers of future growth in COVID-19 transmission rate. Two consecutive upward-trending rare events (i.e., instances where the number of new cases reported is above the upper limit of the 95% credible interval of the predictive posterior) seem to indicate potential for increased transmission in the future. This signal of future growth is perhaps especially strong when anomalies are accompanied by increasing prediction uncertainty, as is the case for Phoenix (Figure 9 , panel B).
We find that the recent increase in rate of transmission of COVID-19 in the Phoenix metropolitan area can be explained by a reduction in the percentage of the population adhering to effective social-distancing practices 6 ( Figure 10 , panel D), from ~49% to ~38%. Relaxation is inferred to have begun around 24-May-2020 ( Figure 10 , panel B). Recent upward trends in the rate of COVID-19 transmission in the Houston, Miami, San Francisco, and Seattle metropolitan areas can also be explained by relaxation of social distancing (Appendix Figure 3 and Appendix Table 1 ). These findings are qualitatively consistent with earlier studies indicating that social distancing is effective at slowing the transmission of COVID-19 (7, 8) , and encouragingly, they suggest that the future course of the pandemic is controllable, especially with accurate recognition of when stronger nonpharmaceutical interventions are needed to slow COVID-19 transmission. 
None.
Dr. Lin works as a Scientist in the Information Sciences Group of the Computer, Computational, and Statistical Sciences Division at Los Alamos National Laboratory. His primary research interest is development and application of advanced data-science methods in modeling of biological systems. -parameter of log-normal distribution -parameter of log-normal distribution Probability parameter of NB( , ) Dispersal parameter of NB( , ) ! 35* (d) 6.6* 7.9* 1.6** 0.42** Constrained*** 4.4* *Estimates of the adjustable parameters ( , ! , , , and ) are region-specific and inferencetime-dependent. Inferences are performed daily. Here, we report the maximum a posteriori (MAP) estimates inferred for New York City using all confirmed COVID-19 case-count data available in the GitHub repository maintained by The New York Times newspaper (11) for the period starting on 21-January-2020 and ending on 21-June-2020 (inclusive dates). Time = 0 corresponds to 0000 hours on 21-January-2020. **Estimates of the fixed parameters and are those of Lauer et al. (12) . These parameter estimates define a log-normal distribution that reproduces the empirical distribution of waiting times for the onset of symptoms after infection with SARS-CoV-2. ***The value of is constrained, i.e., its reporting-time-dependent value is determined by a formula, which is given by Equation (30) in the Appendix. Table 2 . Inferred values of the adjustable parameters of the compartmental model ( ! , , ! , ! , and ), the auxiliary measurement model ( " ), and the associated statistical model for noise in case detection and reporting ( ).
Estimate 0.12 12 *All estimates are region-specific and inference-time-dependent. Inferences are performed daily. Here, we report the maximum a posteriori (MAP) estimates inferred for the New York City MSA using all confirmed COVID-19 case-count data available in the GitHub repository maintained by The New York Times newspaper (11) for the period starting on 21-January-2020 and ending on 21-June-2020 (inclusive dates). Time = 0 corresponds to 0000 hours on 21-January-2020. **The probability parameter of NB( , ) is constrained, i.e., its reporting-timedependent value is determined by a formula, which is given by Equation (30) in the Appendix. apply to all regions of interest except for , the number of distinct social-distancing periods that follow an initial social-distancing period, and Figure 1 . Illustration of shapes that can be produced by the fitting function that we are using to capture trends in regional COVID-19 epidemic curves. The curve-fitting model is formulated such that it has the capacity to reproduce the shape of an epidemic curve having two timescales. Figure 2 . Illustration of the populations and processes considered in a mechanistic compartmental model for the dynamics of COVID-19 transmission. The model accounts for susceptible (S), exposed (E), asymptomatic (A), symptomatic (I), hospitalized (H), recovered (R), and deceased (D) populations. It also accounts for social distancing, which establishes mixing and protected subpopulations, quarantine driven by testing and contact tracing, and self-isolation spurred by symptom awareness. The incubation period is divided into 5 stages, which allows the model to reproduce an empirically determined Erlang distribution of waiting times for the onset of symptoms after infection (12). The exposed population (consisting of individuals incubating virus) includes presymptomatic and asymptomatic individuals. The -populations consist of true asymptomatic individuals in the immune clearance phase. The gray background indicates the populations that contribute to disease transmission. An auxiliary measurement model (Equations (27) and (28) in the Appendix) accounts for imperfect detection and reporting of new cases. with consideration of an initial period of social distancing followed by a distinct period of relatively lax adherence to social-distancing practices ( = 1) for the Phoenix metropolitan statistical area (MSA). In panel (C), the marginal posterior for the socialdistancing setpoint parameter ! inferred in the analysis of (A) is shown. In panel (D), the marginal posteriors for the social-distancing parameters ! and # inferred in the analysis of (B) are shown. Model selection indicates that the two-phase model is to be preferred (Appendix Table 1 ). In this analysis, we used data available from 21-January-2020 to 18-June-2020 (inclusive dates).
Appendix Table 1 . Strength-of-evidence comparison of compartmental models accounting for an initial social-distancing period only ( = 0) and for an initial social-distancing period followed by a distinct second-phase social distancing period ( = 1).
ΔAIC* ΔBIC* ! ,-! (95%)** (24), we interpret ΔAIC>10 to indicate "no support for = 0" and ΔAIC<-10 to indicate "no support for = 1." There are five MSAs for which ΔAIC and ΔBIC are both greater than 10: Houston, Miami, Phoenix, San Francisco, and Seattle. There is one MSA for which ΔAIC and ΔBIC are both less than −10: Boston. **The first entry in each row of this column is the maximum a posteriori (MAP) estimate. The next entry, a pair of numbers within parentheses, indicates the 95% credible interval. In this analysis, we used data available from 21-January-2020 to 26-June-2020 (inclusive dates).
Appendix Figure 1 . Detailed diagram of the populations and processes considered in the mechanistic compartmental model. In this illustration of the compartmental model, the labels attached to arrows indicate the parameters and variables that affect the rates of the processes represented by the arrows. There is a one-to-one correspondence between arrow labels and terms on the right-hand sides of Equations (5)- (21) . The diagram is otherwise the same as that shown in Figure 2 . Figure 2 . We evaluated the potential impact of a one-time mass gathering that causes 50,000 and 5,000 individuals to become newly infected on 30-May-2020 in (A) the New York City metropolitan statistical area (MSA) and (B) the Phoenix MSA, respectively. As can be seen, such an event will step the epidemic curve up without significantly changing the slope of the curve. According to the model, the trend in the slope of the curve is determined by the sustained level of adherence to effective social-distancing practices. We conclude that the recent trajectory of the epidemic curve of the Phoenix MSA cannot be explained by a one-time mass gathering.
Appendix Figure 3 . Predictive inferences conditioned on the compartmental model with either (A)-(E) one or (F)-(J) two distinct social-distancing periods. We consider five metropolitan statistical areas (MSAs), which according to Appendix Table 1 , have epidemic curves better explained by the compartmental model with two social-distancing periods than by the compartmental model with just one. Maximum a posteriori (MAP) estimates for # indicate that the second social distancing period began on 27-May-2020 for Houston, 19-April-2020 for Miami, 24-May-2020 for Phoenix, 12-June-2020 for San Francisco, and 07-June-2020 for Seattle. In this analysis, we used the data available from 21-January-2020 to 26-June-2020 (inclusive dates). Table 1 ). Upward-trending anomalies were detected for each of these MSAs. Los Angeles, Chicago, Dallas, and Washington, DC have epidemic curves better explained by the compartmental model with just one social-distancing period vs. two (Appendix Table 1 ). Upward-trending anomalies were not detected for three of these MSAs. The exception is Chicago. Two anomalies were detected for Chicago near the end of the period considered in the analysis. 
For each metropolitan statistical area (MSA) of interest, we assume that there is an infection curve ( ) describing the number of individuals who become infected at time with SARS-CoV-2 and who will later be detected in local COVID-19 surveillance efforts. Furthermore, we assume that this curve has a shape that can be generated/reproduced by . ( , , ) , the probability density function (PDF) of a gamma distribution Γ( , ). In other words, we assume that ( ) = . ( , , ) , where is a scaling factor that we can identify as the number of individuals who will be detected over the entire course of the local epidemic. The shape of a gamma distribution is flexible and determined by the values of its two parameters: , which is called the shape parameter, and , which is called the scale parameter. The functional form that we assume for ( ) allows the curve-fitting model to reproduce the shape of an epidemic curve having two timescales. Many empirical COVID-19 epidemic curves appear to have two timescales: an initial period during which new case reports increase relatively quickly from day to day followed by a period during which new case reports decrease relatively slowly from day to day.
We do not take the infection curve ( ) to correspond directly to the number of new COVID-19 cases reported on the date encompassing time , because only symptomatic individuals are likely to be detected in COVID-19 surveillance testing (to a first approximation). This situation complicates our model as there is known to be a potentially lengthy, variable delay in the onset of symptoms after infection (12) . We assume that the waiting time − for the onset of COVID-19 symptoms after SARS-CoV-2 infection at time is distributed according to a log-normal distribution. Let us use '/ ( − ; , ) to denote the PDF of the waiting-time distribution modeled by a log-normal distribution with parameters and set to the values estimated by Lauer et al. (12) . Let us use ( 0 , 01# ) to denote the predicted number of new COVID-19 cases reported within a period beginning at time 0 ≡ ! + d and ending at time 01# , where ! > 0 is the start time of the local epidemic. We assume that surveillance testing for SARS-CoV-2 infection starts prior to time ! , and we take time = 0 to correspond to 0000 hours on 21-January-2020, the date on which detection of the first US COVID-19 case was widely reported (3) . Under the aforementioned assumptions, ( 0 , 01# ) is given by a convolution of integral functions. Namely, ( 0 , 01# ) is given by the following expression:
It should be noted that in this expression is a dummy variable of integration.
Equation (1) can be evaluated through numerical quadrature, but this procedure is computationally expensive. To overcome this limitation, we replace the double integral in Equation (1) with a sum, and we calculate ( 0 , 01# ) using the following expression instead of Equation (1):
and
In Equation (2), the 045 terms are weighting functions (i.e., kernels) that account for the variable duration of the incubation period, and the 5 terms represent cumulative numbers of new detectable infections occurring over discrete 1-d periods. In Equation (3), each '/ term denotes a cumulative distribution function (CDF) of a log-normal distribution, and in Equation (4), each . term denotes a CDF of a gamma distribution. In other words, 5 is the cumulative number of individuals infected in the period ( 5 , 51# ) who will eventually be detected, and 045 is the probability that one of these individuals becomes symptomatic and is detected in the period ( 0 , 01# ), where 0 ≥ 5 .
The functional form of our curve-fitting model is defined by Equations (2)-(4), which are derived from Equation (1). As can be seen by inspecting Equation (1), the curve-fitting model has six parameters: , ! (which is hidden in the definition of 0 ), , , , and . As noted earlier, estimates are available for and from Lauer et al. (12) . These parameters characterize the variable duration of the incubation period, which starts at infection and ends at the onset of symptoms. Thus, we take and to have fixed region-independent values. We take the remaining parameters-(a population size/scaling factor), ! (the start time of the local epidemic), and and (the parameters that determine the shape of the infection curve ( ))-to have adjustable region-specific values. In our daily inferences, we consider one additional region-specific adjustable parameter, the dispersal parameter of the likelihood function (see Equation (31) below). The value of this parameter, , is inferred jointly with the values of , ! , , and .
The compartmental model, which is illustrated in Appendix Figure 1 , consists of the following 25 ordinary differential equations (ODEs): 
:% ",* :3 = (^0,6 + 0,7 _ + ' 04#,( − ' 0,( , for = 3, 4, 5 The initial condition is taken to be 6 ( ! ) = ! , 6 ( ! ) = ! = 1, with all other populations ( 7 , #,6 , … , <,6 , #,7 , … , <,7 , ;,( , … , <,( , 6 , 7 , ( , 7 , ( , , , and ) equal to 0. The parameter ! denotes the total region-specific population size. Thus, we assume that the entire population is susceptible at the start of the epidemic at time = ! >0, where time = 0 is 0000 hours on 21-January-2020. The parameter ! , which we always take to be 1, denotes the number of infectious symptomatic individuals at the start of the regional epidemic.
Subscripts attached to state variables are used to denote subpopulations. The subscripts and are attached to variables representing mixing and protected populations, respectively. For example, the variables 6 and 7 denote the population sizes of mixing and protected individuals who are susceptible to infection. Individuals in a protected population practice social distancing; individuals in a mixing population do not. The approach that we have taken to model social distancing is similar to that of Anderson et al. (30) .
The incubation period is divided into five stages. The numerical subscripts 1, 2, 3, 4, and 5 attached to variables indicate progression through these five stages. Exposed individuals in the incubation period, except for those in the first stage, are taken to be infectious. They are also taken to lack symptoms. They are either presymptomatic (i.e., individuals who will later develop symptoms) or asymptomatic (i.e., individuals who will never develop symptoms).
The subscript is attached to variables representing populations of quarantined individuals. The state variable ( is a special case; it accounts for symptomatic individuals who are quarantined as well as individuals who are self-isolating because of symptom awareness.
The parameter ( characterizes the rate at which infected individuals move into quarantine because of testing and contact tracing. The parameter ( characterizes the rate at which symptomatic individuals self-isolate because of symptom awareness. We recognize that susceptible individuals may enter quarantine (through contact tracing) but we assume that the size of the quarantined population is negligible compared to that of the total susceptible population and that susceptible individuals entering quarantine leave quarantine as susceptible individuals.
The parameters and $ < 1 characterize transmission of disease: characterizes the rate of transmission attributable to contacts between two mixing individuals, $ characterizes the rate of transmission attributable to contacts between one mixing and one protected individual, and $ ; characterizes the rate of transmission attributable to contacts between two protected individuals. Infectious individuals taken to contribute to COVID-19 transmission include those in the following pools: ;,6 , … , <,6 and ;,7 , … , <,7 , 6 and 7 , and 6 and 7 . Recall that we do not consider individuals in the first stage of the incubation period (i.e., individuals in # pools) to be infectious. The assumption is that these individuals are not shedding enough virus to be infectious (or detectable in surveillance testing).
The variables #,6 , … , <,6 and #,7 , … , <,7 denote the population sizes of mixing and protected exposed individuals in the five stages of the incubation period. The variables ;,( , … , <,( denote the population sizes of quarantined exposed individuals in the five stages. There is no #,( population, as we assume that individuals in the first stage of the incubation period are unlikely to test positive for SARS-CoV-2 or to be reached in contact tracing efforts before leaving the # state. The parameter ' characterizes disease progression, from one stage of the incubation period to the next and ultimately to an immune clearance phase. Individuals leaving the < pools enter the immune clearance phase, meaning that they become eligible for recovery. An individual leaving an < pool with symptom onset enters an I pool, whereas an individual leaving an < pool without symptom onset enters an A pool. Individuals in pools are considered to have mild disease (with the possibility to progress to severe disease).
The dynamics of social distancing are characterized by three step functions (i.e., piecewise constant functions having only finitely many pieces): 8 , Λ 2 , and 2 . The subscripts attached to these functions denote times: is a particular time, whereas is a set of times, as discussed later. The value of 8 switches from 0 to 1 at time = > ! , the start of an initial social-distancing period. As discussed later, the function Λ = defines a timescale for change in social-distancing practices for one or more distinct periods of social distancing, and the function 2 establishes a setpoint for the fraction of the total population of susceptible and infectious, nonquarantined/non-self-isolated/non-hospitalized individuals adhering to social-distancing practices for one or more distinct periods of social distancing.
The parameter & denotes the fraction of infected individuals who never develop symptoms (i.e., the fraction of all cases that are asymptomatic). The variables 6 and 7 denote the sizes of the populations of mixing and protected individuals who have been infected, progressed through the incubation period, are currently in the immune clearance phase, and will never develop symptoms. The parameter & characterizes the rate at which asymptomatic individuals recover.
The variable tracks recoveries of asymptomatic individuals, symptomatic individuals with mild disease, and symptomatic (hospitalized) individuals with severe disease. All individuals who recover are assumed to have immunity.
The variables 6 and 7 denote the sizes of the populations of mixing and protected symptomatic individuals with mild disease. The parameter + characterizes the rate at which symptomatic individuals with mild disease recover or progress to severe disease. The parameter ) is the fraction of symptomatic individuals who progress to severe disease requiring hospitalization. As a simplification, we assume that all individuals with severe disease are hospitalized. The variable represents the population size of hospitalized individuals, which are taken to be quarantined. Thus, the model does not consider nosocomial transmission.
The parameter * denotes the fraction of (hospitalized) individuals with severe disease who recover. The parameter ) characterizes the hospital discharge rate, i.e., the rate at which hospitalized individuals with severe disease either recover or die. The variable tracks deaths.
The time-dependent terms 6 ( , ) and 7 ( , ) appearing in Equations (5) where = ( % , & ), % is a constant characterizing the relative infectiousness of presymptomatic individuals compared to symptomatic individuals and & is a constant characterizing the relative infectiousness of asymptomatic individuals compared to symptomatic individuals. Recall that we assume that individuals in the first stage of the incubation period (i.e., individuals in either the #,6 or #,7 population) are not infectious. We also assume that the individuals in these populations cannot be quarantined until after transitioning to the ;,6 or ;,7 population (because they are assumed to test negative and because contact tracing is assumed to be too slow to catch individuals in the transient first stage of incubation). Recall that individuals in the 6 , individuals (i.e., individuals who will enter the populations) and asymptomatic individuals (i.e., individuals who will enter the populations).
The time-dependent terms 8 (t), 2 ( ), and Λ 2 ( ) appearing in Equations (5)-(10), Equations (13) and (14) and Equations (16) and (17) are step functions defined as follows:
where > ! is the time at which widespread social distancing initially begins, the integer ≥ 0 is the number of societal (major/widespread) shifts in social-distancing practices after the initial onset of social distancing, each 0 < 1 is a parameter characterizing the quasi-stationary fraction of susceptible individuals practicing social distancing during the ( + 1)th period of social distancing, each 0 is a constant defining a timescale for change in social-distancing practices during the ( + 1)th period of social distancing, = { ! , … , ,1# }, ! ≡ , ,1# ≡ ∞, and 01# > 0 for = 0, … , − 1. The value of 2 ( ) defines a setpoint for the quasi-stationary size of the protected population of susceptible individuals: 2 ( ) × 100% of the total susceptible population. The value of Λ 2 ( ) determines how quickly the setpoint is reached. As indicated in Equations (25) and (26), we only consider step-changes in the values of 2 ( ) and Λ 2 ( ), a simplification. Thus, for a period during which social-distancing practices are intensifying (relaxing), we increase (decrease) the value of 2 ( ) at the start of the period in a step-change and then hold it constant until the next step-change, if any.
To determine how consistent a particular parameterization of the compartmental model is with available COVID-19 surveillance data, we need to define a quantity-a model output-that corresponds to daily reports of the number of new confirmed COVID-19 cases. Case reporting by public health officials is typically daily. We expect that the vast majority of cases are detected because of symptom-driven (vs. random) testing and/or presentation in a clinical setting. Accordingly, as a simplification, we assume that individuals detected in surveillance are symptomatic. To define a model output comparable to the number of new cases reported on a given day, we start by considering the predicted cumulative number of presymptomatic individuals who become symptomatic while evading quarantine (because of contact tracing) until at least the onset of symptoms, which we will denote as @ . According to the model, the time rate of change of @ is given by the following equation:
The right-hand side of this equation gives the rate at which non-quarantined presymptomatic individuals exit the incubation period and enter the immune clearance phase, in which they are symptomatic and therefore taken to be detectable in local surveillance efforts.
Equation (27) and the ODEs of the compartmental model form a coupled system of equations, which can be numerically integrated to obtain trajectories for the state variables and @ , the expected cumulative number of symptomatic cases. From the trajectory for @ , we obtain a prediction for ( 0 , 01# ), the expected number of new COVID-19 cases reported on a given calendar date A , from the following equation:
where " is taken to be an adjustable region-specific parameter characterizing the time-averaged fraction of symptomatic cases detected. Equation (28) completes the formulation of our measurement model. ( 0 , 01# ) is the model output that we compare to 0 , the number of new cases reported on calendar date A .
The parameters of the compartmental model (Equations (5)- (26) ) and the auxiliary measurement model (Equations (27) and (28)) are taken to have either adjustable or fixed values. The adjustable parameter values are estimated (daily) through Bayesian inference on the basis of surveillance data (i.e., reports of newly detected cases). The fixed parameter values are held constant during inference; they are based on non-surveillance data and/or assumptions, which are discussed in the section below. In this section, we simply delineate the parameters with adjustable and fixed values. The compartmental model formulated for a given regional epidemic has a total of 16 + 3( + 1) parameters. The value of is structural; it sets the number socialdistancing periods considered.
The value of n corresponds to the number of periods of distinct social-distancing behaviors that follow an initial period of social distancing, which we take to begin at time = > ! . Here, we take = 0 or 1 for all regional epidemics of interest. Usually, we set = 0. In cases where we set = 1, this setting was motivated by second wave-type dynamics suggested by the surveillance data, which we take to indicate a relaxation of social-distancing practices at time = # > . The parameters of the initial social-distancing period are , ! , and ! . The parameters of the second social-distancing period, if considered, are # , # , and # . Thus, there are 3( + 1) social-distancing parameters, all of which are taken to be adjustable.
In addition to the 3( + 1) social-distancing parameters, we have 16 other parameters. Three of these define the initial condition: ! , 6 ( = ! ) = ! , and 6 ( = ! ) = ! , where ! is the time at which the epidemic begins, ! is taken to be the total population of the region of interest, and ! (the initial number of infected individuals) is always taken to be 1 (an assumption). We take ! to be adjustable and ! and ! to be fixed. The value of ! is set on the basis of population estimates by the US Census Bureau for the metropolitan statistical areas of interest (13) , which are delineated by the US Office of Management and Budget (10) .
There is only one more adjustable parameter of the compartmental model: , which characterizes the rate of disease transmission attributable to contacts among individuals within the mixing population. In the period before the onset of social distancing, from ! to , when 6 / ! ≈ 1, the instantaneous rate of disease transmission is 6 ( , ), where 6 ( , ) is the effective number of infectious individuals at time , a weighted sum of the numbers of symptomatic, presymptomatic, and asymptomatic individuals determined by = ( % , & ). We assume that exposed individuals after the first stage of disease incubation are infectious, as are asymptomatic individuals in the immune clearance phase who have passed through all 5 stages of disease incubation and who will never develop symptoms.
The remaining 12 parameters of the compartmental model, which are taken to have fixed, region-independent values, are as follows: $ , % , & , ' , ( , ( , & , ) , & , + , * and ) . Our estimates for these parameters are discussed in the section immediately below. It should be noted that settings for * and ) do not affect predictions of new cases because these parameters characterize recovery/morbidity of hospitalized individuals. The parameter * is the fraction of hospitalized individuals who recover, and the parameter ) characterizes the hospital discharge rate. Although nosocomial disease transmission is a significant concern, we assume that hospitalized individuals are effectively quarantined such that the overall rate of disease transmission in a given region is insensitive to the number of hospitalized individuals in that region.
Here, we summarize the rationale/justification for each of our estimates for the values of the following 12 parameters of the compartmental model: $ , % , & , ' , ( , ( , & , ) , & , + , * , and ) . These estimates are assumed to apply to all regions, i.e., we take these parameters to have region-independent values.
The parameter $ characterizes the effects of social distancing on disease transmission. Without social distancing, all contacts responsible for disease transmission are between mixing individuals (i.e., between individuals in the 6 and 6 pools) and the rate of transmission is characterized by . With social distancing, there are contacts involving 1 individual in a mixing population and 1 individual in a protected population (e.g., between individuals in the 6 and 7 pools or in the 6 and 7 pools) and also contacts involving 2 individuals in protected populations (e.g., between individuals in the 7 and 7 pools). In the model, the rates of transmission associated with these types of contacts are characterized by $ and $ ; , respectively. We can be confident that social distancing is protective (i.e., $ < 1) but there is little information available to suggest the magnitude of the effect. We arbitrarily set $ = 0.1, which can be interpreted to mean that a susceptible individual practicing social distancing has a 10-fold smaller chance of becoming infected than a susceptible individual that is not practicing social distancing. In exploratory analyses, wherein we allowed $ to be a free parameter, we found that its inferred value is positively correlated with the extent of social distancing, which is determined by the relevant social-distancing setpoint parameter (viz., ! ). Thus, we interpret the inferred quasi-stationary value of 7 to be an effective population size. If our estimate for $ is too high (i.e., we underestimate the protective effect of social distancing), the effective size will be larger than the true size. Conversely, if our estimate for $ is too low, the effective size will be smaller than the true size.
The parameters % and & characterize the relative infectiousness of individuals without symptoms during the incubation period and the immune clearance phase, respectively. Infectiousness is compared to that of a symptomatic individual. Using a one-step real-time reverse transcriptase-polymerase chain reaction (rRT-PCR) assay to quantify viral RNA abundance in nasopharyngeal and oropharyngeal samples, Arons et al. (14) determined rRT-PCR threshold cycle (Ct) values for 17 symptomatic and 24 presymptomatic individuals. The former group of individuals had typical symptoms, and the latter group of individuals lacked symptoms at the time of testing but later developed symptoms (within 1 week after testing). At the time of testing, the median Ct values for symptomatic and presymptomatic individuals were 24.8 and 23.1, respectively. (NB: Ct value is inversely proportional to abundance.) Based on these results and an assumption that infectiousness is proportional to viral load, we estimate that % = 1.1. An estimate for % greater than 1 is consistent with the findings of He et al. (25) , who inferred that viral load is maximal 0.7 d before the onset of symptoms from an analysis of temporal viral load data and information available about infector-infectee transmission pairs. Over a period of 19 d, Nguyen et al. (15) performed daily rRT-PCR assays for viral RNA in nasopharyngeal samples from 17 symptomatic and 13 asymptomatic individuals. Ngyuen et al. (15) developed a curvefitting model for each group to characterize their viral decay kinetics. These models indicate that the mean Ct value for symptomatic individuals was roughly 90% of the mean Ct value for asymptomatic individuals over the first week of the study, after which most individuals tested negative or had a Ct value near the threshold of detection, 40. Thus, we estimate that & = 0.9.
The parameter ' characterizes the duration of the incubation period. In the model, the incubation period is divided into 5 stages (for reasons explained shortly). The waiting time for completion of all 5 stages is described by an Erlang distribution with a shape parameter = 5 and a scale parameter = 1/ ' . Lauer et al. (12) estimated times of exposure and symptom onset for 181 confirmed cases and found that the median time between SARS-CoV-2 infection and onset of COVID-19 symptoms is 5.1 d. Lauer et al. (12) also found that the empirical distribution of waiting times is fit by an Erlang distribution with = 6 and = 0.88 d. This latter finding suggests that the empirical waiting time distribution can be reproduced by dividing the incubation period into 6 stages and setting ' = 1.14 d 4# . However, an Erlang distribution with = 5 and = 1.06 d has a nearly identical shape. Because simulation costs are reduced by dividing the incubation period into 5 instead of 6 stages, we considered 5 stages in the model. The distribution of waiting times estimated by Lauer et al. (12) is reproduced by our model when we set ' = 0.94 d 4# .
The parameters ( and ( characterize testing-driven quarantine and symptom-driven selfisolation. We assume that testing is random. Thus, the number of infected individuals moving into quarantine per d is the number of infected individuals subject to quarantine times the fraction of the total population tested per d times a multiplier capturing the effect of contact tracing. We take the multiplier to be average household size, 2.5 (US Census Bureau). Thus, based on approximately 500,000 tests per d in the US (https://covidtracking.com/data/us-daily) and a total population of 330 million (https://www.census.gov/popclock/), we estimate ( = 0.0038 d 4# . We assume ( = 0.4 d 4# . With this setting, the median waiting time from onset of symptoms to initiation of self-isolation is approximately 40 h. A faster timescale for selfisolation is probably not realistic despite general awareness of the COVID-19 pandemic, because as considered in the study of Böhmer et al. (26) , for any given individual, there may be a prodromal phase of ~1 d marked by non-COVID-19-specific symptoms other than fever and cough. The parameter ) is the fraction of symptomatic individuals progressing to severe disease. We set ) such that our model predicts a uniform infection fatality rate (IFR) consistent with that determined by Perez-Saez et al. (18) from serological survey results and death incidence reports: 0.0064 (~0.64%). For a discussion of other IFR estimates, which tend to be similar, see Grewelle and De Leo (29) . According to our model, IFR is given by (1 − & ) ) (1 − * ). This quantity is the fraction of all infected individuals predicted to develop symptoms and then to progress to severe disease (and hospitalization) and finally a fatal outcome. Thus, based on our estimates for & (0.44) and * (0.79) and the empirical IFR (0.0064), we set ) = !.!!D?
!.<D ×!.;# ≈ 0.054.
The parameter & characterizes the duration of infectiousness of asymptomatic individuals in the immune response phase. For each of 89 asymptomatic individuals, Sakurai et al. (17) reported the time between the first positive PCR test for SARS-CoV-2 and the first of two serial negative PCR tests. The mean duration of this period was ~9.1 d. We assume that this period coincides with the period of infectiousness and that this period encompasses both the incubation period and the immune response phase. According to our model, the mean duration of the incubation period is 5/ ' for both presymptomatic and asymptomatic individuals. Based on our earlier estimate that ' = 0.94 d 4# , the mean duration of the incubation period is estimated as 5.3 d. Accordingly, the mean duration of the immune clearance phase for asymptomatic individuals is estimated as 9. If ) ≪ 1, the parameter + characterizes the duration of infectiousness of individuals who develop mild COVID-19 symptoms (i.e., symptoms not severe enough to require hospitalization). Wölfel et al. (20) attempted to isolate live virus from clinical throat swab and sputum samples collected from 9 patients at multiple time points after the onset of mild COVID-19 symptoms. Roughly 67%, 38%, and 0% of attempts to isolate virus were successful at 6, 8, and 10 d after infection, respectively. Assuming that a negative culture coincides with loss of infectiousness, we estimate that + = −
The parameters * and ) characterize the hospital stays of the severely ill. These parameters affect predictions of COVID-19-caused deaths and hospital resource utilization but do not affect the predicted transmission dynamics, because we assume that hospitalized patients are effectively quarantined and do not contribute significantly to disease transmission, i.e., there is no ) term in 6 or 7 (see Equations (22) and (23)). The parameter * is the fraction of hospitalized patients who recover, and the parameter ) characterizes the rate at which patients are discharged (as either recovered or dead). Richardson et al. (19) reported that the overall median length of hospital stay for 2,634 discharged patients (alive or dead) was 4. 
We assume that the likelihood of a set of adjustable parameter values M given a report of 0 new cases on calendar date 0 , which we will denote as ℒ 0 ( M ; 0 ), is given by the following equation:
where 0 is a non-negative integer (the number of new cases reported), is an integer indicating the date 0 or the period ( 0 , 01# ); nbinom( 0 ; , 0 ) is the probability mass function of the negative binomial distribution NB( , 0 ), which has two parameters, 0 ∈ [0,1] and > 0; and M is a model-dependent ordered set of feasible (i.e., allowable) values for the adjustable model parameters (e.g., , ! , , and in the case of the curve-fitting model) augmented with a feasible value for . Recall that 0 ≡ ! + d, where ! > 0 is 0000 hours of ! , the start date of the local epidemic. We take the dispersion parameter of NB( , 0 ) to be date/time-independent and infer the value of jointly with the values of the model parameters. Because the deterministic models are designed to capture the mean behavior of the random reporting numbers, we have to impose the constraint ( 0 , 01# ) = [NB( , 0 )] = (1 − 0 )/ 0 , which leads to the date/time-dependent parameter 0 of the negative binomial distribution:
If + 1 daily case reports are available, from date ! to date Q , we assume that each likelihood ℒ 0 ( M ; 0 ) given by Equations (29) and (30) 
We chose the Bayesian inference framework to parametrize the models with uncertainty quantification. In Bayesian inference, given a set of data , the probability of each set of the parameters, denoted in M is constrained by the Bayes formula 
Here, the ℙ{ M } is the prior parameter distribution, which represents our belief of how the model parameters should distribute in the parameter space Ω, and ℙ{ | M } is the likelihood of the parameter set M given the dataset , that is, ℒ ( In contrast of many modeling analyses which focus on identifying the parameter distributions, we are interested in projections of the models, whose parameters are inferred by past data, into the future. To this end, we evaluate the model with a probabilistic parameter set distributed by the obtained posterior distribution ℙ{ M | }. Formally, we denote the prediction of the confirmed cases between future day 0 and 01# by our deterministic model with a set of parameters M by ( 0 , 01# ; M ). Recall that this deterministic prediction represents the mean of the fundamentally random new confirmed cases reported in a future interval ( 0 , 01# ). If there was only parametric uncertainty which propagates through the deterministic model, the confirmed cases reported in a future interval ( 0 , 01# ) would be distributed according to
However, there is also observation noise, which we model by a negative binomial distribution. The observation noise also needs to be injected into the prediction to quantify the full uncertainty. The full prediction accounting for parametric uncertainty is a random variable distributed according to 
In practice, the above random variable is resampled from the posterior chain derived from the MCMC sampling. We denote the MCMC posterior chain by q M (#) , M (;) … M (/) Ž. We sample the posterior chain and denote the resampled parameter set by M S and the deterministic prediction of that resampled parameter in interval ( 0 , 01# ) by ( 0 , 01# ; M S ). Then, we generate a negative binomial random number with the first parameter of the negative binomial distribution set as S /( S + ( 01# , 0 ; M S )) where S is the resampled which is also a free parameter in M S and is inferred in the MCMC. We repeat the resampling procedures and use the generated samples to compute the percentile of the past history and future prediction.
Our aim is to perform the Bayesian inference daily as soon as a new regional confirmed case number is reported. Although the Bayesian framework allows a sequential analysis, that the previous derived posterior distribution is used as a prior and the new inference problem involves only one new data point, in practice, such an analysis is difficult if the posterior distribution cannot be emulated or interpolated from the discrete posterior chain. Our analysis shows that in some regions, the posterior is far from Gaussian, making accurate interpolation or emulation difficult. Thus, we are not adopting this workflow, and instead, we perform the inference with all the data points collected up to the time of inference. Nevertheless, we warm-start the MCMC chain from the maximum a posteriori estimator estimated from the previous derived chain, and we also use the previously derived chain for estimating the optimal covariance matrix for the proposal of the normal symmetric random-walk Metropolis sampler. This approach allows an online learning of the optimal proposal which significantly reduces the mixing time.
Because the variability of the data due to the regional and temporal differences, it is difficult to identify a universal sampling strategy (the proposal kernel). Thus, we adopted an adaptive Metropolis algorithm, specifically Algorithm 4 in Andrieu and Thoms (21) to accommodate the regional and temporal differences. For all the model parameters, we assume their priors are uniformly distributed, whether they are proper or improper. Denote the 0:00 of the calendar date of the first confirmed case in a specific region by TAUVW and the total population of that region by ! . For the curve-fitting model, we assume the parameters are bounded by ∈ (0, ! ), ! ∈ ( TAUV3 − 21, TAUV3 ), ∈ (0, ∞), ∈ (0, ∞), ∈ (0, ∞). For the compartmental model, we assume that the parameters are bounded by ! ∈ (0, ∞), ∈ ( ! , ∞), ∈ (0, ∞), ∈ (0,10), 7 ∈ (0,1), " ∈ (0,1). We assume ∈ (0, ∞) for both models. We adopted rejection-based sampling to assure the parameter values are sampled in the hypercube.
We start inference with an isotropic proposal kernel, that randomly perturbs the parameter values by independent Gaussian proposals whose standard deviations are set to be 5% of the parameter values. We carry out the standard Metropolis-Hastings algorithm for 5 × 10 ? iterations first to identify local minimum. Then, we turned on the adaptive Metropolis algorithm to calculate the covariance matrix on-the-fly, for another 5 × 10 ? iteration, when we turned on an on-the-fly learning for optimal proposed increment, i.e., in Algorithm 4 of Andrieu and Thoms (21) . Because the weight for learning the proposed increment decays 1/iteration (21) , the proposed increment stabilizes after about 10 > more iterations. We began to collect the statistics from 1.5 × 10 < iteration, until the simulation finishes at 6 × 10 < iterations.
Except for the first time of the procedure (i.e., online learning and day-to-day operation), we warm start the simulation from the previously obtained best-fit (maximum a posteriori (MAP) estimator) and with the previously obtained covariance matrix and proposed increment. We carry out standard Metropolis-Hastings algorithms for 2.5 × 10 ? iterations first to identify a local minimum, noting that it is often not far away from the previously identified MAP. We then turn on the adaptive MCMC algorithm to calculate the covariance matrix on-the-fly, again for another 5 × 10 ? iteration. We then use another 2.5 × 10 ? iterations to calculate the optimal proposed