Domains of convergence for the EM algorithm: a cautionary tale in a location estimation problem

The EM algorithm is a popular method for maximizing a likelihood in the presence of incomplete data. When the likelihood has multiple local maxima, the parameter space can be partitioned into domains of convergence, one for each local maximum. In this paper we investigate these domains for the location family generated by the t-distribution. We show that, perhaps somewhat surprisingly, these domains need not be connected sets. As an extreme case we give an example of a domain which consists of an infinite union of disjoint open intervals. Thus the convergence behaviour of the EM algorithm can be quite sensitive to the starting point.     

Maximum Likelihood Inference for Multivariate Frailty Models Using an Automated Monte Carlo EM Algorithm

We present a maximum likelihood estimation procedure for the multivariate frailty model. The estimation is based on a Monte Carlo EM algorithm. The expectation step is approximated by averaging over random samples drawn from the posterior distribution of the frailties using rejection sampling. The maximization step reduces to a standard partial likelihood maximization. We also propose a simple rule based on the relative change in the parameter estimates to decide on sample size in each iteration and a stopping time for the algorithm. An important new concept is acquiring absolute convergence of the algorithm through sample size determination and an efficient sampling technique. The method is illustrated using a rat carcinogenesis dataset and data on vase lifetimes of cut roses. The estimation results are compared with approximate inference based on penalized partial likelihood using these two examples. Unlike the penalized partial likelihood estimation, the proposed full maximum likelihood estimation method accounts for all the uncertainty while estimating standard errors for the parameters.     

Two‐Sample Bootstrap Hypothesis Tests for Three‐Dimensional Labelled Landmark Data

Abstract. We investigate resampling methodologies for testing the null hypothesis that two samples of labelled landmark data in three dimensions come from populations with a common mean reflection shape or mean reflection size‐and‐shape. The investigation includes comparisons between (i) two different test statistics that are functions of the projection onto tangent space of the data, namely the James statistic and an empirical likelihood statistic; (ii) bootstrap and permutation procedures; and (iii) three methods for resampling under the null hypothesis, namely translating in tangent space, resampling using weights determined by empirical likelihood and using a novel method to transform the original sample entirely within refection shape space. We present results of extensive numerical simulations, on which basis we recommend a bootstrap test procedure that we expect will work well in practise. We demonstrate the procedure using a data set of human faces, to test whether humans in different age groups have a common mean face shape.     

Robustness of a multiple ranking procedure: a monte carlo experiment illustrating design and analysis techniques

This case study demonstrates statistical design and analysis techniques applicable to any Monte Carlo or simulation experiment, namely a 2^7?3 experimental design, antithetic variates, sample size determination, analysis of variance, regression analysis, and simultaneous inference. The example is a Monte Carlo investigation of the robustness of Bechhofer and Blumenthal’s multiple ranking procedure (MRP). The investigation shows that their procedure works often, but not always. Factors that make it break down, are identified.     

Approximate Bayesian computation with composite score functions

Both approximate Bayesian computation (ABC) and composite likelihood methods are useful for Bayesian and frequentist inference, respectively, when the likelihood function is intractable. We propose to use composite likelihood score functions as summary statistics in ABC in order to obtain accurate approximations to the posterior distribution. This is motivated by the use of the score function of the full likelihood, and extended to general unbiased estimating functions in complex models. Moreover, we show that if the composite score is suitably standardised, the resulting ABC procedure is invariant to reparameterisations and automatically adjusts the curvature of the composite likelihood, and of the corresponding posterior distribution. The method is illustrated through examples with simulated data, and an application to modelling of spatial extreme rainfall data is discussed.     

Analytical uses of Kalman filtering in econometrics — A survey

This paper surveys the different uses of Kalman filtering in the estimation of statistical (econometric) models. The Kalman filter will be portrayed as (i) a natural generalization of exponential smoothing with a time-dependent smoothing factor, (ii) a recursive estimation technique for a variety of econometric models amenable to a state space formulation in particular for econometric models with time varying coefficients (iii) an instrument for the recursive calculation of the likelihood of the (constant) state space coefficients (iv) a means of helping to implement the scoring⁻ and EM-method for iteratively maximizing this likelihood (v) an analytical tool in asymptotic estimation theory. The concluding section points to the importance of Kalman filtering for alternatives to maximum⁻ likelihood estimation of state space parameters.     

Weak Convergence of the Regularization Path in Penalized M‐Estimation

