You do not have any saved searches

For selected items:
  • Asymptotic Global Behavior for Stochastic Approximation and Diffusions with Slowly Decreasing Noise Effects: Global Minimization via Monte Carlo

    Abstract

    The asymptotic behavior of the systems $X_{n + 1} = X_n + a_n b( {X_n ,\xi _n } ) + a_n \sigma ( X_n )\psi_n $ and $dy = \bar b( y )dt + \sqrt {a( t )} \sigma ( y )dw$ is studied, where $\{ {\psi _n } \}$ is i.i.d. Gaussian, $\{ \xi _n \}$ is a (correlated) bounded sequence of random variables and $a_n \approx A_0/\log (A_1 + n )$. Without $\{ \xi _n \}$, such algorithms are versions of the “simulated annealing” method for global optimization. When the objective function values can only be sampled via Monte Carlo, the discrete algorithm is a combination of stochastic approximation and simulated annealing. Our forms are appropriate. The $\{ \psi _n \}$ are the “annealing” variables, and $\{ \xi _n \}$ is the sampling noise. For large $A_0 $, a full asymptotic analysis is presented, via the theory of large deviations: Mean escape time (after arbitrary time n) from neighborhoods of stable sets of the algorithm, mean transition times (after arbitrary time n) from a neighborhood of one stable set to another, approximate asymptotic invariant measures, and location of the values of $\{ X_n \}$ or $y( \cdot )$, the case where $Eb( x,\xi ) = \bar b( x )$ is the (negative) of a gradient of a function $\bar B( x )$, and application to global function minimization via Monte Carlo methods.
  • A Simple Unpredictable Pseudo-Random Number Generator

    Abstract

    Two closely-related pseudo-random sequence generators are presented: The ${1 / P}$generator, with input P a prime, outputs the quotient digits obtained on dividing 1 by P. The $x^2 \bmod N$generator with inputs N, $x_0 $ (where $N = P \cdot Q$ is a product of distinct primes, each congruent to 3 mod 4, and $x_0 $ is a quadratic residue $\bmod N$), outputs $b_0 b_1 b_2 \cdots $ where $b_i = {\operatorname{parity}}(x_i )$ and $x_{i + 1} = x_i^2 \bmod N$.From short seeds each generator efficiently produces long well-distributed sequences. Moreover, both generators have computationally hard problems at their core. The first generator’s sequences, however, are completely predictable (from any small segment of $2|P| + 1$ consecutive digits one can infer the “seed,” P, and continue the sequence backwards and forwards), whereas the second, under a certain intractability assumption, is unpredictable in a precise sense. The second generator has additional interesting properties: from knowledge of $x_0 $ and N but notP or Q, one can generate the sequence forwards, but, under the above-mentioned intractability assumption, one can not generate the sequence backwards. From the additional knowledge of P and Q, one can generate the sequence backwards; one can even “jump” about from any point in the sequence to any other. Because of these properties, the $x^2 \bmod N$generator promises many interesting applications, e.g., to public-key cryptography. To use these generators in practice, an analysis is needed of various properties of these sequences such as their periods. This analysis is begun here.
  • From Micro to Macro Dynamics via a New Closure Approximation to the FENE Model of Polymeric Fluids

    Abstract

    We present a new closure approximation needed for deriving effective macroscopic moment equations from the microscopic finite-extensible-nonlinear-elastic kinetic theory modeling viscoelastic polymeric fluids. The closure is based on restricting the otherwise general probability distribution functions (PDFs) to a class of smooth distributions motivated by perturbing the equilibrium PDF\@. The simplified system coupling the moment equations and the Navier--Stokes equations still possesses an energy law analogous to the original micro-macro system. Some theoretical analysis and numerical experiments are presented to ensure the validity of the moment-closure system, and to illustrate the excellent agreement of the simplified model with the original system solved using a Monte Carlo approach, for a certain regime of physical parameters.
  • Solving Linear Equations by Monte Carlo Simulation

    Abstract

    A new Monte Carlo estimator for solving the matrix equation x=Hx+b is presented, and theoretical results comparing this estimator with the traditional terminal and collision estimators are given. We then make a detailed investigation of the average complexity of the Monte Carlo estimators when several popular random variable generation techniques are used, and we show that the average complexity can be as low as $C\Vert z\Vert _{1}n+N$, where n is the number of random walks generated, N is the size of the matrix H, z is the solution of an associated matrix equation, and C is a small constant. As a consequence of the complexity results, we observe how scaling of matrices, a well-known technique in deterministic methods, can increase the efficiency of the Monte Carlo method.One advantage of the Monte Carlo method is the efficiency at which it can be parallelized. The algorithms we discuss can provide fast and approximate solutions to systems of linear equations in massively parallel computing environments. Surprisingly, sequential (or adaptive) Monte Carlo methods can even be competitive in single-processorcomputing environments. We present numerical results and compare our Monte Carlo algorithms with Krylov subspace methods for some test matrices.
  • The Linear Process Deferment Algorithm: A new technique for solving population balance equations

    Abstract

    In this paper a new stochastic algorithm for the solution of population balance equations is presented. The population balance equations have the form of extended Smoluchowski equations which include linear and source terms. The new algorithm, called the linear process deferment algorithm (LPDA), is used for solving a detailed model describing the formation of soot in premixed laminar flames. A measure theoretic formulation of a stochastic jump process is developed and the corresponding generator presented. The numerical properties of the algorithm are studied in detail and compared to the direct simulation algorithm and various splitting techniques. LPDA is designed for all kinds of population balance problems where nonlinear processes cannot be neglected but are dominated in rate by linear ones. In those cases the LPDA is seen to reduce run times for a population balance simulation by a factor of up to 1000 with a negligible loss of accuracy.
  • A Hierarchical Multiscale Approach to Protein Structure Prediction: Production of Low‐Resolution Packing Arrangements of Helices and Refinement of the Best Models with a United‐Residue Force Field

    Abstract

    A hierarchical, two‐stage approach to ab initio protein structure prediction is presented and applied to four α‐helical proteins. In the first stage, a bank of low‐resolution models is generated using a highly simplified protein representation and energy function, coupled with a Conformation‐Family Monte Carlo (CFMC) search for the energy minimum. For helical proteins, this procedure (referred to as REPACK) produces a set of plausible packed arrangements of the helices, given their positions in the amino acid sequence. Secondary structure prediction methods such as JPRED can be used to provide the secondary structure assignment. In the second stage, these packing arrangements are used as starting points for a new search method (Local Search), based on the Monte Carlo‐with‐Minimization (MCM) algorithm and a united‐residue (UNRES) energy function. The focus of the Local Search is mainly on improving loop conformations and side‐chain positions, with minor modifications to the overall packing of the helices. By reducing the size of the conformational space that must be sampled with the UNRES energy function, which is much more expensive to compute than the REPACK energy function, this prediction scheme can be applied to much larger proteins than were tractable in the past with other UNRES‐based search methods. It was applied successfully to a 224–residue protein (target T0198, PDB code 1SUM) in the sixth community‐wide blind‐prediction experiment on the Critical Assessment of techniques for protein Structure Prediction (CASP6).
  • Quantization Based Filtering Method Using First Order Approximation

    Abstract

    The quantization based filtering method (see [G. Pagès and H. Pham, Bernoulli, 11 (2005), pp. 893–932; G. Pagès, H. Pham, and J. Printems, Optimal quantization methods and applications to numerical problems in finance, in Handbook of Computational and Numerical Methods in Finance, S. T. Rachev, ed., Birkhäuser, Boston, 2004, pp. 253–297]) is a grid based approximation method to solve nonlinear filtering problems with discrete time observations. It relies on off-line preprocessing of some signal grids in order to construct fast recursive schemes for filter approximation. We give here an improvement of this method by taking advantage of the stationary quantizer property. The key ingredient is the use of vanishing correction terms to describe schemes based on piecewise linear approximations. Convergence results are given and numerical results are presented for the particular cases of linear Gaussian model and stochastic volatility models.
  • Dynamical Spatial Warping: A Novel Method for the Conformational Sampling of Biophysical Structure

    Abstract

    The difficulties encountered in sampling of systems with rough energy landscapes using present methodology significantly limit the impact of simulation on molecular biology, in particular protein folding and design. Here, we present a major methodological development based on a promising new technique, the reference potential spatial warping algorithm (REPSWA) [Z. Zhu et al., Phys. Rev. Lett., 88 (2002), pp. 100201–100204], and present applications to several realistic systems. REPSWA works by introducing a variable transformation in the classical partition function that reduces the volume of phase space associated with a priori known barrier regions while increasing that associated with attractive basins. In this way, the partition function is preserved so that enhanced sampling is achieved without the need for reweighting phase-space averages. Here, a new class of transformations, designed to overcome the barriers induced by intermolecular/nonbonded interactions, whose locations are not known a priori, is introduced. The new transformations are designed to work in synergy with transformations originally introduced for overcoming intramolecular barriers. The new transformation adapts to the fluctuating local environment and is able to handle barriers that arise “on the fly.” Thus, the new method is referred to as dynamic contact REPSWA (DC-REPSWA). In addition, combining hybrid Monte Carlo (HMC) with DC-REPSWA allows more aggressive sampling to take place. The combined DC-REPSWA-HMC method and its variants are shown to substantially enhance conformational sampling in long molecular chains composed of interacting single beads and beads with branches. The latter topologies characterize the united residue and united side chain representation of protein structures.
  • Fast Monte Carlo Simulation Methods for Biological Reaction-Diffusion Systems in Solution and on Surfaces

    Abstract

    Many important physiological processes operate at time and space scales far beyond those accessible to atom-realistic simulations, and yet discrete stochastic rather than continuum methods may best represent finite numbers of molecules interacting in complex cellular spaces. We describe and validate new tools and algorithms developed for a new version of the MCell simulation program (MCell3), which supports generalized Monte Carlo modeling of diffusion and chemical reaction in solution, on surfaces representing membranes, and combinations thereof. A new syntax for describing the spatial directionality of surface reactions is introduced, along with optimizations and algorithms that can substantially reduce computational costs (e.g., event scheduling, variable time and space steps). Examples for simple reactions in simple spaces are validated by comparison to analytic solutions. Thus we show how spatially realistic Monte Carlo simulations of biological systems can be far more cost-effective than often is assumed, and provide a level of accuracy and insight beyond that of continuum methods.
  • A Sample Approximation Approach for Optimization with Probabilistic Constraints

    Abstract

    We study approximations of optimization problems with probabilistic constraints in which the original distribution of the underlying random vector is replaced with an empirical distribution obtained from a random sample. We show that such a sample approximation problem with a risk level larger than the required risk level will yield a lower bound to the true optimal value with probability approaching one exponentially fast. This leads to an a priori estimate of the sample size required to have high confidence that the sample approximation will yield a lower bound. We then provide conditions under which solving a sample approximation problem with a risk level smaller than the required risk level will yield feasible solutions to the original problem with high probability. Once again, we obtain a priori estimates on the sample size required to obtain high confidence that the sample approximation problem will yield a feasible solution to the original problem. Finally, we present numerical illustrations of how these results can be used to obtain feasible solutions and optimality bounds for optimization problems with probabilistic constraints.
  • A Quasi-Monte Carlo Approach to Particle Simulation of the Heat Equation

    Abstract

    The convergence of the Monte Carlo method for numerical integration can often be improved by replacing random numbers with more uniformly distributed numbers known as quasi-random. In this paper the convergence of Monte Carlo particle simulation is studied when these quasi-random sequences are used. For the one-dimensional heat equation discretized in both space and time, convergence is proved for a quasi-random simulation using reordering of the particles according to their position. Experimental results are presented for the spatially continuous heat equation in one and two dimensions. The results indicate that a significant improvement in both magnitude of error and convergence rate can be achieved over standard Monte Carlo simulations for certain low-dimensional problems.
  • Accurate Uncertainty Quantification Using Inaccurate Computational Models

    Abstract

    This paper proposes a novel uncertainty quantification framework for computationally demanding systems characterized by a large vector of non-Gaussian uncertainties. It combines state-of-the-art techniques in advanced Monte Carlo sampling with Bayesian formulations. The key departure from existing works is the use of inexpensive, approximate computational models in a rigorous manner. Such models can readily be derived by coarsening the discretization size in the solution of the governing PDEs, increasing the time step when integration of ODEs is performed, using fewer iterations if a nonlinear solver is employed or making use of lower-order models. It is shown that even in cases where the inaccurate models provide very poor approximations of the high-fidelity response, statistics of the latter can be quantified accurately with significant reductions in the computational effort. Multiple approximate models can be used, and rigorous confidence bounds of the estimates produced are provided at all stages.
  • A Fast, Easily Implemented Method for Sampling from Decreasing or Symmetric Unimodal Density Functions

    Abstract

    The fastest computer methods for sampling from a given density are those based on a mixture of a fast and slow part. This paper describes a new method of this type, suitable for any decreasing or symmetric unimodal density. It differs from others in that it is faster and more easily implemented, thereby providing a standard procedure for developing both the fast and the slow part for many given densities. It is called the ziggurat method, after the shape of a single, convenient density that provides for both the fast and the slow parts of the generating process. Examples are given for REXP and RNOR, subroutines that generate exponential and normal variates that, as assembler routines, are nearly twice as fast as the previous best assembler routines, and that, as Fortran subroutines, approach the limiting possible speed: the time for Fortran subroutine linkage conventions plus the time to generate one uniform variate.
  • Generating Correlation Matrices

    Abstract

    This paper describes a variety of methods for generating random correlation matrices, with emphasis on choice of random variables and distributions so as to provide matrices with given structure, expected values or eigenvalues.
  • A Monte Carlo Method for Scalar Reaction Diffusion Equations

    Abstract

    A probabilistic method is presented to solve reaction diffusion equations. A random walk is combined with creation and destruction of elements. The method is applied to Nagumo’s equation. Numerical results are given demonstrating convergence of the method. The stochastic process also gives a direct probability interpretation of the equation which may be useful for analysis.
  • Generation of Random Orthogonal Matrices

    Abstract

    In order to generate a random orthogonal matrix distributed according to Haar measure over the orthogonal group it is natural to start with a matrix of normal random variables and then factor it by the singular value decomposition. A more efficient method is obtained by using Householder transformations. We propose another alternative based on the product of ${{n(n - 1)}/2}$ orthogonal matrices, each of which represents an angle of rotation. Some numerical comparisons of alternative methods are made.
  • Sampling from Gauss Rules

    Abstract

    Approximating multidimensional integrals via product quadrature rules becomes increasingly intractable as the dimension increases. Hammersley [Ann. New York Acad. Sci., 86 (1960), pp. 844–874] suggested sampling from product quadrature rules and Tsuda [Numer. Math., 20 (1973), pp. 377–391] considered this method using Féjer rules. In this paper we consider this approach using Gauss rules. Results are obtained concerning the variance of this form of sampling relative to sampling from the continuous distributions represented by the weight functions. It is shown that this approach can lead to variance reduction and its use is discussed in several examples.
  • Discrete Event Simulations and Parallel Processing: Statistical Properties

    Abstract

    This paper addresses statistical issues that arise when discrete event, or Monte Carlo, simulations are run on parallel processing computers. In particular, the statistical properties of estimators obtained by running parallel independent replications on a multiple processor computing system are considered. Because of the effects of parallelism, care must be taken in order to obtain estimators with the proper statistical properties. A variety of estimators are considered. The convergence properties, including strong laws, central limit theorems and bias expansions of these estimators are derived. It is shown that some of the more obvious estimators are guaranteed to converge to the wrong quantity as the number of processors increases. Strong laws and central limit theorems for completion times of the estimators are also given. The application of results from reliability and scheduling theory yields bounds on expected completion times under a variety of distributional assumptions. Based on these results, it does not appear possible to obtain a strongly consistent estimate in finite expected time as the number of processors increases, unless the computational time to complete a single replication is bounded.
  • Monte Carlo, Control Variates, and Stochastic Ordering

    Abstract

    Using prior information about the stochastic orderings between a phenomenon of interest and other ancillary phenomena whose population parameters are known in a Monte Carlo experiment, this paper derives an unbiased control variate estimator and an exact $100 \times (1 - \alpha )$ percent confidence interval for a proportion and, more generally, for a distribution function for finite sample size K. These results substantially improve on the more traditional control variate variance-reduction method for which, in the case of finite K, neither an unbiased point estimator is generally available nor, except in the case of multinormal data, is an exact $100 \times (1 - \alpha )$ percent confidence interval available.The paper also derives upper bounds on the estimator’s variance and coefficient of variation and a lower bound on the achievable variance reduction in terms of the known parameters for the ancillary phenomena. Computing these bounds prior to sampling provides useful guidance rarely available with other variance reduction techniques. In the case of estimating a discrete distribution function, the paper also shows how to derive a variance-reduced estimator of the population mean from the estimated distribution function, and gives the variance of this sample mean to order ${1 / K}$. An example illustrates the technique for estimating the distribution of maximal flow in a flow-network problem with random arc capacities.
  • Comparing Averaged-Out Utilities of Probability Trees Having Random Parameters

    Abstract

    This paper develops a procedure to approximate the marginal distributions of each of the following performance measures defined on a forest of probability trees with random parameters: (a) the averaged-out utility of each individual tree in the forest, and (b) the difference between the averaged-out utilities of each pair of trees in the forest. The parameters of a probability tree are its branching probabilities and node utilities, and they may be specified as mutually independent random variables with given distributions; moreover, any of these parameters may be duplicated at several places in the forest to represent parameter dependencies within and between trees. The approximation procedure incorporates an efficient,numerically robust algorithm for computing up to the first four moments of the selected performance measure. The desired distribution is approximated by an Edgeworth expansion based on the computed moments. With respect to accuracy and computational complexity, the approximation procedure is compared to estimation techniques based on simulation and normal distribution theory. To illustrate how the analysis of probability trees can directly incorporate uncertainty or random variation in the values of tree parameters, the approximation procedure is applied to the comparison of alternative medical protocols for managing patients following myocardial infarction (heart attack). In addition to its applications in medical decision making, the procedure developed in this paper should also have applications in risk analysis, reliability, project planning, marketing, and other fields in which probability trees are used to model various processes of interest.