Simulation from Wishart Distributions with Eigenvalue Constraints

Abstract

This article provides an efficient algorithm for generating a random matrix according to a Wishart distribution, but with eigenvalues constrained to be less than a given vector of positive values. The procedure of Odell and Feiveson provides a guide, but the modifications here ensure that the diagonal elements of a candidate matrix are less than the corresponding elements of the constraint vector, thus greatly improving the chances that the matrix will be acceptable. The Normal hierarchical model with vector outcomes and the multivariate random effects model provide motivating applications.     

Error-free matrix symmetrizers and equivalent symmetric matrices

A symmetrizer of the matrix A is a symmetric solution X that satisfies the matrix equation XA=AX. An exact matrix symmetrizer is computed by obtaining a general algorithm and superimposing a modified multiple modulus residue arithmetic on this algorithm. A procedure based on computing a symmetrizer to obtain a symmetric matrix, called here an equivalent symmetric matrix, whose eigenvalues are the same as those of a given real nonsymmetric matrix is presented.Supported by CSIR.     

On the Spectral Distribution of Gaussian Random Matrices

We consider the empirical spectral distribution （ESD） of a random matrix from the Gaussian Unitary Ensemble. Based on the Plancherel-Rotaeh approximation formula for Hermite polynomials, we prove that the expected empirical spectral distribution converges at the rate of O（n^-1） to the Wigner distribution function uniformly on every compact intervals [u,v] within the limiting support （-1, 1）. Furthermore, the variance of the ESD for such an interval is proved to be （πn）^-2 logn asymptotically which surprisingly enough, does not depend on the details （e.g. length or location） of the interval, This property allows us to determine completely the covariance function between the values of the ESD on two intervals.     

Iterative Weighted Semiparametric Least Squares Estimation in Repeated Measurement Partially Linear Regression Models

Consider a repeated measurement partially linear regression model with an unknown vector parameter β, an unknown function g(.), and unknown heteroscedastic error variances. In order to improve the semiparametric generalized least squares estimator (SGLSE) of β, we propose an iterative weighted semiparametric least squares estimator (IWSLSE) and show that it improves upon the SGLSE in terms of asymptotic covariance matrix. An adaptive procedure is given to determine the number of iterations. We also show that when the number of replicates is less than or equal to two, the IWSLSE can not improve upon the SGLSE. These results are generalizations of those in [2] to the case of semiparametric regressions.     

Estimating tail probabilities of heavy tailed distributions with asymptotically zero relative error

Efficient estimation of tail probabilities involving heavy tailed random variables is amongst the most challenging problems in Monte-Carlo simulation. In the last few years, applied probabilists have achieved considerable success in developing efficient algorithms for some such simple but fundamental tail probabilities. Usually, unbiased importance sampling estimators of such tail probabilities are developed and it is proved that these estimators are asymptotically efficient or even possess the desirable bounded relative error property. In this paper, as an illustration, we consider a simple tail probability involving geometric sums of heavy tailed random variables. This is useful in estimating the probability of large delays in M/G/1 queues. In this setting we develop an unbiased estimator whose relative error decreases to zero asymptotically. The key idea is to decompose the probability of interest into a known dominant component and an unknown small component. Simulation then focuses on estimating the latter ‘residual’ probability. Here we show that the existing conditioning methods or importance sampling methods are not effective in estimating the residual probability while an appropriate combination of the two estimates it with bounded relative error. As a further illustration of the proposed ideas, we apply them to develop an estimator for the probability of large delays in stochastic activity networks that has an asymptotically zero relative error.      

Numerical analysis of(s, S) inventory systems with repeated attempts

This paper deals with a continuous review (s,S) inventory system where arriving demands finding the system out of stock, leave the service area and repeat their request after some random time. This assumption introduces a natural alternative to classical approaches based either on lost demand models or on backlogged models. The stochastic model formulation is based on a bidimensional Markov process which is numerically solved to investigate the essential operating characteristics of the system. An optimal design problem is also considered. AMS subject classification: 90B05 90B22     

A posteriori error estimate for computing tr(f(A)) by using the Lanczos method

An outstanding problem when computing a function of a matrix, f(A), by using a Krylov method is to accurately estimate errors when convergence is slow. Apart from the case of the exponential function that has been extensively studied in the past, there are no well‐established solutions to the problem. Often, the quantity of interest in applications is not the matrix f(A) itself but rather the matrix–vector products or bilinear forms. When the computation related to f(A) is a building block of a larger problem (e.g., approximately computing its trace), a consequence of the lack of reliable error estimates is that the accuracy of the computed result is unknown. In this paper, we consider the problem of computing tr(f(A)) for a symmetric positive‐definite matrix A by using the Lanczos method and make two contributions: (a) an error estimate for the bilinear form associated with f(A) and (b) an error estimate for the trace of f(A). We demonstrate the practical usefulness of these estimates for large matrices and, in particular, show that the trace error estimate is indicative of the number of accurate digits. As an application, we compute the log determinant of a covariance matrix in Gaussian process analysis and underline the importance of error tolerance as a stopping criterion as a means of bounding the number of Lanczos steps to achieve a desired accuracy.     

