A discrepancy bound for deterministic acceptance-rejection samplers beyond $$N^{-1/2}$$ in dimension 1

In this paper we consider an acceptance-rejection (AR) sampler based on deterministic driver sequences. We prove that the discrepancy of an N element sample set generated in this way is bounded by \(\mathcal {O} (N^{-2/3}\log N)\), provided that the target density is twice continuously differentiable with non-vanishing curvature and the AR sampler uses the driver sequence \(\mathcal {K}_M= \{( j \alpha , j \beta ) ~~ mod~~1 \mid j = 1,\ldots ,M\},\) where \(\alpha ,\beta \) are real algebraic numbers such that \(1,\alpha ,\beta \) is a basis of a number field over \(\mathbb {Q}\) of degree 3. For the driver sequence \(\mathcal {F}_k= \{ ({j}/{F_k}, \{{jF_{k-1}}/{F_k}\} ) \mid j=1,\ldots , F_k\},\) where \(F_k\) is the k-th Fibonacci number and \(\{x\}=x-\lfloor x \rfloor \) is the fractional part of a non-negative real number x, we can remove the \(\log \) factor to improve the convergence rate to \(\mathcal {O}(N^{-2/3})\), where again N is the number of samples we accepted. We also introduce a criterion for measuring the goodness of driver sequences. The proposed approach is numerically tested by calculating the star-discrepancy of samples generated for some target densities using \(\mathcal {K}_M\) and \(\mathcal {F}_k\) as driver sequences. These results confirm that achieving a convergence rate beyond \(N^{-1/2}\) is possible in practice using \(\mathcal {K}_M\) and \(\mathcal {F}_k\) as driver sequences in the acceptance-rejection sampler.     

Wavelet regression estimations with strong mixing data

Using a wavelet basis, we establish in this paper upper bounds of wavelet estimation on \( L^{p}({\mathbb {R}}^{d}) \) risk of regression functions with strong mixing data for \( 1\le p<\infty \). In contrast to the independent case, these upper bounds have different analytic formulae for \(p\in [1, 2]\) and \(p\in (2, +\infty )\). For \(p=2\), it turns out that our result reduces to a theorem of Chaubey et al. (J Nonparametr Stat 25:53–71, 2013); and for \(d=1\) and \(p=2\), it becomes the corresponding theorem of Chaubey and Shirazi (Commun Stat Theory Methods 44:885–899, 2015).     

Rate of uniform consistency for a class of mode regression on functional stationary ergodic data

The aim of this paper is to study the asymptotic properties of a class of kernel conditional mode estimates whenever functional stationary ergodic data are considered. To be more precise on the matter, in the ergodic data setting, we consider a random elements (X, Z) taking values in some semi-metric abstract space \(E\times F\). For a real function \(\varphi \) defined on the space F and \(x\in E\), we consider the conditional mode of the real random variable \(\varphi (Z)\) given the event “\(X=x\)”. While estimating the conditional mode function, say \(\theta _\varphi (x)\), using the well-known kernel estimator, we establish the strong consistency with rate of this estimate uniformly over Vapnik–Chervonenkis classes of functions \(\varphi \). Notice that the ergodic setting offers a more general framework than the usual mixing structure. Two applications to energy data are provided to illustrate some examples of the proposed approach in time series forecasting framework. The first one consists in forecasting the daily peak of electricity demand in France (measured in Giga-Watt). Whereas the second one deals with the short-term forecasting of the electrical energy (measured in Giga-Watt per Hour) that may be consumed over some time intervals that cover the peak demand.     

Multilevel particle filters: normalizing constant estimation

In this article, we introduce two new estimates of the normalizing constant (or marginal likelihood) for partially observed diffusion (POD) processes, with discrete observations. One estimate is biased but non-negative and the other is unbiased but not almost surely non-negative. Our method uses the multilevel particle filter of Jasra et al. (Multilevel particle lter, arXiv:1510.04977, 2015). We show that, under assumptions, for Euler discretized PODs and a given \(\varepsilon >0\) in order to obtain a mean square error (MSE) of \({\mathcal {O}}(\varepsilon ^2)\) one requires a work of \({\mathcal {O}}(\varepsilon ^{-2.5})\) for our new estimates versus a standard particle filter that requires a work of \({\mathcal {O}}(\varepsilon ^{-3})\). Our theoretical results are supported by numerical simulations.     