Abstract. We consider a function defined as the pointwise minimization of a doubly index random process. We are interested in the weak convergence of the minimizer in the space of bounded functions. Such convergence results can be applied in the context of penalized M‐estimation, that is, when the random process to minimize is expressed as a goodness‐of‐fit term plus a penalty term multiplied by a penalty weight. This weight is called the regularization parameter and the minimizing function the regularization path. The regularization path can be seen as a collection of estimators indexed by the regularization parameter. We obtain a consistency result and a central limit theorem for the regularization path in a functional sense. Various examples are provided, including the ?¹‐regularization path for general linear models, the ?¹‐ or ?²‐regularization path of the least absolute deviation regression and the Akaike information criterion.     

Profile likelihood approaches for semiparametric copula and frailty models for clustered survival data

ABSTRACT

In clustered survival data, the dependence among individual survival times within a cluster has usually been described using copula models and frailty models. In this paper we propose a profile likelihood approach for semiparametric copula models with different cluster sizes. We also propose a likelihood ratio method based on profile likelihood for testing the absence of association parameter (i.e. test of independence) under the copula models, leading to the boundary problem of the parameter space. For this purpose, we show via simulation study that the proposed likelihood ratio method using an asymptotic chi-square mixture distribution performs well as sample size increases. We compare the behaviors of the two models using the profile likelihood approach under a semiparametric setting. The proposed method is demonstrated using two well-known data sets.     

Tail dependence of skew t-copulas

We examine tail behavior of skew t-copula in the bivariate case. The tail dependence coefficient is calculated for different skewing parameter values and compared with the corresponding coefficient for the t-copula. It is shown that depending on skewing parameter values, the tail dependence coefficient can differ considerably from the tail dependence of the t-copula. The speed of convergence of the estimator of tail dependence coefficient to its theoretical value is examined in a simulation experiment. Method of moments and maximum likelihood method are compared by simulation either. In the considered cases, maximum likelihood method converged faster to the theoretical value.     

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 new ranked set sampling protocol for the signed rank test

In this paper, we construct a new ranked set sampling protocol that maximizes the Pitman asymptotic efficiency of the signed rank test. The new sampling design is a function of the set size and independent order statistics. If the set size is odd and the underlying distribution is symmetric and unimodal, then the new sampling protocol quantifies only the middle observation. On the other hand, if the set size is even, the new sampling design quantifies the two middle observations. This data collection procedure for use in the signed rank test outperforms the data collection procedure in the standard ranked set sample. We show that the exact null distribution of the signed rank statistic W_RSS⁺ based on a data set generated by the new ranked set sample design for odd set sizes is the same as the null distribution of the simple random sample signed rank statistic W_SRS⁺ based on the same number of measured observations. For even set sizes, the exact null distribution of W_RSS⁺ is simulated.     

Inferring large graphs using $$\ell _1$$-penalized likelihood

We address the issue of recovering the structure of large sparse directed acyclic graphs from noisy observations of the system. We propose a novel procedure based on a specific formulation of the \(\ell _1\)-norm regularized maximum likelihood, which decomposes the graph estimation into two optimization sub-problems: topological structure and node order learning. We provide convergence inequalities for the graph estimator, as well as an algorithm to solve the induced optimization problem, in the form of a convex program embedded in a genetic algorithm. We apply our method to various data sets (including data from the DREAM4 challenge) and show that it compares favorably to state-of-the-art methods. This algorithm is available on CRAN as the R package GADAG.     

Selecting the normal population with the smallest variance: A restricted subset selection rule

Consider k( ? 2) normal populations whose means are all known or unknown and whose variances are unknown. Let σ²_[1] ? ??? ? σ_[k]² denote the ordered variances. Our goal is to select a non empty subset of the k populations whose size is at most m(1 ? m ? k ? 1) so that the population associated with the smallest variance (called the best population) is included in the selected subset with a guaranteed minimum probability P* whenever σ²_[2]/σ_[1]² ? δ* > 1, where P* and δ* are specified in advance of the experiment. Based on samples of size n from each of the populations, we propose and investigate a procedure called R_BCP. We also derive some asymptotic results for our procedure. Some comparisons with an earlier available procedure are presented in terms of the average subset sizes for selected slippage configurations based on simulations. The results are illustrated by an example.     

Is a Brownian Motion Skew?

We study the asymptotic behaviour of the maximum likelihood estimator corresponding to the observation of a trajectory of a skew Brownian motion, through a uniform time discretization. We characterize the speed of convergence and the limiting distribution when the step size goes to zero, which in this case are non‐classical, under the null hypothesis of the skew Brownian motion being an usual Brownian motion. This allows to design a test on the skewness parameter. We show that numerical simulations can be easily performed to estimate the skewness parameter and provide an application in Biology.     

Decomposition for multivariate extremal processes

