In this paper, we use the class of Wasserstein metrics to study asymptotic properties of posterior distributions. Our first goal is to provide sufficient conditions for posterior consistency. In addition to the well-known Schwartz’s Kullback–Leibler condition on the prior, the true distribution and most probability measures in the support of the prior are required to possess moments up to an order which is determined by the order of the Wasserstein metric. We further investigate convergence rates of the posterior distributions for which we need stronger moment conditions. The required tail conditions are sharp in the sense that the posterior distribution may be inconsistent or contract slowly to the true distribution without these conditions. Our study involves techniques that build on recent advances on Wasserstein convergence of empirical measures. We apply the results to some examples including a Dirichlet process mixture prior and conduct a simulation study for further illustration.

In network data analysis, multivariate spatial autoregression (MSAR) models may be used to analyze the autocorrelation among multiple responses. With large-scale networks, the estimation for MSAR on the entire network is computationally expensive. In this case, the subsampling method could be adopted. This approach selects a sample of nodes and then uses the estimate based on the sample to approximate the estimate on the full data. However, traditional sampling methods cannot obtain unbiased parameter estimates. Considering the second-order friend information of sampled nodes, we propose the crawling subsampling (CS) method for a general framework. Thereafter, based on the sampled data only, we construct the least-squares objective function. Under certain conditions, the computational complexity of optimizing the objective function is linear with the sample size $n_s$. The identification condition for the parameters on the sampled network is theoretically provided. The sample size order requirement is provided, and the asymptotic properties of the least-squares estimators are investigated. The numerical results for the simulated and real data are presented to demonstrate the performance of the proposed CS method and least-squares estimator for the MSAR model.

We consider survival data from a population with cured subjects in the presence of mismeasured covariates. We use the mixture cure model to account for the individuals that will never experience the event and at the same time distinguish between the effect of the covariates on the cure probabilities and on survival times. In particular, for practical applications, it seems of interest to assume a logistic form of the incidence and a Cox proportional hazards model for the latency. To correct the estimators for the bias introduced by the measurement error, we use the simex algorithm, which is a very general simulation based method. It essentially estimates this bias by introducing additional error to the data and then recovers bias corrected estimators through an extrapolation approach. The estimators are shown to be consistent and asymptotically normally distributed when the true extrapolation function is known. We investigate their finite sample performance through a simulation study and apply the proposed method to analyse the effect of the prostate specific antigen (PSA) on patients with prostate cancer.

Kernel ridge regression (KRR) is a well-known and popular nonparametric regression approach with many desirable properties, including minimax rate-optimality in estimating functions that belong to common reproducing kernel Hilbert spaces (RKHS). The approach, however, is computationally intensive for large data sets, due to the need to operate on a dense $n\times n$ kernel matrix, where n is the sample size. Recently, various approximation schemes for solving KRR have been considered, and some analyzed. Some approaches such as Nyström approximation and sketching have been shown to preserve the rate optimality of KRR. In this paper, we consider the simplest approximation, namely, spectrally truncating the kernel matrix to its largest $r<n$ eigenvalues. We derive an exact expression for the maximum risk of this truncated KRR, over the unit ball of the RKHS. This result can be used to study the exact trade-off between the level of spectral truncation and the regularization parameter. We show that, as long as the RKHS is infinite-dimensional, there is a threshold on r, above which, the spectrally-truncated KRR surprisingly outperforms the full KRR in terms of the minimax risk, where the minimum is taken over the regularization parameter. This strengthens the existing results on approximation schemes, by showing that not only one does not lose in terms of the rates, truncation can in fact improve the performance, for all finite samples (above the threshold). Moreover, we show that the implicit regularization achieved by spectral truncation is not a substitute for Hilbert norm regularization. Both are needed to achieve the best performance.

In the common time series model $X_{i,n}=\mathrm{\mu }(i∕n)+{\mathit{\epsilon }}_{i,n}$ with non-stationary errors we consider the problem of detecting a significant deviation of the mean function μ from a benchmark $g(\mathrm{\mu })$ (such as the initial value $\mathrm{\mu }(0)$ or the average trend ${\int }_0^1\mathrm{\mu }(t)dt$). The problem is motivated by a more realistic modelling of change point analysis, where one is interested in identifying relevant deviations in a smoothly varying sequence of means ${(\mathrm{\mu }(i∕n))}_{i=1,\dots ,n}$ and cannot assume that the sequence is piecewise constant. A test for this type of hypotheses is developed using an appropriate estimator for the integrated squared deviation of the mean function and the threshold. By a new concept of self-normalization adapted to non-stationary processes an asymptotically pivotal test for the hypothesis of a relevant deviation is constructed. The results are illustrated by means of a simulation study and a data example.