Stochastically optimal bootstrap sample size for shrinkage-type statistics

In nonregular problems where the conventional \(n\) out of \(n\) bootstrap is inconsistent, the \(m\) out of \(n\) bootstrap provides a useful remedy to restore consistency. Conventionally, optimal choice of the bootstrap sample size \(m\) is taken to be the minimiser of a frequentist error measure, estimation of which has posed a major difficulty hindering practical application of the \(m\) out of \(n\) bootstrap method. Relatively little attention has been paid to a stronger, stochastic, version of the optimal bootstrap sample size, defined as the minimiser of an error measure calculated directly from the observed sample. Motivated by this stronger notion of optimality, we develop procedures for calculating the stochastically optimal value of \(m\). Our procedures are shown to work under special forms of Edgeworth-type expansions which are typically satisfied by statistics of the shrinkage type. Theoretical and empirical properties of our methods are illustrated with three examples, namely the James–Stein estimator, the ridge regression estimator and the post-model-selection regression estimator.     

Conditional density estimation using the local Gaussian correlation

Let \(\mathbf {X} = (X_1,\ldots ,X_p)\) be a stochastic vector having joint density function \(f_{\mathbf {X}}(\mathbf {x})\) with partitions \(\mathbf {X}_1 = (X_1,\ldots ,X_k)\) and \(\mathbf {X}_2 = (X_{k+1},\ldots ,X_p)\). A new method for estimating the conditional density function of \(\mathbf {X}_1\) given \(\mathbf {X}_2\) is presented. It is based on locally Gaussian approximations, but simplified in order to tackle the curse of dimensionality in multivariate applications, where both response and explanatory variables can be vectors. We compare our method to some available competitors, and the error of approximation is shown to be small in a series of examples using real and simulated data, and the estimator is shown to be particularly robust against noise caused by independent variables. We also present examples of practical applications of our conditional density estimator in the analysis of time series. Typical values for k in our examples are 1 and 2, and we include simulation experiments with values of p up to 6. Large sample theory is established under a strong mixing condition.     

Bayes and maximum likelihood for $$L^1$$-Wasserstein deconvolution of Laplace mixtures

We consider the problem of recovering a distribution function on the real line from observations additively contaminated with errors following the standard Laplace distribution. Assuming that the latent distribution is completely unknown leads to a nonparametric deconvolution problem. We begin by studying the rates of convergence relative to the \(L^2\)-norm and the Hellinger metric for the direct problem of estimating the sampling density, which is a mixture of Laplace densities with a possibly unbounded set of locations: the rate of convergence for the Bayes’ density estimator corresponding to a Dirichlet process prior over the space of all mixing distributions on the real line matches, up to a logarithmic factor, with the \(n^{-3/8}\log ^{1/8}n\) rate for the maximum likelihood estimator. Then, appealing to an inversion inequality translating the \(L^2\)-norm and the Hellinger distance between general kernel mixtures, with a kernel density having polynomially decaying Fourier transform, into any \(L^p\)-Wasserstein distance, \(p\ge 1\), between the corresponding mixing distributions, provided their Laplace transforms are finite in some neighborhood of zero, we derive the rates of convergence in the \(L^1\)-Wasserstein metric for the Bayes’ and maximum likelihood estimators of the mixing distribution. Merging in the \(L^1\)-Wasserstein distance between Bayes and maximum likelihood follows as a by-product, along with an assessment on the stochastic order of the discrepancy between the two estimation procedures.     

Hypotheses testing about the drift parameter in linear stochastic differential equation driven by stable processes

