Variance of a standard vector Monte Carlo estimate in the theory of polarized radiative transfer

The spectral radius ρ of the matrix integral operator defining the covariance matrix of a standard vector Monte Carlo estimate in the polarized radiative transfer problem is examined. The theory of positive operators is used to analytically calculate ρ = ρ₀ for transfer through an infinite homogeneous medium. For a bounded medium, it is shown that ρ is approximately equal to ρ₀ times the spectral radius of the operator corresponding to radiative transfer without polarization. This is shown numerically by estimating the iterations of the corresponding resolvent and approximately analytically by using a perturbation of a special functional.     

Monte Carlo study of time asymptotics of the polarized radiation intensity

The parameters of time asymptotics of the polarized radiation intensity are estimated. Precision Monte Carlo estimates of these parameters are derived for finite medium layers by iterating the resolvent of the corresponding transfer operator with a given scattering matrix and by evaluating parametric time derivatives. The computations are performed for two versions of the problem: with a Rayleigh scattering matrix and an aerosol scattering matrix. It is shown that the asymptotics of the radiation intensity are affected by polarization, except for the spatially homogeneous problem, for which the results are obtained analytically.     

Modifications of weighted Monte Carlo algorithms for nonlinear kinetic equations

Test problems for the nonlinear Boltzmann and Smoluchowski kinetic equations are used to analyze the efficiency of various versions of weighted importance modeling as applied to the evolution of multiparticle ensembles. For coagulation problems, a considerable gain in computational costs is achieved via the approximate importance modeling of the “free path” of the ensemble combined with the importance modeling of the index of a pair of interacting particles. A weighted modification of the modeling of the initial velocity distribution was found to be the most efficient for model solutions to the Boltzmann equation. The technique developed can be useful as applied to real-life coagulation and relaxation problems for which the model problems considered give approximate solutions.     

The Monte Carlo methods for solving the vector and Stochastic Helmholtz equations

Dedicated to the jubilee of the 70th birthday of my teacher Academician G. I. Marchuk.     

Dual representation of the mean square of the Monte Carlo vector estimator

Doklady Mathematics -     

An analysis of a Monte Carlo algorithm for estimating the permanent

Karmarkar, Karp, Lipton, Lovász, and Luby proposed a Monte Carlo algorithm for approximating the permanent of a non-negativen×n matrix, which is based on an easily computed, unbiased estimator. It is not difficult to construct 0,1-matrices for which the variance of this estimator is very large, so that an exponential number of trials is necessary to obtain a reliable approximation that is within a constant factor of the correct value. Nevertheless, the same authors conjectured that for a random 0,1-matrix the variance of the estimator is typically small. The conjecture is shown to be true; indeed, for almost every 0,1-matrixA, just O(nw(n)e ^-2) trials suffice to obtain a reliable approximation to the permanent ofA within a factor 1±ɛ of the correct value. Here ω(n) is any function tending to infinity asn→∞. This result extends to random 0,1-matrices with density at leastn ^−1/2ω(n). It is also shown that polynomially many trials suffice to approximate the permanent of any dense 0,1-matrix, i.e., one in which every row- and column-sum is at least (1/2+α)n, for some constant α>0. The degree of the polynomial bounding the number of trials is a function of α, and increases as α→0. Supported by NSF grant CCR-9225008. The work described here was partly carried out while the author was visiting Princeton University as a guest of DIMACS (Center for Discrete Mathematics and Computer Science).     

Monte Carlo analysis of convertible bonds with reset clauses

This paper analyzes some features of non-callable convertible bonds with reset clauses via both analytic and Monte Carlo simulation approaches. Assume that the underlying stock receives no dividends and that it has credit risk of the issuer. We mean by reset that the conversion price is adjusted downwards if the underlying stock price does not exceed pre-specified prices. Reset convertibles are usually issued when the outlook for the issuer is unfavorable. The price of any convertible bonds can be approximately viewed as a sum of values of an otherwise identical non-convertible bond plus an embedded option to convert the bond into the underlying stock. In this paper, we first develop an exact formula for the conversion option value of the European riskless convertible in the classical Black–Scholes–Merton framework. It is shown by Monte Carlo simulation that conversion option value estimates of the American risky convertible are located in a certain region defined by this formula. From estimates of the conversion probability, it is also shown that there exists an optimal reset time in the latter half of the trading interval.     