This paper is on the normal approximation of singular subspaces when the noise matrix has i.i.d. entries. Our contributions are three-fold. First, we derive an explicit representation formula of the empirical spectral projectors. The formula is neat and holds for deterministic matrix perturbations. Second, we calculate the expected projection distance between the empirical singular subspaces and true singular subspaces. Our method allows obtaining arbitrary k-th order approximation of the expected projection distance. Third, we prove the non-asymptotical normal approximation of the projection distance with different levels of bias corrections. By the $\lceil log(d_1+d_2)\rceil$-th order bias corrections, the asymptotical normality holds under optimal signal-to-noise ratio (SNR) condition where $d_1$ and $d_2$ denote the matrix sizes. In addition, it shows that higher order approximations are unnecessary when $|d_1-d_2|=O({(d_1+d_2)}^{1∕2})$. Finally, we provide comprehensive simulation results to merit our theoretic discoveries.

Unlike the existing results, our approach is non-asymptotical and the convergence rates are established. Our method allows the rank r to diverge as fast as $o({(d_1+d_2)}^{1∕3})$. Moreover, our method requires no eigen-gap condition (except the SNR) and no constraints between $d_1$ and $d_2$.

We introduce noncausal counting processes, defined by time-reversing an INAR(1) process, a non-INAR(1) Markov affine counting process, or a random coefficient INAR(1) [RCINAR(1)] process. The noncausal processes are shown to be generically time irreversible and their calendar time dynamic properties are unreplicable by existing causal models. In particular, they allow for locally bubble-like explosion, while at the same time preserving stationarity. Many of these processes have also closed form calendar time conditional predictive distribution, and allow for a simple queuing interpretation, similar as their causal counterparts.

In this paper, a general stochastic optimization procedure is studied, unifying several variants of the stochastic gradient descent such as, among others, the stochastic heavy ball method, the Stochastic Nesterov Accelerated Gradient algorithm (S-NAG), and the widely used Adam algorithm. The algorithm is seen as a noisy Euler discretization of a non-autonomous ordinary differential equation, recently introduced by Belotto da Silva and Gazeau, which is analyzed in depth. Assuming that the objective function is non-convex and differentiable, the stability and the almost sure convergence of the iterates to the set of critical points are established. A noteworthy special case is the convergence proof of S-NAG in a non-convex setting. Under some assumptions, the convergence rate is provided under the form of a Central Limit Theorem. Finally, the non-convergence of the algorithm to undesired critical points, such as local maxima or saddle points, is established. Here, the main ingredient is a new avoidance of traps result for non-autonomous settings, which is of independent interest.

In this correction, we provide an acknowledgment for the source of several images in the paper: J. Hendricks, C. Neumann, C.P. Saunders (2021) “Quantification of the weight of fingerprint evidence using a ROC-based Approximate Bayesian Computation algorithm for model selection”, Electronic Journal of Statistics (15)1, pp. 1228–1262.

Incorporating information about the target distribution in proposal mechanisms generally produces efficient Markov chain Monte Carlo algorithms (or at least, algorithms that are more efficient than uninformed counterparts). For instance, it has proved successful to incorporate gradient information in fixed-dimensional algorithms, as seen with algorithms such as Hamiltonian Monte Carlo. In trans-dimensional algorithms, Green (2003) recommended to sample the parameter proposals during model switches from normal distributions with informative means and covariance matrices. These proposal distributions can be viewed as asymptotic approximations to the parameter distributions, where the limit is with regard to the sample size. Models are typically proposed using uninformed uniform distributions. In this paper, we build on the approach of Zanella (2020) for discrete spaces to incorporate information about neighbouring models. We rely on approximations to posterior model probabilities that are asymptotically exact. We prove that, in some scenarios, the samplers combining this approach with that of Green (2003) behave like ideal ones that use the exact model probabilities and sample from the correct parameter distributions, in the large-sample regime. We show that the implementation of the proposed samplers is straightforward in some cases. The methodology is applied to a real-data example. The code is available online.

