dataesr/scientific-paragraphs-categorization huggingface dataset:
ACCELERATING MCMC BY BREAKING THE PROBLEM INTO PIECES The explosion in the collection and analysis of "big" datasets in recent years has brought new challenges to the MCMC algorithms that are used for Bayesian inference. When examining whether or not a new proposed sample is accepted at the accept-reject step, an MCMC algorithm such as the Metropolis-Hastings version needs to sweep over the whole data set, at each and every iteration, for the evaluation of the likelihood function. MCMC algorithms are then difficult to scale up, which strongly hinders their application in big data settings. In some cases, the datasets may be too large to fit on a single machine. It may also be that confidentiality measures impose different databases to stand on separate networks, with the possible added burden of encrypted data  (Aslett et al.
Sub-posterior methods are motivated by the independent product equation: (1) π(θ) ∝ k i=1   π 0 (θ) 1/k ∈X i p(x |θ)   = k i=1 π i (θ) , and they target the densities π i (θ) (up to a constant) in their respective MCMC steps. They thus bypass communication costs  (Scott et al. , 2016) , by running MCMC samplers independently on each batch, and they most often increase MCMC mixing rates (in effective samples sizes produced by second), given that the sub-posterior distributions π i (θ) are based on smaller datasets. For instance,  Scott et al. (2016)  combine the samples from the sub-posteriors, π i (θ), by a Gaussian reweighting. Neiswanger et al. (2013)  estimate the sub-posteriors π i (θ) by non-parametric and semi-parametric methods, and they run additional MCMC samplers on the product of these estimators towards approximating the true posterior π(θ).
ACKNOWLEDGEMENTS Christian P. Robert  is grateful to  Gareth Roberts ,  Mike Betancourt , and  Julien Stoehr  for helpful discussions. He is currently supported by an  Institut Universitaire de France   2016-2021  senior grant. Changye Wu  is currently a Ph. D. candidate at  Université Paris-Dauphine  and supported by a grant of the  Chinese Government (CSC). Víctor Elvira acknowledges support from the  Agence Nationale de la Recherche of France  under  PISCES  project ( ANR-17-CE40-0031-01 ), the  Fulbright program , and the  Marie Curie Fellowship  (FP7/2007-2013) under REA grant agreement n. PCOFUND-GA-2013-609102, through the  PRESTIGE program. The authors are quite grateful to a reviewer for his or her detailed coverage of an earlier version of the paper, which contributed to significant improvements in the presentation and coverage of the topic.
Introduction: The global epidemiology of many infectious diseases is changing, but little attention has been paid to whether the timing of seasonal influenza epidemics changed in recent years. This study investigated whether the timing of the peak of influenza epidemics has changed in countries of the World Health Organization (WHO) European Region between 1996 and 2016. Methods: Surveillance data were obtained from the WHO FluNet database. For each country and season (July to June of the next year), the peak was defined as the week with the highest 3-week moving average for reported cases. Linear regression models were used to test for temporal trends in the timing of the epidemic peak in each country and to determine whether this differed geographically. Results: More than 600,000 influenza cases were included from 38 countries of the WHO European Region.
Data Influenza virological surveillance data were obtained from the publicly available web-based database FluNet, which is coordinated by the WHO  [18]. Information on the weekly number of laboratory-confirmed cases of influenza (overall and by virus type, subtype and lineage) is entered into the FluNet database by the national influenza centres and other influenza reference laboratories of 113 countries participating in the Global Influenza Surveillance and Response System. On 6 November 2016, influenza surveillance data were downloaded for the period between week 27/1996 (starting 1 July 1996) and week 26/2016 (starting 27 June 2016), henceforth referred to as the study period, for all 53 countries in the WHO European Region  [19]. We excluded from the dataset the 2009/10 influenza season and influenza A(H1N1)pdm2009 influenza cases reported between April and June 2009 (for all countries), seasons with fewer than 20 weeks of reporting or fewer than 100 influenza cases overall (only for countries where this applied), and countries that had data for fewer than five influenza seasons.
Sensitivity analysis The impact of very early or late epidemic peaks was assessed by excluding one season at a time from the analyses. The impact of geographical outliers was assessed by first identifying outliers and influential points (i. e. points that significantly influence the output of Model 2) using studentised residuals, Cook's D and difference in fits (DFFITS)  [23, 24] ) and then by excluding the identified countries from analyses. Finally, to assess the impact of removing countries  with limited data, all analyses were repeated excluding countries that had data for less than seven or less than 10 seasons (the basic analysis already excluded countries with less than five seasons, see above). In order to test the robustness of results, all analyses were repeated using the 3-week moving average of the positive detection rate (defined as the ratio of the number of reported cases over the number of respiratory specimens processed in the same week) instead of the 3-week moving average of the number of reported cases.
Acknowledgements We would like to thank  EU/EEA  member states' laboratories reporting their influenza surveillance data to The European Surveillance System (TESSy) database, hosted at the  European Centre for Disease Prevention and Control (ECDC) , from which data are shared with the global WHO FluNet database. Dr Phillip Leventhal  (4Clinics, Paris, France) provided scientific editing, which was paid for by  Sanofi Pasteur. Source of funding: The study is supported by a research grant by  Sanofi Pasteur : the ' Global Epidemiology of Influenza B' research project.
3 N 8. 0 W 10. 8 17 303 8,956 Romania 46. 0 N 25. 0 E 21. 7 14 335 5,729 Russian Federation 60. 0 N 100. 0 E 142. 4 11 4,629 79,376 Serbia 44. 0 N 21. 0 E 7. 2 6 334 1,809 Slovakia 48. 4 N 19. 3 E 5. 4 5 279 1,657 Slovenia 46. 1 N 14. 5 E 2 13 425 6,924 Spain 40. 0 N 4. 0 W 48. 1 11 4,697 41,185 Sweden 62. 0 N 15. 0 E 9. 8 15 2,452 57,236 Switzerland 47. 0 N 8. 0 E 8. 1 19 361 20,771 Turkey 39. 0 N 35. 0 E 79. 4 7 1,056 14,725 Ukraine 49. 0 N 32. 0 E 44. 4 6 620 6,586 United Kingdom 54. 0 N 2. 0 W 64. 1 17 1,079 61,507 [35]ial dataset included 819,332 influenza cases from 676 seasons in 49 countries of the WHO European Region; no data were available for Andorra, Cyprus, Monaco and San Marino. After applying the exclusion E: east; N: north; S: south; W: west; WHO: World Health Organization.
Sourceb Most recent estimates[35]. criteria, the final dataset used for the analysis included 635,046 influenza cases from 400 seasons in 38 countries (Table1). The median number of seasons per country was 10 and ranged from five (Belarus, Bulgaria, Lithuania, Moldova and Slovakia) to 19 (Finland and Switzerland). The median number of influenza cases reported per season was 601 (interquartile range: 267-1,414).
0 N 8. 0 E 19 2. 12 0. 26 to 3. 98 0. 028 25. 4% Turkey 39. 0 N 35. 0 E 7 -1. 85 -15. 66 to 11. 96 0. 745 2. 3% Ukraine 49. 0 N 32. 0 E 6 -9. 63 -17. 68 to -1. 58 0. 029 73. 4% United Kingdom 54. 0 N 2. 0 W 17 2. 26 -0. 80 to 5. 32 0. 136 14. 2% CI: confidence interval; E: east; N: north; S: south; W: west; WHO: World Health Organization. Source: WHO FluNet database[18]for July 1996 to June 2016. Data from July 2009 to June 2010 were excluded.
WHO FluNet database[18]for July 1996 to June 2016. Data from July 2009 to June 2010 were excluded. a Latitude and longitude are for the centroid or (if not available) the centre point of each country. b A beta coefficient above (below) zero means that the temporal shift of the influenza peak was greater (smaller) by moving from east to west. Dependent variable Temporal shift of the influenza epidemic peak (days/season) Independent variable Country longitude a Model output beta b 95% CI p value R 2 All countries included (n = 38) 0. 077 0. 034 to 0. 121 0. 001 26. 6% One season removed at a time (range) 0. 056 0. 091 0. 018 to 0. 095 0. 037 to 0. 145 0. 005 0. 002 24. 0% 27. 0% Countries with ≥7 seasons of data (n = 27) 0. 062 0. 039 to 0. 085 < 0. 001 53. 5% Countries with ≥10 seasons of data (n = 20) 0.
1) can be used. The benchmark can be constructed for various decision criteria and then be adapted to evaluate any design under investigation. Wang and Ivanova (2015)  is to find the dose with the mean response closest to the target response γ. The criterion of choosing the dose which maximises the probability of the average level of toxicity µ j to be in the ε neighbourhood of γ was considered. Let g j (•|y j ) be a probability density function of µ j given the data y j. Then, the decision criterion takes the form T (y j ) = γ+ε γ-ε g j (v|y j )dv. (2. 3) The TD is the dose for which the criterion T (y j ) is maximised. Following the original framework, γ = 0. 1 and ε = 0. 01 are chosen. Using the complete information generated in Table  1  and the density function of Normal distribution with corresponding mean and variance parameters yields: T (y 1 ) = 0.
, n, j = 1,. , m and store y j = (y 1j ,. , y nj ). 4. Compute T (y j ) for all j = 1,. , m, find dose J for which T (y J ) is maximised (minimised) and set Z s = J. 5. Repeat steps 2-4 for s = 1,. , S simulated trials 6. Use Z(j) = S s=1 I (Z s = j) /S as the selection proportion of dose d j , j = 1,. , m The proposed benchmark can be applied to a wide range of distributions as it requires the quantile information only, which is available for many distributions in various statis-tical software (for example, qbinom, qnorm, qexp , etc in R (R Core Team, 2015)). Note that the probability integral transform can be also applied to discrete random variables in which case the quantile transformation F -1 j (•) is given explicitly. It is easy to see that using the F -1 j (•) corresponding to a Bernoulli random variable in Algorithm 1 results in the original benchmark construction proposed by O'Quigley and others  (2002).
j (•|y (2) j ) are probability density functions of an efficacy response and of a toxicity probability given the data y (1) j , y (2) j , respectively, and θ (1) , θ (2) are controlling probabilities. This decision criterion is used to construct the benchmark in this setting. Applied to the benchmark, the integrals in (2. 5) are computed using density functions of Beta distribution and Normal distribution for toxicity and efficacy outcomes, respectively. Using summary statistics given in Table  2  and controlling probabilities θ (1) = θ (2) = 0. 50, values of the criterion are T (y (1) 1 , y (2) 1 ) = 0. 30; T (y (1) 2 , y (2) 2 ) = 0. 79; T (y  (1)   3 , y   (2)  3 ) = 1. 28; T (y (1) 4 , y (2) 4 ) = 0;. The criterion is maximised for dose level d 3 which is selected as the TD in this single trial.
Acknowledgments This project has received funding from the  European Union 's  Horizon 2020 research and innovation programme  under the  Marie Sklodowska-Curie  grant agreement No  633567. Xavier Paoletti  is partially funded by the  Institut National du Cancer (French NCI ) grant SHS-2015 Optidose immuno project. This report is independent research arising in part from  Prof Jaki 's  Senior Research Fellowship  ( NIHR-SRF-2015-08-001 ) supported by the  National Institute for Health Research. The views expressed in this publication are those of the authors and not necessarily those of the  NHS , the  National Institute for Health Research  or the Department of Health.
This paper introduces a framework for speeding up Bayesian inference conducted in presence of large datasets. We design a Markov chain whose transition kernel uses an unknown fraction of fixed size of the available data that is randomly refreshed throughout the algorithm. Inspired by the Approximate Bayesian Computation (ABC) literature, the subsampling process is guided by the fidelity to the observed data, as measured by summary statistics. The resulting algorithm, Informed Sub-Sampling MCMC (ISS-MCMC), is a generic and flexible approach which, contrary to existing scalable methodologies, preserves the simplicity of the Metropolis-Hastings algorithm. Even though exactness is lost, i. e the chain distribution approximates the posterior, we study and quantify theoretically this bias and show on a diverse set of examples that it yields excellent performances when the computational budget is limited.
Introduction The development of statistical methodology that scale to large datasets represents a significant research frontier in modern statistics. This paper presents a generic and flexible approach to directly address this challenge when a Bayesian strategy is followed. Given a set of observed data (Y 1 ,. , Y N ), a specified prior distribution p and a likelihood function f , estimating parameters θ ∈ Θ of the model proceeds via exploration of the posterior distribution π defined on (Θ, B(Θ)) by π(dθ | Y 1 ,. , Y N ) ∝ f (Y 1 ,. , Y N | θ)p(dθ). (1) In this paper, we propose Informed Sub-Sampling MCMC (ISS-MCMC), a novel methodology which aims to make the best use of a computational resource available for a given computational run-time, while still preserving the celebrated simplicity of the standard M-H sampler.
We show that setting the similarity measure as the squared L2 distance between the full dataset and subsample maximum likelihood estimators is supported by theoretical arguments. Section 2 presents a striking real data example which we hope will help the reader to understand the problem we address and motivate the solution we propose, without going into further technical details at this stage. In Section 3, we provide theoretical results concerning exponential-family models, which we illustrate through a probit example. This section allows us to justify our motivations supporting the Informed Sub-Sampling general methodology which is rigorously presented in Section 4. In Section 5, we study the transition kernel of our algorithm and show that it yields a Markov chain targeting, marginally, an approximation of π.
An introductory example We showcase the principles of our approach on a first real data example. The problem at hand is to infer some template shapes of handwritten digits from the MNIST database (http://yann. lecun. com/exdb/mnist/). Example 1. The data Y 1 , Y 2 ,. are modelled by a deformable template model  (Allassonnière et al. , 2007). Each data Y i is a 15 × 15 pixel image representing an handwritten digit whose conditional distribution given its class J(i) ∈ (0, 1,. , 9) is a random deformation of the template shape, parameterized by a d = 256 dimensional vector θ J(i). Assuming small deformations, the model is similar to a standard regression problem: Y i = φ(θ J(i) ) + σ 2 i , (2) where Y i is regarded as a vector R 225 , φ : R 256 → R 225 is some deterministic mapping and σ > 0 is the standard deviation of the additive noise i ∼ N (0 225 , Id 225 ).
Optimal subsets for the Kullback-Leibler divergence between π and πn Recall that for two measures π and π defined on the same measurable space (Θ, B(Θ)), the Kullback-Leibler (KL) divergence between π and π is defined as: KL(π, π) = E π log π(θ) π(θ). (6) Although not a proper distance between probability measures defined on the same space, KL(π, π) is used as a similarity criterion between π and π. It can be interpreted in information theory as a measure of the information lost when π is used to approximate π, which is our primary concern here. We now state the main result of this section. Proposition 1. For any subset U ∈ U n , define the vector of difference of sufficient statistics between the whole dataset and the subset Y U as ∆ n (U ) = N k=1 S(Y k ) -(N/n) k∈U S(Y k ). (7) Then, the following inequality holds: KL {π, πn (U )} ≤ B(Y, U ) , (8) where B(Y, U ) = log E π exp { E π (g(θ)) -g(θ) ∆ n (U ) } (9) and • is the L2 norm.
Weighting the subsamples Consider the distribution ν n, on the discrete space U n defined for all ≥ 0 by: ν n, (U ) ∝ exp -∆ n (U ) 2 , for all U ∈ U n. ( 12 ) ν n, assigns a weight to any subset according to their representativeness with respect to the full dataset. When = 0, ν n, is uniform on U n while when → ∞, ν n, is uniform on the set of subset(s) that minimize(s) U → ∆ n (U ). Proposition 2 suggests that for exponential models, the optimal inference based on subsamples of size n is obtained by picking the subposterior π n (U ) (5) using the distribution U ∼ ν n, with → ∞. We now state Proposition 3. This result is important even though somewhat obscure at this stage. Indeed, we will show that it is a necessary condition for the method we introduce in Section 4 to converge. In fact, moving away to general models (i.
The pdf of the posterior distribution π and any distribution πn (U ) ∈ F n writes respectively as π(θ | Y 1:N ) ∝ p(θ) (1 -α(θ)) N α(θ) 1 -α(θ) N k=1 Y k , πn (θ | Y U ) ∝ p(θ) (1 -α(θ)) N α(θ) 1 -α(θ) (N/n) k∈U Y k , where p is a prior density on θ. Again, in this example, the posterior density is easy to evaluate pointwise, even when N is extremely large, as it only requires to sum over all the binary variables Y 1 ,. , Y N. As a consequence, samples from π can routinely be obtained by a standard M-H algorithm and similarly for any distribution πn (U ) ∈ F n. We simulated N = 10, 000 simulated data Y 1 ,. , Y N from (  14 ), with true parameter θ * = 1. We used the prior distribution p = N (0, 10). In this probit model, S is simply the identity function, implying that ∆ n (U ) gives the absolute value of the difference between the scaled proportion of 1 and 0's between the full dataset and the subset Y U.
M of E{f (Y | θ)/f (Y U | θ) N/n } for θ ∈ (0, 1. 5), for n = 1, 000. When = 0, νn, is uniform (i. e. an identical weight is assigned to all the subsamples) and as a consequence M ≡ ∞. Conversely, when 0, the mass of νn, spreads over the best subsamples Y U , U ∈ Un (i. e. those minimizing ∆n(U )) and the bound M is smaller than e 2. Indeed, by assigning a weight νn, (U ) ∝ exp{-∆n(U ) 2 } those subsamples Y U , U ∈ Un that have a large ∆n(U ) will yield a negligible contribution to the expectation, hence preventing from divergence. Figure  2  reports the density functions of π and several other distributions πn (U ) ∈ F n , for n = 100 and n = 1, 000, with different values for the quantity ∆ n (U ) (7). This plot, as well as the quantitative result of Table  1  are consistent with the statement of Corollary 1: when learning from a subsample of n data, one should work with a subset U featuring a perfect match with the full dataset, i.
g. Nunes and Balding, 2010; Csilléry et al. , 2010; Marin et al. , 2012; Fearnhead and Prangle, 2012). More details on validation of summary statistics are discussed in Section 5. 3. 2. Because the statistics used to assess the representativeness of a subsample Y U w. r. t. the full dataset Y are only summary and not sufficient, the results of Section 3 are no longer valid. In particular, should an optimal subset U * minimising a distance between S(Y U ) and S(Y ) exist, inferring π through the approximation πn (U * ) is in no sense optimal. In fact, as shown in several examples of Section 6, this approximation is usually poor. In such a setting, it is reasonable to consider extending the set of subsamples of interest to a pool of good subsamples. This naturally suggests using the distribution ν n, (12) to discriminate between the subsamples, replacing sufficient by summary statistics and relaxing the assumption → ∞, in order to account for a collection of good subsamples.
Informed Sub-Sampling MCMC To avoid any confusion, we denote by { θi , i ∈ N} the sequence of parameters generated by the Informed Sub-Sampling Markov chain, by contrast to the Markov chain {θ i , i ∈ N} produced by the Metropolis-Hastings algorithm (Alg. 1). The pool of good subsamples used in the Informed Sub-Sampling inference is treated as a sequence of missing data U 1 , U 2 ,. and is thus simulated by our algorithm. More precisely, ISS-MCMC produces a Markov chain {( θi , U i ), i ∈ N} on the extended space Θ × U n. Inspired by the analysis of Section 3, the sequence of subsamples {U i , i ∈ N} is randomly updated in a way that favours those subsets whose summary statistics vector is close to that of the full dataset. Let R be a symmetric transition kernel on (U n , U n ), a transition ( θi , U i ) → ( θi+1 , U i+1 ) consists in the two following steps: i-(a) propose a new subset variable U ∼ R(U i , • ) (b) set U i+1 = U with probability β(U i , U ) = 1 ∧ b(U i , U ) , b(U i , U ) = exp ∆ n (U i ) 2 -∆ n (U ) 2 (19) and U i+1 = U i with probability 1 -β(U i , U ).
A4. Summary Statistics There exists a constant γ n < ∞, such that for all (θ, U ) ∈ Θ × U n |log f (Y | θ) -(N/n) log f (Y U | θ)| ≤ γ n N S(Y ) -S(Y U ). ( 26 ) Assumption A. 4 imposes a condition simultaneously on the model f and the summary statistics S. In particular, it assumes that for any θ ∈ Θ, the variation of the scaled likelihood of the subsamples Y U (U ∈ U n ) around f (Y | θ) is controlled by the distance between the full dataset Y and the subsample Y U , as measured through their summary statistics. This is a strong assumption which is unlikely to hold if Θ is not a compact set. It implies that even in absence of sufficient statistics, a result similar to Proposition 3 exists. One can also note that when n → N , the constant γ n in Eq. (  26 ) goes to zero.
In situations where the maximum likelihood estimator θ * (Y 1:n ) is easy and quick to evaluate numerically, we recommend setting S(Y 1:n ) = θ * (Y 1:n ). In the case of independent observations of a well-specified model, setting the summary statistics as the maximum likelihood estimate is justified by the following Proposition which implies that Assumption A. 4 holds, asymptotically, up to a constant. Proposition 7. We assume that the whole dataset comprises N = ρn independent observations and there exists some θ 0 ∈ Θ such that Y i ∼ f ( • | θ 0 ). Let θ * be the MLE of Y 1 ,. , Y N and θ * U be the MLE of the subsample Y U (U ∈ U n ). Then, there exists a constant β, a metric • θ0 on Θ and a non-decreasing subsequence {σ n } n∈N , (σ n ∈ N) such that for all U ⊂ {1, 2,. , ρσ n } with |U | = σ n , we have for p-almost all θ in an neighborhood of θ 0 : log f (Y 1:ρσn | θ) -ρ log f (Y U | θ) ≤ H n (Y, θ) + β + ρn 2 θ * U -θ * θ0 , (34) where plim n→∞ H n (Y, θ) P θ 0 = 0.
The proof is detailed in Appendix B and follows from a careful application of a Bernstein-von Mises theorem. Note that an extension of Proposition 7 to cases where the observations are not independent may exist provided that a Bernstein-von Mises theorem holds for the model at hand, which is the case for dependent observations if the likelihood model satisfies local asymptotic normality conditions Le  Cam (1953 Cam ( , 1986)). Note that since Assumption A. 4 is mostly used in Propositions 4 and 5 to guarantee that the log-likelihood ratio between the likelihood and the scaled likelihood of a subsample is bounded, the constant β in Proposition 7 is not a major concern. In addition, it is straightforward to see that this constant vanishes when the subset size grows faster than the full dataset, i.