In this paper, we consider the problem of hypotheses testing about the drift parameter \(\theta \) in the process \(\text {d}Y^{\delta }_{t} = \theta \dot{f}(t)Y^{\delta }_{t}\text {d}t + b(t)\text {d}L^{\delta }_{t}\) driven by symmetric \(\delta \)-stable Lévy process \(L^{\delta }_{t}\) with \(\dot{f}(t)\) being the derivative of a known increasing function f(t) and b(t) being known as well. We consider the hypotheses testing \(H_{0}: \theta \le 0\) and \(K_{0}: \theta =0\) against the alternatives \(H_{1}: \theta >0\) and \(K_{1}: \theta \ne 0\), respectively. For these hypotheses, we propose inverse methods, which are motivated by sequential approach, based on the first hitting time of the observed process (or its absolute value) to a pre-specified boundary or two boundaries until some given time. The applicability of these methods is illustrated. For the case \(Y^{\delta }_{0}=0\), we are able to calculate the values of boundaries and finite observed times more directly. We are able to show the consistencies of proposed tests for \(Y^{\delta }_{0}\ge 0\) with \(\delta \in (1,2]\) and for \(Y^{\delta }_{0}=0\) with \(\delta \in (0,2]\) under quite mild conditions.     

Some lower bounds of centered $$L_2$$-discrepancy of $$2^{s-k}$$2s-k designs and their complementary designs

The indicator function is an effective tool in studying factorial designs. This paper presents some lower bounds of centered \(L_2\)-discrepancy through indicator function. Some new lower bounds of centered \(L_2\)-discrepancy for \(2^{s-k}\) designs and their complementary designs are given. Numerical results show that our lower bounds are tight and better than the existing results.     

A search of maximum generalized resolution quaternary-code designs via integer linear programming

Quaternary-code (QC) designs, an attractive class of nonregular fractional factorial designs, have received much attention due to their theoretical elegance and practical applicability. Some recent works of QC designs revealed their good properties over their regular counterparts under commonly used criteria. We develop an optimization tool that can maximize the generalized resolution of a QC design of a given size. The problem can be recast as an integer linear programming (ILP) problem through a linear simplification that combines the \(k\)- and \(a\)-equations, even though the generalized resolution does not linearly depend on the aliasing indexes. The ILP surprisingly improves a class of \((1/16)\)th-fraction QC designs with higher generalized resolutions. It also applies to obtain some \((1/64)\)th-fraction QC designs with maximum generalized resolutions, and these QC designs generally have higher generalized resolutions than the regular designs of the same size.     

Computing the log concave NPMLE for interval censored data

In analyzing interval censored data, a non-parametric estimator is often desired due to difficulties in assessing model fits. Because of this, the non-parametric maximum likelihood estimator (NPMLE) is often the default estimator. However, the estimates for values of interest of the survival function, such as the quantiles, have very large standard errors due to the jagged form of the estimator. By forcing the estimator to be constrained to the class of log concave functions, the estimator is ensured to have a smooth survival estimate which has much better operating characteristics than the unconstrained NPMLE, without needing to specify a parametric family or smoothing parameter. In this paper, we first prove that the likelihood can be maximized under a finite set of parameters under mild conditions, although the log likelihood function is not strictly concave. We then present an efficient algorithm for computing a local maximum of the likelihood function. Using our fast new algorithm, we present evidence from simulated current status data suggesting that the rate of convergence of the log-concave estimator is faster (between \(n^{2/5}\) and \(n^{1/2}\)) than the unconstrained NPMLE (between \(n^{1/3}\) and \(n^{1/2}\)).     

Random projections for Bayesian regression

This article deals with random projections applied as a data reduction technique for Bayesian regression analysis. We show sufficient conditions under which the entire d-dimensional distribution is approximately preserved under random projections by reducing the number of data points from n to \(k\in O({\text {poly}}(d/\varepsilon ))\) in the case \(n\gg d\). Under mild assumptions, we prove that evaluating a Gaussian likelihood function based on the projected data instead of the original data yields a \((1+O(\varepsilon ))\)-approximation in terms of the \(\ell _2\) Wasserstein distance. Our main result shows that the posterior distribution of Bayesian linear regression is approximated up to a small error depending on only an \(\varepsilon \)-fraction of its defining parameters. This holds when using arbitrary Gaussian priors or the degenerate case of uniform distributions over \(\mathbb {R}^d\) for \(\beta \). Our empirical evaluations involve different simulated settings of Bayesian linear regression. Our experiments underline that the proposed method is able to recover the regression model up to small error while considerably reducing the total running time.     