This article presents a novel method for prediction of stationary functional time series, in particular for trajectories that share a similar pattern but display variable phases. The limitation of most of the existing prediction methodologies for functional time series is that they only consider vertical variation (amplitude, scale, or vertical shift). To overcome this limitation, we develop a shape-preserving (SP) prediction method that incorporates both vertical and horizontal variation. One major advantage of our proposed method is the ability to preserve the shape of functions. Moreover, our proposed SP method does not involve unnatural transformations and can be easily implemented using existing software packages. The utility of the SP method is demonstrated in the analysis of non-metanic hydrocarbons (NMHC) concentration. The analysis demonstrates that the prediction by the SP method captures the common pattern better than the existing prediction methods and also provides competitive prediction accuracy.

This paper derives a new randomization procedure by conditioning Efron’s (1971) biased coin design to a prespecified final balance. The new procedure remains a function of the original bias parameter which now controls the probability of intermediate balance rather than final balance. As the sample size increases, the design’s selection bias and intermediate balance are similar to those of the original biased coin, but unlike the original biased coin it always guarantees final balance. It is also shown that the permuted block design for equal allocation is a special case of the new procedure when used in blocks. The latter can substitute the permuted blocks with the added benefit of reducing the expected number of deterministic assignments. The new design is also noteworthy since it shows that a randomization procedure with new properties can be obtained by conditioning an existing one to a subset in its allocation space. New relationships among existing designs can be established in the process, further elucidating the protean nature of randomization.

Many studies have been conducted on flows of probability measures, often in terms of gradient flows. We utilize a generalized notion of derivatives with respect to time to model the instantaneous evolution of empirically observed one-dimensional distributions that vary over time and develop consistent estimates for these derivatives. Employing local Fréchet regression and working in local tangent spaces with regard to the Wasserstein metric, we derive the rate of convergence of the proposed estimators. The resulting time dynamics are illustrated with time-varying distribution data that include yearly income distributions and the evolution of mortality over calendar years.

We investigate a connection between spatial statistics and statistical physics to obtain new covariance functions with direct physical interpretation for spatial random fields. These covariance functions are based on the exponential Boltzmann-Gibbs representation and use an energy functional to represent interactions between the values of the random field at different points in space. This formulation results in closed-form generalized covariance functions, which display infinite variance in Euclidean spaces of dimension larger than one. We propose regularization schemes in real and reciprocal (spectral) space that lead to well-behaved covariance structures. The real-space regularization parameter allows a continuous interpolation between the Boltzmann-Gibbs covariance and the exponential covariance. We also propose discretized approximations on regular grids, and we show that they represent reparametrized versions of the well-known Besag and Leroux lattice models. We then discuss parameter estimation and spatial prediction for the regularized Boltzmann-Gibbs covariance model in two dimensions. We recommend using the pairwise difference likelihood that combines satisfactory estimation performance and good scalability with many observation points. The predictive performance of the regularized covariance function is assessed by means of cross-validation statistics. Irregularly-spaced samples from the Walker Lake dataset are used, and spatial prediction is conducted by means of ordinary kriging. The regularized Boltzmann-Gibbs covariance yields improved predictive performance compared to the exponential covariance model.

The efficiency of a Markov sampler based on the underdamped Langevin diffusion is studied for high dimensional targets with convex and smooth potentials. We consider a classical second-order integrator which requires only one gradient computation per iteration. Contrary to previous works on similar samplers, a dimension-free contraction of Wasserstein distances and convergence rate for the total variance distance are proven for the discrete time chain itself. Non-asymptotic Wasserstein and total variation efficiency bounds and concentration inequalities are obtained for both the Metropolis adjusted and unadjusted chains. In particular, for the unadjusted chain, in terms of the dimension d and the desired accuracy ε, the Wasserstein efficiency bounds are of order $\sqrt{d}∕\mathit{\epsilon }$ in the general case, $\sqrt{d∕\mathit{\epsilon }}$ if the Hessian of the potential is Lipschitz, and $d^{1∕4}∕\sqrt{\mathit{\epsilon }}$ in the case of a separable target, in accordance with known results for other kinetic Langevin or HMC schemes.