Adaptive Monte Carlo methods for matrix equations with applications

This paper discusses empirical studies with both the adaptive correlated sequential sampling method and the adaptive importance sampling method which can be used in solving matrix and integral equations. Both methods achieve geometric convergence (provided the number of random walks per stage is large enough) in the sense: e_ν≤cλ^ν, where e_ν is the error at stage ν, λ∈(0,1) is a constant, c>0 is also a constant. Thus, both methods converge much faster than the conventional Monte Carlo method. Our extensive numerical test results show that the adaptive importance sampling method converges faster than the adaptive correlated sequential sampling method, even with many fewer random walks per stage for the same problem. The methods can be applied to problems involving large scale matrix equations with non-sparse coefficient matrices. We also provide an application of the adaptive importance sampling method to the numerical solution of integral equations, where the integral equations are converted into matrix equations (with order up to 8192×8192) after discretization. By using Niederreiter’s sequence, instead of a pseudo-random sequence when generating the nodal point set used in discretizing the phase space Γ, we find that the average absolute errors or relative errors at nodal points can be reduced by a factor of more than one hundred.     

Estimation of the Failure Rate in a Partially Accelerated Life Test: A Sequential Approach

Abstract

Consider the situation in which a group of units are put on a partially accelerated life test. It is assumed that the lifelengths of the units are independent and exponentially distributed random variables with common failure rate θ, and that θ is the value of a random variable having a gamma distribution. A two‐stage sequential procedure for estimating θ under the squared error loss is proposed. In the first stage, the units are put on the test under normal stress up to time t, where t is determined as a stopping time that minimizes the expected loss plus cost of running the test. In the second stage, the stress is raised to a higher level for those units that did not fail by time t and held constant until they all fail. The accumulated data are then used to estimate θ with the Bayes estimator.     

Pivoting strategy for rank-one modification ofLDM t-like factorization

The rank-one modification algorithm of theLDM ^t factorization was given by Bennett [1]. His method, however, could break down even when the matrix is nonsingular and well-conditioned. We introduce a pivoting strategy for avoiding possible break-down as well as for suppressing error growth in the modification process. The method is based on a symbolic formula of the rank-one modification of the factorization of a possibly singular nonsymmetric matrix. A new symbolic formula is also obtained for the inverses of the factor matrices. Repeated application of our method produces theLDM ^t-like product form factorization of a matrix. A numerical example is given to illustrate our pivoting method. An incomplete factorization algorithm is also introduced for updating positive definite matrix useful in quasi-Newton methods, in which the Fletcher and Powell algorithm [2] and the Gill, Murray and Saunders algorithm [4] are usually used.This paper is presented at the Japan SIAM Annual Meeting held at University of Tokyo, Japan, October 7–9, 1991.     

Generation of “Similar” Images from a Given Discrete Image

Abstract

A discrete image of several colors is viewed as a discrete random field obtained by clipping or quantizing a Gaussian random field at several levels. Given a discrete image, parameters of the unobserved original Gaussian random field are estimated. Discrete images, statistically similar to the original image, are then obtained by generating different realizations of the Gaussian field and clipping them. To overcome the computational difficulties, the block Toeplitz covariance matrix of the Gaussian field is embedded into a block circulant matrix which is diagonalized by the fast Fourier transform. The Gibbs sampler is used to apply the stochastic EM algorithm for the estimation of the field's parameters.     

Low-lying zeros of number field L-functions

TextOne of the most important statistics in studying the zeros of L-functions is the 1-level density, which measures the concentration of zeros near the central point. Fouvry and Iwaniec (2003) [FI] proved that the 1-level density for L-functions attached to imaginary quadratic fields agrees with results predicted by random matrix theory. In this paper, we show a similar agreement with random matrix theory occurring in more general sequences of number fields. We first show that the main term agrees with random matrix theory, and similar to all other families studied to date, is independent of the arithmetic of the fields. We then derive the first lower order term of the 1-level density, and see the arithmetic enter.VideoFor a video summary of this paper, please click here or visit http://www.youtube.com/watch?v=zpb-gu3G8i0.     

Random orthogonal matrix simulation

This paper introduces a method for simulating multivariate samples that have exact means, covariances, skewness and kurtosis. We introduce a new class of rectangular orthogonal matrix which is fundamental to the methodology and we call these matrices L matrices. They may be deterministic, parametric or data specific in nature. The target moments determine the L matrix then infinitely many random samples with the same exact moments may be generated by multiplying the L matrix by arbitrary random orthogonal matrices. This methodology is thus termed “ROM simulation”. Considering certain elementary types of random orthogonal matrices we demonstrate that they generate samples with different characteristics. ROM simulation has applications to many problems that are resolved using standard Monte Carlo methods. But no parametric assumptions are required (unless parametric L matrices are used) so there is no sampling error caused by the discrete approximation of a continuous distribution, which is a major source of error in standard Monte Carlo simulations. For illustration, we apply ROM simulation to determine the value-at-risk of a stock portfolio.     