The probability distribution of an extremal process in R^d with independent max-increments is completely determined by its distribution function. The df of an extremal process is similar to the cdf of a random vector. It is a monotone function on (0, ∞) × R^d with values in the interval [0,1]. On the other hand the probability distribution of an extremal process is a probability measure on the space of sample functions. That is the space of all increasing right continuous functions y: (0, ∞) → R^d with the topology of weak convergence. A sequence of extremal processes converges in law if the probability distributions converge weakly. This is shown to be equivalent to weak convergence of the df's.

An extremal process Y: [0, ∞) → R^d is generated by a point process on the space [0, ∞) × [-∞, ∞)^d and has a decomposition Y = X v Z as the maximum of two independent extremal processes with the same lower curve as the original process. The process X is the continuous part and Z contains the fixed discontinuities of the process Y. For a real valued extremal process the decomposition is unique: for a multivariate extremal process uniqueness breaks down due to blotting.     

Estimating nonlinear additive models with nonstationarities and correlated errors

In this paper, we study a nonparametric additive regression model suitable for a wide range of time series applications. Our model includes a periodic component, a deterministic time trend, various component functions of stochastic explanatory variables, and an AR(p) error process that accounts for serial correlation in the regression error. We propose an estimation procedure for the nonparametric component functions and the parameters of the error process based on smooth backfitting and quasimaximum likelihood methods. Our theory establishes convergence rates and the asymptotic normality of our estimators. Moreover, we are able to derive an oracle‐type result for the estimators of the AR parameters: Under fairly mild conditions, the limiting distribution of our parameter estimators is the same as when the nonparametric component functions are known. Finally, we illustrate our estimation procedure by applying it to a sample of climate and ozone data collected on the Antarctic Peninsula.     

Convergence Rate of Strong Consistency of the Maximum Likelihood Estimator in Exponential Family Nonlinear Models

This article proposes some regularity conditions. On the basis of the proposed regularity conditions, we show the strong consistency of the maximum likelihood estimator (MLE) in exponential family nonlinear models (EFNM) and give its convergence rate. In an important case, we obtain the convergence rate O(n ^?1/2(log log n)^1/2)—the rate as that in the Law of the Iterated Logarithm (LIL) for iid partial sums and thus cannot be improved anymore.     

On empirical Bayes simultaneous selection procedures for comparing normal populations with a standard

In this paper, we derive statistical selection procedures to partition k normal populations into ‘good’ or ‘bad’ ones, respectively, using the nonparametric empirical Bayes approach. The relative regret risk of a selection procedure is used as a measure of its performance. We establish the asymptotic optimality of the proposed empirical Bayes selection procedures and investigate the associated rates of convergence. Under a very mild condition, the proposed empirical Bayes selection procedures are shown to have rates of convergence of order close to O(k^−1/2) where k is the number of populations involved in the selection problem. With further strong assumptions, the empirical Bayes selection procedures have rates of convergence of order O(k^{−α(r−1)/(2r+1)}), where 1<α<2 and r is an integer greater than 2.     

Stochastic proximal-gradient algorithms for penalized mixed models

Motivated by penalized likelihood maximization in complex models, we study optimization problems where neither the function to optimize nor its gradient has an explicit expression, but its gradient can be approximated by a Monte Carlo technique. We propose a new algorithm based on a stochastic approximation of the proximal-gradient (PG) algorithm. This new algorithm, named stochastic approximation PG (SAPG) is the combination of a stochastic gradient descent step which—roughly speaking—computes a smoothed approximation of the gradient along the iterations, and a proximal step. The choice of the step size and of the Monte Carlo batch size for the stochastic gradient descent step in SAPG is discussed. Our convergence results cover the cases of biased and unbiased Monte Carlo approximations. While the convergence analysis of some classical Monte Carlo approximation of the gradient is already addressed in the literature (see Atchadé et al. in J Mach Learn Res 18(10):1–33, 2017), the convergence analysis of SAPG is new. Practical implementation is discussed, and guidelines to tune the algorithm are given. The two algorithms are compared on a linear mixed effect model as a toy example. A more challenging application is proposed on nonlinear mixed effect models in high dimension with a pharmacokinetic data set including genomic covariates. To our best knowledge, our work provides the first convergence result of a numerical method designed to solve penalized maximum likelihood in a nonlinear mixed effect model.

     

Empirical likelihood-based serial correlation testing in partially varying coefficient single-index models

ABSTRACT

Partially varying coefficient single-index models (PVCSIM) are a class of semiparametric regression models. One important assumption is that the model error is independently and identically distributed, which may contradict with the reality in many applications. For example, in the economical and financial applications, the observations may be serially correlated over time. Based on the empirical likelihood technique, we propose a procedure for testing the serial correlation of random error in PVCSIM. Under some regular conditions, we show that the proposed empirical likelihood ratio statistic asymptotically follows a standard χ² distribution. We also present some numerical studies to illustrate the performance of our proposed testing procedure.