In this work, we incorporate matrix projections into the reduced rank regression method, and then develop reduced rank regression estimators based on random projection and orthogonal projection in high-dimensional multivariate linear regression model. We propose a consistent estimator of the rank of the coefficient matrix and achieve prediction performance bounds for the proposed estimators based on mean squared errors. Finally, some simulation studies and a real data analysis are carried out to demonstrate that the proposed methods possess good stability, prediction performance and rank consistency compared to some other existing methods.

This manuscript presents an approach to perform generalized linear regression with multiple high dimensional covariance matrices as the outcome. In many areas of study, such as resting-state functional magnetic resonance imaging (fMRI) studies, this type of regression can be utilized to characterize variation in the covariance matrices across units. Model parameters are estimated by maximizing a likelihood formulation of a generalized linear model, conditioning on a well-conditioned linear shrinkage estimator for multiple covariance matrices, where the shrinkage coefficients are proposed to be shared across matrices. Theoretical studies demonstrate that the proposed covariance matrix estimator is optimal achieving the uniformly minimum quadratic loss asymptotically among all linear combinations of the identity matrix and the sample covariance matrix. Under certain regularity conditions, the proposed estimator of the model parameters is consistent. The superior performance of the proposed approach over existing methods is illustrated through simulation studies. Implemented to a resting-state fMRI study acquired from the Alzheimer’s Disease Neuroimaging Initiative, the proposed approach identified a brain network within which functional connectivity is significantly associated with Apolipoprotein E ε4, a strong genetic marker for Alzheimer’s disease.

The Gaussian model equips strong properties that facilitate studying and interpreting graphical models. Specifically it reduces conditional independence and the study of positive association to determining partial correlations and their signs. When Gaussianity does not hold partial correlation graphs are a useful relaxation of graphical models, but it is not clear what information they contain (besides the obvious lack of linear association). We study elliptical and transelliptical distributions as middle-ground between the Gaussian and other families that are more flexible but either do not embed strong properties or do not lead to simple interpretation. We characterize the meaning of zero partial correlations in elliptical and elliptical copula models and show that they retain much of the dependence structure from the Gaussian case. Regarding positive dependence, we prove impossibility results to learn certain positive (trans)elliptical graphical models, including that an elliptical distribution that is multivariate totally positive of order two for all dimensions must be essentially Gaussian. We then show how to interpret positive partial correlations as a relaxation, and obtain important properties related to faithfulness and Simpson’s paradox. We illustrate the transelliptical model potential to study tail dependence in S&P500 data, and of positivity to improve regularized inference.

We propose a difference-based nonparametric methodology for the estimation and inference of the time-varying auto-covariance functions of a locally stationary time series when it is contaminated by a complex trend with both abrupt and smooth changes. Simultaneous confidence bands (SCB) with asymptotically correct coverage probabilities are constructed for the auto-covariance functions under complex trend. A simulation-assisted bootstrapping method is proposed for the practical construction of the SCB. Detailed simulation and a real data example round out our presentation.

Suppose X is a multivariate diffusion process that is observed discretely in time. At each observation time, a transformation of the state of the process is observed with noise. The smoothing problem consists of recovering the path of the process, consistent with the observations. We derive a novel Markov Chain Monte Carlo algorithm to sample from the exact smoothing distribution. The resulting algorithm is called the Backward Filtering Forward Guiding (BFFG) algorithm. We extend the algorithm to include parameter estimation. The proposed method relies on guided proposals introduced in [53]. We illustrate its efficiency in a number of challenging problems.

We investigate estimation in Poisson regression model when the count response is right-censored and the censoring indicators are missing at random. We propose several estimators based on the regression calibration, multiple imputation and augmented inverse probability weighting methods. Under appropriate regularity conditions, we prove the consistency of our estimators and we derive their asymptotic distributions. Simulation experiments are carried out to investigate the finite sample behaviour and relative performance of the proposed estimates. These estimates are illustrated on a real data set.

Hypothesis testing of structure in covariance matrices is of significant importance, but faces great challenges in high-dimensional settings. Although consistent frequentist one-sample covariance tests have been proposed, there is a lack of simple, computationally scalable, and theoretically sound Bayesian testing methods for large covariance matrices. Motivated by this gap and by the need for tests that are powerful against sparse alternatives, we propose a novel testing framework based on the maximum pairwise Bayes factor. Our initial focus is on one-sample covariance testing; the proposed test can optimally distinguish null and alternative hypotheses in a frequentist asymptotic sense. We then propose diagonal tests and a scalable covariance graph selection procedure that are shown to be consistent. A simulation study evaluates the proposed approach relative to competitors. We illustrate advantages of our graph selection method on a gene expression data set.