Exponential Mixture Representation of Geometric Stable Distributions

We show that every strictly geometric stable (GS) random variable can be represented as a product of an exponentially distributed random variable and an independent random variable with an explicit density and distribution function. An immediate application of the representation is a straightforward simulation method of GS random variables. Our result generalizes previous representations for the special cases of Mittag-Leffler and symmetric Linnik distributions.     

A reinforcement-depletion urn model: A contiguity case

Summary An urn contains balls ofs different colors. The problem of the reinforcement of a specified color and random depletion of balls has been considered by Bernard (1977,Bull. Math. Biol.,39, 463–470) and Shenton (1981,Bull. Math. Biol.,43, 327–340), (1983,Bull. Math. Biol.,45, 1–9). Here we consider a special relation between a reinforcement and depletion, leading to a hypergeometric distribution. Research sponsored in part by the Applied Mathematical Sciences Research Program, Office of Energy Research, U.S. Department of Energy under contract DE-AC05-840R21400 with the Martin Marietta Energy Systems, Inc.     

The Limiting Spectra of Girko’s Block-Matrix

To analyze the limiting spectral distribution of some random block-matrices, Girko (Random Oper. Stoch. Equ. 8(2), 189–194, 2000) uses a system of canonical equations from (An Introduction to Statistical Analysis of Random Arrays. VSP, Utrecht, 1998). In this paper, we use the method of moments to give an integral form for the almost sure limiting spectral distribution of such matrices.     

A shorter proof of Kanter’s Bessel function concentration bound

We give a shorter proof of Kanter’s (J. Multivariate Anal. 6, 222–236, 1976) sharp Bessel function bound for concentrations of sums of independent symmetric random vectors. We provide sharp upper bounds for the sum of modified Bessel functions I₀(x) + I₁(x), which might be of independent interest. Corollaries improve concentration or smoothness bounds for sums of independent random variables due to Čekanavičius & Roos (Lith. Math. J. 46, 54–91, 2006); Roos (Bernoulli, 11, 533–557, 2005), Barbour & Xia (ESAIM Probab. Stat. 3, 131–150, 1999), and Le Cam (Asymptotic Methods in Statistical Decision Theory. Springer, Berlin Heidelberg New York, 1986).      

Second-order nonlinear least squares estimation

The ordinary least squares estimation is based on minimization of the squared distance of the response variable to its conditional mean given the predictor variable. We extend this method by including in the criterion function the distance of the squared response variable to its second conditional moment. It is shown that this “second-order” least squares estimator is asymptotically more efficient than the ordinary least squares estimator if the third moment of the random error is nonzero, and both estimators have the same asymptotic covariance matrix if the error distribution is symmetric. Simulation studies show that the variance reduction of the new estimator can be as high as 50% for sample sizes lower than 100. As a by-product, the joint asymptotic covariance matrix of the ordinary least squares estimators for the regression parameter and for the random error variance is also derived, which is only available in the literature for very special cases, e.g. that random error has a normal distribution. The results apply to both linear and nonlinear regression models, where the random error distributions are not necessarily known.     

Steepest Descent Method with Random Step Lengths

The paper studies the steepest descent method applied to the minimization of a twice continuously differentiable function. Under certain conditions, the random choice of the step length parameter, independent of the actual iteration, generates a process that is almost surely R-convergent for quadratic functions. The convergence properties of this random procedure are characterized based on the mean value function related to the distribution of the step length parameter. The distribution of the random step length, which guarantees the maximum asymptotic convergence rate independent of the detailed properties of the Hessian matrix of the minimized function, is found, and its uniqueness is proved. The asymptotic convergence rate of this optimally created random procedure is equal to the convergence rate of the Chebyshev polynomials method. Under practical conditions, the efficiency of the suggested random steepest descent method is degraded by numeric noise, particularly for ill-conditioned problems; furthermore, the asymptotic convergence rate is not achieved due to the finiteness of the realized calculations. The suggested random procedure is also applied to the minimization of a general non-quadratic function. An algorithm needed to estimate relevant bounds for the Hessian matrix spectrum is created. In certain cases, the random procedure may surpass the conjugate gradient method. Interesting results are achieved when minimizing functions having a large number of local minima. Preliminary results of numerical experiments show that some modifications of the presented basic method may significantly improve its properties.     

Maxima and minima of homogeneous Gaussian random fields over continuous time and uniform grids

ABSTRACT

In this paper, for centred homogeneous Gaussian random fields the joint limiting distributions of normalized maxima and minima over continuous time and uniform grids are investigated. It is shown that maxima and minima are asymptotic dependent for strongly dependent homogeneous Gaussian random field with the choice of sparse grid, Pickands' grid or dense grid, while for the weakly dependent Gaussian random field maxima and minima are asymptotically independent.