Automated selection of <Emphasis Type="Italic">r</Emphasis> for the <Emphasis Type="Italic">r</Emphasis> largest order statistics approach with adjustment for sequential testing

The r largest order statistics approach is widely used in extreme value analysis because it may use more information from the data than just the block maxima. In practice, the choice of r is critical. If r is too large, bias can occur; if too small, the variance of the estimator can be high. The limiting distribution of the r largest order statistics, denoted by GEV\(_r\), extends that of the block maxima. Two specification tests are proposed to select r sequentially. The first is a score test for the GEV\(_r\) distribution. Due to the special characteristics of the GEV\(_r\) distribution, the classical chi-square asymptotics cannot be used. The simplest approach is to use the parametric bootstrap, which is straightforward to implement but computationally expensive. An alternative fast weighted bootstrap or multiplier procedure is developed for computational efficiency. The second test uses the difference in estimated entropy between the GEV\(_r\) and GEV\(_{r-1}\) models, applied to the r largest order statistics and the \(r-1\) largest order statistics, respectively. The asymptotic distribution of the difference statistic is derived. In a large scale simulation study, both tests held their size and had substantial power to detect various misspecification schemes. A new approach to address the issue of multiple, sequential hypotheses testing is adapted to this setting to control the false discovery rate or familywise error rate. The utility of the procedures is demonstrated with extreme sea level and precipitation data.     

Point process-based Monte Carlo estimation

This paper addresses the issue of estimating the expectation of a real-valued random variable of the form \(X = g(\mathbf {U})\) where g is a deterministic function and \(\mathbf {U}\) can be a random finite- or infinite-dimensional vector. Using recent results on rare event simulation, we propose a unified framework for dealing with both probability and mean estimation for such random variables, i.e. linking algorithms such as Tootsie Pop Algorithm or Last Particle Algorithm with nested sampling. Especially, it extends nested sampling as follows: first the random variable X does not need to be bounded any more: it gives the principle of an ideal estimator with an infinite number of terms that is unbiased and always better than a classical Monte Carlo estimator—in particular it has a finite variance as soon as there exists \(k \in \mathbb {R}> 1\) such that \({\text {E}}\left[ X^k \right] < \infty \). Moreover we address the issue of nested sampling termination and show that a random truncation of the sum can preserve unbiasedness while increasing the variance only by a factor up to 2 compared to the ideal case. We also build an unbiased estimator with fixed computational budget which supports a Central Limit Theorem and discuss parallel implementation of nested sampling, which can dramatically reduce its running time. Finally we extensively study the case where X is heavy-tailed.     

On periodic time-varying bilinear processes: structure and asymptotic inference

This paper is devoted to the bilinear time series models with periodic-varying coefficients \(\left( { PBL}\right) \). So, firstly conditions ensuring the existence of periodic stationary solutions of the \({ PBL}\) and the existence of higher-order moments of such solutions are given. A distribution free approach to the parameter estimation of \({ PBL}\) is presented. The proposed method relies on minimum distance estimator based on the first and second order empirical moments of the observed process. Consistency and asymptotic normality of the estimator are discussed. Examples and Monte Carlo simulation results illustrate the practical relevancy of our general theoretical results are presented.     

Investigation of the widely applicable Bayesian information criterion