Envelope method was recently proposed as a method to reduce the dimension of responses in multivariate regressions. However, when there exists missing data, the envelope method using the complete case observations may lead to biased and inefficient results. In this paper, we generalize the envelope estimation when the predictors and/or the responses are missing at random. Specifically, we incorporate the envelope structure in the expectation-maximization (EM) algorithm. As the parameters under the envelope method are not pointwise identifiable, the EM algorithm for the envelope method was not straightforward and requires a special decomposition. Our method is guaranteed to be more efficient, or at least as efficient as, the standard EM algorithm. Moreover, our method has the potential to outperform the full data MLE. We give asymptotic properties of our method under both normal and non-normal cases. The efficiency gain over the standard EM is confirmed in simulation studies and in an application to the Chronic Renal Insufficiency Cohort (CRIC) study.

We study approximation theorems for the Euler characteristic of the Vietoris-Rips and Čech filtration. The filtration is obtained from a Poisson or binomial sampling scheme in the critical regime. We apply our results to the smooth bootstrap of the Euler characteristic and determine its rate of convergence in the Kantorovich-Wasserstein distance and in the Kolmogorov distance.

We consider the problem of estimating the parameters a Gaussian Mixture Model with K components of known weights, all with an identity covariance matrix. We make two contributions. First, at the population level, we present a sharper analysis of the local convergence of EM and gradient EM, compared to previous works. Assuming a separation of $\mathrm{\Omega }(\sqrt{logK})$, we prove convergence of both methods to the global optima from an initialization region larger than those of previous works. Specifically, the initial guess of each component can be as far as (almost) half its distance to the nearest Gaussian. This is essentially the largest possible contraction region. Our second contribution are improved sample size requirements for accurate estimation by EM and gradient EM. In previous works, the required number of samples had a quadratic dependence on the maximal separation between the K components, and the resulting error estimate increased linearly with this maximal separation. In this manuscript we show that both quantities depend only logarithmically on the maximal separation.

Estimation of sparse, high-dimensional precision matrices is an important and challenging problem. Existing methods all assume that observations can be made precisely but, in practice, this often is not the case; for example, the instruments used to measure the response may have limited precision. The present paper incorporates measurement error in the context of estimating a sparse, high-dimensional precision matrix. In particular, for a Gaussian graphical model with data corrupted by Gaussian measurement error with unknown variance, we establish a general result which gives sufficient conditions under which the posterior contraction rates that hold in the no-measurement-error case carry over to the measurement-error case. Interestingly, this result does not require that the measurement error variance be small. We apply our general result to several cases with well-known prior distributions for sparse precision matrices and also to a case with a newly-constructed prior for precision matrices with a sparse factor-loading form. Two different simulation studies highlight the empirical benefits of accounting for the measurement error as opposed to ignoring it, even when that measurement error is relatively small.

We investigate the mean squared error and the asymptotic normality for a class of recursive kernel estimators based on a sample of spatially dependent observations. Our main result provides sufficient conditions for a spatial version of a recursive estimator introduced by Hall and Patil (1994) to satisfy a central limit theorem. The results are stated for strongly mixing random fields in the sense of Rosenblatt (1956) and for weakly dependent random fields in the sense of Wu (2005).

We consider the problem of constructing confidence intervals for the median of a response $Y\in \mathbb{R}$ conditional on features $X\in {\mathbb{R}}^d$ in a situation where we are not willing to make any assumption whatsoever on the underlying distribution of the data $(X,Y)$. We propose a method based upon ideas from conformal prediction and establish a theoretical guarantee of coverage while also going over particular distributions where its performance is sharp. Additionally, we prove an equivalence between confidence intervals for the conditional median and confidence intervals for the response variable, resulting in a lower bound on the length of any possible conditional median confidence interval. This lower bound is independent of sample size and holds for all distributions with no point masses.