Further analysis of multilevel Monte Carlo methods for elliptic PDEs with random coefficients

We consider the application of multilevel Monte Carlo methods to elliptic PDEs with random coefficients. We focus on models of the random coefficient that lack uniform ellipticity and boundedness with respect to the random parameter, and that only have limited spatial regularity. We extend the finite element error analysis for this type of equation, carried out in Charrier et al. (SIAM J Numer Anal, 2013), to more difficult problems, posed on non-smooth domains and with discontinuities in the coefficient. For this wider class of model problem, we prove convergence of the multilevel Monte Carlo algorithm for estimating any bounded, linear functional and any continuously Fréchet differentiable non-linear functional of the solution. We further improve the performance of the multilevel estimator by introducing level dependent truncations of the Karhunen–Loève expansion of the random coefficient. Numerical results complete the paper.     

Implementations of the Monte Carlo EM Algorithm

The Monte Carlo EM (MCEM) algorithm is a modification of the EM algorithm where the expectation in the E-step is computed numerically through Monte Carlo simulations. The most exible and generally applicable approach to obtaining a Monte Carlo sample in each iteration of an MCEM algorithm is through Markov chain Monte Carlo (MCMC) routines such as the Gibbs and Metropolis–Hastings samplers. Although MCMC estimation presents a tractable solution to problems where the E-step is not available in closed form, two issues arise when implementing this MCEM routine: (1) how do we minimize the computational cost in obtaining an MCMC sample? and (2) how do we choose the Monte Carlo sample size? We address the first question through an application of importance sampling whereby samples drawn during previous EM iterations are recycled rather than running an MCMC sampler each MCEM iteration. The second question is addressed through an application of regenerative simulation. We obtain approximate independent and identical samples by subsampling the generated MCMC sample during different renewal periods. Standard central limit theorems may thus be used to gauge Monte Carlo error. In particular, we apply an automated rule for increasing the Monte Carlo sample size when the Monte Carlo error overwhelms the EM estimate at any given iteration. We illustrate our MCEM algorithm through analyses of two datasets fit by generalized linear mixed models. As a part of these applications, we demonstrate the improvement in computational cost and efficiency of our routine over alternative MCEM strategies.     

M5 model tree and Monte Carlo simulation for efficient structural reliability analysis

Practically, the performance of many engineering problems can be defined using a complex implicit limit state function. Approximation of the accurate failure probability is very time-consuming and inefficient based on Monte Carlo simulation (MCS) for complex performance functions. M5 model tree (M5Tree) model is robust approach for simulation and prediction phenomena, which provides ability to dealing with complex implicit problems by dividing them into smaller problems. By improving the efficiency of reliability method using accurate approximated failure probability, an efficient reliability method using the MCS and M5Tree is proposed to calibrate the performance function and estimate the failure probability, respectively. The superiorities including simplicity and accuracy of M5Tree meta-model are investigated to evaluate the actual performance function through five nonlinear complex mathematical and structural reliability problems. The proposed reliability method-based MCS and M5Tree improved the computational efforts for evaluating the performance function in reliability analysis. The M5Tree significantly increased the efficiency of reliability analysis with accurate failure probability.     

A Monte Carlo Estimation of the Entropy for Markov Chains

We introduce an estimate of the entropy of the marginal density p ^t of a (eventually inhomogeneous) Markov chain at time t≥1. This estimate is based on a double Monte Carlo integration over simulated i.i.d. copies of the Markov chain, whose transition density kernel is supposed to be known. The technique is extended to compute the external entropy , where the p ₁ ^t s are the successive marginal densities of another Markov process at time t. We prove, under mild conditions, weak consistency and asymptotic normality of both estimators. The strong consistency is also obtained under stronger assumptions. These estimators can be used to study by simulation the convergence of p ^t to its stationary distribution. Potential applications for this work are presented: (1) a diagnostic by simulation of the stability property of a Markovian dynamical system with respect to various initial conditions; (2) a study of the rate in the Central Limit Theorem for i.i.d. random variables. Simulated examples are provided as illustration.      

Monte Carlo analysis of estimation methods for the prediction of customer response patterns in direct marketing