The widely applicable Bayesian information criterion (WBIC) is a simple and fast approximation to the model evidence that has received little practical consideration. WBIC uses the fact that the log evidence can be written as an expectation, with respect to a powered posterior proportional to the likelihood raised to a power \(t^*\in {(0,1)}\), of the log deviance. Finding this temperature value \(t^*\) is generally an intractable problem. We find that for a particular tractable statistical model that the mean squared error of an optimally-tuned version of WBIC with correct temperature \(t^*\) is lower than an optimally-tuned version of thermodynamic integration (power posteriors). However in practice WBIC uses the a canonical choice of \(t=1/\log (n)\). Here we investigate the performance of WBIC in practice, for a range of statistical models, both regular models and singular models such as latent variable models or those with a hierarchical structure for which BIC cannot provide an adequate solution. Our findings are that, generally WBIC performs adequately when one uses informative priors, but it can systematically overestimate the evidence, particularly for small sample sizes.     

Finie mixures of canonical fundamenal skew $$$$-disribuions

This paper introduces a finite mixture of canonical fundamental skew \(t\) (CFUST) distributions for a model-based approach to clustering where the clusters are asymmetric and possibly long-tailed (in: Lee and McLachlan, arXiv:1401.8182 [statME], 2014b). The family of CFUST distributions includes the restricted multivariate skew \(t\) and unrestricted multivariate skew \(t\) distributions as special cases. In recent years, a few versions of the multivariate skew \(t\) (MST) mixture model have been put forward, together with various EM-type algorithms for parameter estimation. These formulations adopted either a restricted or unrestricted characterization for their MST densities. In this paper, we examine a natural generalization of these developments, employing the CFUST distribution as the parametric family for the component distributions, and point out that the restricted and unrestricted characterizations can be unified under this general formulation. We show that an exact implementation of the EM algorithm can be achieved for the CFUST distribution and mixtures of this distribution, and present some new analytical results for a conditional expectation involved in the E-step.     

Semi-parametric bivariate polychotomous ordinal regression

A pair of polychotomous random variables \((Y_1,Y_2)^\top =:{\varvec{Y}}\), where each \(Y_j\) has a totally ordered support, is studied within a penalized generalized linear model framework. We deal with a triangular generating process for \({\varvec{Y}}\), a structure that has been employed in the literature to control for the presence of residual confounding. Differently from previous works, however, the proposed model allows for a semi-parametric estimation of the covariate-response relationships. In this way, the risk of model mis-specification stemming from the imposition of fixed-order polynomial functional forms is also reduced. The proposed estimation methods and related inferential results are finally applied to study the effect of education on alcohol consumption among young adults in the UK.     

A comparison of efficient approximations for a weighted sum of chi-squared random variables

In many applications, the cumulative distribution function (cdf) \(F_{Q_N}\) of a positively weighted sum of N i.i.d. chi-squared random variables \(Q_N\) is required. Although there is no known closed-form solution for \(F_{Q_N}\), there are many good approximations. When computational efficiency is not an issue, Imhof’s method provides a good solution. However, when both the accuracy of the approximation and the speed of its computation are a concern, there is no clear preferred choice. Previous comparisons between approximate methods could be considered insufficient. Furthermore, in streaming data applications where the computation needs to be both sequential and efficient, only a few of the available methods may be suitable. Streaming data problems are becoming ubiquitous and provide the motivation for this paper. We develop a framework to enable a much more extensive comparison between approximate methods for computing the cdf of weighted sums of an arbitrary random variable. Utilising this framework, a new and comprehensive analysis of four efficient approximate methods for computing \(F_{Q_N}\) is performed. This analysis procedure is much more thorough and statistically valid than previous approaches described in the literature. A surprising result of this analysis is that the accuracy of these approximate methods increases with N.     

Optimal regular graph designs

A typical problem in optimal design theory is finding an experimental design that is optimal with respect to some criteria in a class of designs. The most popular criteria include the A- and D-criteria. Regular graph designs occur in many optimality results, and if the number of blocks is large enough, an A-optimal (or D-optimal) design is among them (if any exist). To explore the landscape of designs with a large number of blocks, we introduce extensions of regular graph designs. These are constructed by adding the blocks of a balanced incomplete block design repeatedly to the original design. We present the results of an exact computer search for the best regular graph designs and the best extended regular graph designs with up to 20 treatments v, block size \(k \le 10\) and replication r \(\le 10\) and \(r(k-1)-(v-1)\lfloor r(k-1)/(v-1)\rfloor \le 9\).