The $\mathrm{ROC}$ curve is the gold standard for measuring the performance of a test/scoring statistic regarding its capacity to discriminate between two statistical populations in a wide variety of applications, ranging from anomaly detection in signal processing to information retrieval, through medical diagnosis. Most practical performance measures used in scoring/ranking applications such as the $\mathrm{AUC}$, the local $\mathrm{AUC}$, the p-norm push, the DCG and others, can be viewed as summaries of the $\mathrm{ROC}$ curve. In this paper, the fact that most of these empirical criteria can be expressed as two-sample linear rank statistics is highlighted and concentration inequalities for collections of such random variables, referred to as two-sample rank processes here, are proved, when indexed by VC classes of scoring functions. Based on these nonasymptotic bounds, the generalization capacity of empirical maximizers of a wide class of ranking performance criteria is next investigated from a theoretical perspective. It is also supported by empirical evidence through convincing numerical experiments.

In this work, we study a new recursive stochastic algorithm for the joint estimation of quantile and superquantile of an unknown distribution. The novelty of this algorithm is to use the Cesaro averaging of the quantile estimation inside the recursive approximation of the superquantile. We provide some sharp non-asymptotic bounds on the quadratic risk of the superquantile estimator for different step size sequences. We also prove new non-asymptotic $L^p$-controls on the Robbins Monro algorithm for quantile estimation and its averaged version. Finally, we derive a central limit theorem of our joint procedure using the diffusion approximation point of view hidden behind our stochastic algorithm.

Recursive max-linear vectors model causal dependence between node variables by a structural equation model, expressing each node variable as a max-linear function of its parental nodes in a directed acyclic graph (DAG) and some exogenous innovation. For such a model, there exists a unique minimum DAG, represented by the Kleene star matrix of its edge weight matrix, which identifies the model and can be estimated. For a more realistic statistical modeling we introduce some random observational noise. A probabilistic analysis of this new noisy model reveals that the unique minimum DAG representing the distribution of the non-noisy model remains unchanged and identifiable. Moreover, the distribution of the minimum ratio estimators of the model parameters at their left limits are completely determined by the distribution of the noise variables up to a positive constant. Under a regular variation condition on the noise variables we prove that the estimated Kleene star matrix converges to a matrix of independent Weibull entries after proper centering and scaling.

In this paper, we study the nonparametric linear model, when the error process is a dependent Gaussian process. We focus on the estimation of the mean vector via a model selection approach. We first give the general theoretical form of the penalty function, ensuring that the penalized estimator among a collection of models satisfies an oracle inequality. Then we derive a penalty shape involving the spectral radius of the covariance matrix of the errors, which can be chosen proportional to the dimension when the error process is stationary and short range dependent. However, this penalty can be too rough in some cases, in particular when the error process is long range dependent. In a second part, we focus on the fixed-design regression model assuming that the error process is a stationary Gaussian process. We propose a model selection procedure in order to estimate the mean function via piecewise polynomials on a regular partition, when the error process is either short range dependent, long range dependent or anti-persistent. We present different kinds of penalties, depending on the memory of the process. For each case, an adaptive estimator is built, and the rates of convergence are computed. Thanks to several sets of simulations, we study the performance of these different penalties for all types of errors (short memory, long memory and anti-persistent errors). Finally, we give an application of our method to the well-known Nile data, which clearly shows that the type of dependence of the error process must be taken into account.

A sample covariance matrix S of completely observed data is the key statistic in a large variety of multivariate statistical procedures, such as structured covariance/precision matrix estimation, principal component analysis, and testing of equality of mean vectors. However, when the data are partially observed, the sample covariance matrix from the available data is biased and does not provide valid multivariate procedures. To correct the bias, a simple adjustment method called inverse probability weighting (IPW) has been used in previous research, yielding the IPW estimator. The estimator can play the role of S in the missing data context, thus replacing S in off-the-shelf multivariate procedures such as the graphical lasso algorithm. However, theoretical properties (e.g. concentration) of the IPW estimator have been only established in earlier work under very simple missing structures; every variable of each sample is independently subject to missingness with equal probability. We investigate the deviation of the IPW estimator when observations are partially observed under general missing dependency. We prove the optimal convergence rate $O_p(\sqrt{logp∕n})$ of the IPW estimator based on the element-wise maximum norm, even when two unrealistic assumptions (known mean and/or missing probabilities) frequently assumed to be known in the past work are relaxed. The optimal rate is especially crucial in estimating a precision matrix, because of the “meta-theorem” [26] that claims the rate of the IPW estimator governs that of the resulting precision matrix estimator. In the simulation study, we discuss one of practically important issues, non-positive semi-definiteness of the IPW estimator, and compare the estimator with imputation methods.