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.     

An urn model to construct an efficient test procedure for response adaptive designs

We study the statistical performance of different tests for comparing the mean effect of two treatments. Given a reference classical test \({\mathcal {T}}_0\), we determine which sample size and proportion allocation guarantee to a test \({\mathcal {T}}\), based on response-adaptive design, to be better than \({\mathcal {T}}_0\), in terms of (a) higher power and (b) fewer subjects assigned to the inferior treatment. The adoption of a response-adaptive design to implement the random allocation procedure is necessary to ensure that both (a) and (b) are satisfied. In particular, we propose to use a Modified Randomly Reinforced Urn design and we show how to perform the model parameters selection for the purpose of this paper. Then, the opportunity of relaxing some assumptions on treatment response distributions is presented. Results of simulation studies on the test performance are reported and a real case study is analyzed.     

On comparison of dispersion matrices of estimators under a constrained linear model

We introduce some new mathematical tools in the analysis of dispersion matrices of the two well-known OLSEs and BLUEs under general linear models with parameter restrictions. We first establish some formulas for calculating the ranks and inertias of the differences of OLSEs’ and BLUEs’ dispersion matrices of parametric functions under the general linear model \({\mathscr {M}}= \{\mathbf{y}, \ \mathbf{X }\pmb {\beta }, \ \pmb {\Sigma }\}\) and the constrained model \({\mathscr {M}}_r = \{\mathbf{y}, \, \mathbf{X }\pmb {\beta }\, | \, \mathbf{A }\pmb {\beta }= \mathbf{b}, \ \pmb {\Sigma }\}\), where \(\mathbf{A }\pmb {\beta }= \mathbf{b}\) is a consistent linear matrix equation for the unknown parameter vector \(\pmb {\beta }\) to satisfy. As applications, we derive necessary and sufficient conditions for many equalities and inequalities of OLSEs’ and BLUEs’ dispersion matrices to hold under \({\mathscr {M}}\) and \({\mathscr {M}}_r\).     

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.     

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.     

New recursive estimators of the time-average variance constant

Estimation of the time-average variance constant (TAVC) of a stationary process plays a fundamental role in statistical inference for the mean of a stochastic process. Wu (2009) proposed an efficient algorithm to recursively compute the TAVC with \(O(1)\) memory and computational complexity. In this paper, we propose two new recursive TAVC estimators that can compute TAVC estimate with \(O(1)\) computational complexity. One of them is uniformly better than Wu’s estimator in terms of asymptotic mean squared error (MSE) at a cost of slightly higher memory complexity. The other preserves the \(O(1)\) memory complexity and is better then Wu’s estimator in most situations. Moreover, the first estimator is nearly optimal in the sense that its asymptotic MSE is \(2^{10/3}3^{-2} \fallingdotseq 1.12\) times that of the optimal off-line TAVC estimator.     

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).     

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.     

Transformation approaches of linear random-effects models

Assume that a linear random-effects model \(\mathbf{y}= \mathbf{X}\varvec{\beta }+ \varvec{\varepsilon }= \mathbf{X}(\mathbf{A}\varvec{\alpha }+ \varvec{\gamma }) + \varvec{\varepsilon }\) is transformed as \(\mathbf{T}\mathbf{y}= \mathbf{T}\mathbf{X}\varvec{\beta }+ \mathbf{T}\varvec{\varepsilon }= \mathbf{T}\mathbf{X}(\mathbf{A}\varvec{\alpha }+ \varvec{\gamma }) + \mathbf{T}\varvec{\varepsilon }\) by pre-multiplying a given matrix \(\mathbf{T}\) of arbitrary rank. These two models are not necessarily equivalent unless \(\mathbf{T}\) is of full column rank, and we have to work with this derived model in many situations. Because predictors/estimators of the parameter spaces under the two models are not necessarily the same, it is primary work to compare predictors/estimators in the two models and to establish possible links between the inference results obtained from two models. This paper presents a general algebraic approach to the problem of comparing best linear unbiased predictors (BLUPs) of parameter spaces in an original linear random-effects model and its transformations, and provides a group of fundamental and comprehensive results on mathematical and statistical properties of the BLUPs. In particular, we construct many equalities for the BLUPs under an original linear random-effects model and its transformations, and obtain necessary and sufficient conditions for the equalities to hold.     

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.     

Efficient $$\hbox {SMC}^2$$ schemes for stochastic kinetic models

Fitting stochastic kinetic models represented by Markov jump processes within the Bayesian paradigm is complicated by the intractability of the observed-data likelihood. There has therefore been considerable attention given to the design of pseudo-marginal Markov chain Monte Carlo algorithms for such models. However, these methods are typically computationally intensive, often require careful tuning and must be restarted from scratch upon receipt of new observations. Sequential Monte Carlo (SMC) methods on the other hand aim to efficiently reuse posterior samples at each time point. Despite their appeal, applying SMC schemes in scenarios with both dynamic states and static parameters is made difficult by the problem of particle degeneracy. A principled approach for overcoming this problem is to move each parameter particle through a Metropolis-Hastings kernel that leaves the target invariant. This rejuvenation step is key to a recently proposed \(\hbox {SMC}^2\) algorithm, which can be seen as the pseudo-marginal analogue of an idealised scheme known as iterated batch importance sampling. Computing the parameter weights in \(\hbox {SMC}^2\) requires running a particle filter over dynamic states to unbiasedly estimate the intractable observed-data likelihood up to the current time point. In this paper, we propose to use an auxiliary particle filter inside the \(\hbox {SMC}^2\) scheme. Our method uses two recently proposed constructs for sampling conditioned jump processes, and we find that the resulting inference schemes typically require fewer state particles than when using a simple bootstrap filter. Using two applications, we compare the performance of the proposed approach with various competing methods, including two global MCMC schemes.     

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 characterization of geometric distribution based on weak records

Let \({\{X_n, n\geq 1\}}\) be a sequence of independent and identically distributed non-degenerated random variables with common cumulative distribution function F. Suppose X ₁ is concentrated on 0, 1, . . . , N ≤ ∞ and P(X ₁ = 1) > 0. Let \({X_{U_w(n)}}\) be the n-th upper weak record value. In this paper we show that for any fixed m ≥ 2, X ₁ has Geometric distribution if and only if \({X_{U_{w}(m)}\mathop=\limits^d X_1+\cdots+X_m ,}\) where \({\underline{\underline{d}}}\) denotes equality in distribution. Our result is a generalization of the case m = 2 obtained by Ahsanullah (J Stat Theory Appl 8(1):5–16, 2009).     

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.     

A blocked Gibbs sampler for NGG-mixture models via a priori truncation

We define a new class of random probability measures, approximating the well-known normalized generalized gamma (NGG) process. Our new process is defined from the representation of NGG processes as discrete measures where the weights are obtained by normalization of the jumps of Poisson processes and the support consists of independent identically distributed location points, however considering only jumps larger than a threshold \(\varepsilon \). Therefore, the number of jumps of the new process, called \(\varepsilon \)-NGG process, is a.s. finite. A prior distribution for \(\varepsilon \) can be elicited. We assume such a process as the mixing measure in a mixture model for density and cluster estimation, and build an efficient Gibbs sampler scheme to simulate from the posterior. Finally, we discuss applications and performance of the model to two popular datasets, as well as comparison with competitor algorithms, the slice sampler and a posteriori truncation.     

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.     

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}\)).     

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.     

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.