In direct marketing, customers are usually asked to take a specific action, and their responses are recorded over time and stored in a database. Based on the response data, we can estimate the number of customers who will ultimately respond, the number of responses anticipated to receive by a certain period of time, and the like. The goal of this article is to derive and propose several estimation methods and compare their performances in a Monte Carlo simulation. The response patterns can be described by a simple geometric function, which relates the number of responses to elapsed time. The “maximum likelihood” estimator appears to be the most effective method of estimating the parameters of this function. As we have more sample observations, the maximum likelihood estimates also converge to the true parameter values rapidly.     

Weight Monte Carlo algorithms for estimation and parametric analysis of the solution to the kinetic coagulation equation

The Smoluchowski equation with linear coagulation coefficients depending on two parameters is considered. We construct weight algorithms for estimating various linear functionals in an ensemble that is governed by the equation under study. The algorithms constructed allow us to estimate the functionals for various parameters, as well as parametric derivatives by using the same set of trajectories. Moreover, we construct the value algorithms and analyze their efficiency for estimating the total monomer concentration, as well as the total monomer and dimer concentration in the ensemble. The computational cost is considerably reduced via the approximate value simulation of the time between interactions combined with the value simulation of the interacting pair number.     

Smoothed Monte Carlo estimators for the time-in-the-red in risk processes

We consider a modified version of the de Finetti model in insurance risk theory in which, when surpluses become negative the company has the possibility of borrowing, and thus continue its operation. For this model we examine the problem of estimating the time-in-the red over a finite horizon via simulation. We propose a smoothed estimator based on a conditioning argument which is very simple to implement as well as particularly efficient, especially when the claim distribution is heavy tailed. We establish unbiasedness for this estimator and show that its variance is lower than the naïve estimator based on counts. Finally we present a number of simulation results showing that the smoothed estimator has variance which is often significantly lower than that of the naïve Monte-Carlo estimator.     

A Monte Carlo method for the Broadwell model with relaxation

A Monte Carlo method is proposed for numerical solution of the Broadwell model. Developing a probabilistic interpretation of the equations, the transport and collision parts are treated separately in the method. Particles are advected according to their velocities and collisions are performed between randomly chosen particles. We numerically test the algorithm for a variety of examples. In particular we are interested in situations which generate structures that have nonsmooth fronts. Our simulations show that this Monte Carlo method is capable of capturing the nonlinear regime in presence of shocks and interactions. (© 2005 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Parallel Speed-Up of Monte Carlo Methods for Global Optimization

We introduce the notion of expected hitting time to a goal as a measure of the convergence rate of a Monte Carlo optimization method. The techniques developed apply to simulated annealing, genetic algorithms, and other stochastic search schemes. The expected hitting time can itself be calculated from the more fundamental complementary hitting time distribution (CHTD) which completely characterizes a Monte Carlo method. The CHTD is asymptotically a geometric series, (1/s)/(1 − λ), characterized by two parameters, s, λ, related to the search process in a simple way. The main utility of the CHTD is in comparing Monte Carlo algorithms. In particular we show that independent, identical Monte Carlo algorithms run in parallel, IIP parallelism, and exhibit superlinear speedup. We give conditions under which this occurs and note that equally likely search is linearly sped up. Further we observe that a serial Monte Carlo search can have an infinite expected hitting time, but the same algorithm when parallelized can have a finite expected hitting time. One consequence of the observed superlinear speedup is an improved uniprocessor algorithm by the technique of in-code parallelism.     

A Monte Carlo simulation approach to the characterization of uncertainties in cancer staging and radiation treatment decisions

Radiation treatment (RT) for cancer is a critical medical procedure that occurs in a complex environment that is subject to uncertainties and errors. We employed a simulation (a variant of Monte Carlo) model that followed a cohort of hypothetical breast cancer patients to estimate the probability of incorrect staging and treatment decisions. As inputs, we used a combination of literature information and expert judgement. Input variables were defined as probability distributions within the model. Uncertainties were propagated via simulation. Sensitivity and value-of-information analyses were then conducted to quantify the effect of variable uncertainty on the model outputs. We found a small but non-trivial probability that patients would be incorrectly staged and thus be subjected to inappropriate treatment. Some routinely used tests for staging and metastasis detection have very limited informational value. This work has implications for the methods used in cancer staging and subsequent risk assessment of treatment errors.