Florian Puchhammer
@bcamath.org
Modelling and Simulation in Life and Material Sciences
Basque Center for Applied Mathematics
RESEARCH INTERESTS
- Monte Carlo and quasi-Monte Carlo Methods
- Numerical mathematics
- Stochastic simulation
- Statistics
- Markov chain Monte Carlo and Hamiltonian Monte Carlo
- Probability
- Operations research
- Mathematical programming
Scopus Publications
Scopus Publications
Pierre L’Ecuyer, Florian Puchhammer, and Amal Ben Abdellah
Institute for Operations Research and the Management Sciences (INFORMS)
Estimating the unknown density from which a given independent sample originates is more difficult than estimating the mean in the sense that, for the best popular nonparametric density estimators, the mean integrated square error converges more slowly than at the canonical rate of [Formula: see text]. When the sample is generated from a simulation model and we have control over how this is done, we can do better. We examine an approach in which conditional Monte Carlo yields, under certain conditions, a random conditional density that is an unbiased estimator of the true density at any point. By averaging independent replications, we obtain a density estimator that converges at a faster rate than the usual ones. Moreover, combining this new type of estimator with randomized quasi–Monte Carlo to generate the samples typically brings a larger improvement on the error and convergence rate than for the usual estimators because the new estimator is smoother as a function of the underlying uniform random numbers. Summary of Contribution: Stochastic simulation is commonly used to estimate the mathematical expectation of some output random variable X together with a confidence interval for this expectation. But the simulations usually provide information to do much more, such as estimating the entire distribution (or density) of X. Histograms are routinely provided by standard simulation software, but they are very primitive density estimators. Kernel density estimators perform better, but they are trickier to use, have bias, and their mean square error converges more slowly than the canonical rate of O(1/n) with n independent samples. In this paper, we explain how to construct unbiased density estimators that converge at the canonical rate and even much faster when combined with randomized quasi–Monte Carlo. The key idea is to use conditional Monte Carlo to hide appropriate information and obtain a computable (random) conditional density, which acts (under certain conditions) as an unbiased density estimator. Moreover, this sample density is typically smoother than the classic density estimators as a function of the underlying uniform random numbers, so it can get along much better with randomized quasi–Monte Carlo methods. This offers an opportunity to further improve the O(1/n) rate. We observe rates near O(1/n2) on some examples, and we give conditions under which this type of rate provably holds. The proposed approach is simple, easy to implement, and extremely effective, so it provides a significant addition to the stochastic simulation toolbox.
Florian Puchhammer and Pierre L'Ecuyer
IEEE
We consider the problem of estimating the density of a random variable X which is the output of a simulation model. We show how an unbiased density estimator can be constructed via the classical likelihood ratio derivative estimation method proposed over 35 years ago by Glynn, Rubinstein, and others. We then extend this density estimation method to cover situations where it does not apply directly. What we obtain is closely related to the generalized likelihood ratio method proposed recently by Peng and his co-authors, although the assumptions differ. We compare the methods and assumptions on some examples.
Pierre L’Ecuyer and Florian Puchhammer
Springer International Publishing
Pierre L’Ecuyer, Pierre Marion, Maxime Godin, and Florian Puchhammer
Springer International Publishing
Florian Puchhammer, Amal Ben Abdellah, and Pierre L’Ecuyer
Springer Science and Business Media LLC
We explore the use of Array-RQMC, a randomized quasi-Monte Carlo method designed for the simulation of Markov chains, to reduce the variance when simulating stochastic biological or chemical reaction networks with [Formula: see text]-leaping. The task is to estimate the expectation of a function of molecule copy numbers at a given future time T by the sample average over n sample paths, and the goal is to reduce the variance of this sample-average estimator. We find that when the method is properly applied, variance reductions by factors in the thousands can be obtained. These factors are much larger than those observed previously by other authors who tried RQMC methods for the same examples. Array-RQMC simulates an array of realizations of the Markov chain and requires a sorting function to reorder these chains according to their states, after each step. The choice of sorting function is a key ingredient for the efficiency of the method, although in our experiments, Array-RQMC was never worse than ordinary Monte Carlo, regardless of the sorting method. The expected number of reactions of each type per step also has an impact on the efficiency gain.
Amal Ben Abdellah, Pierre L'Ecuyer, Art B. Owen, and Florian Puchhammer
Society for Industrial & Applied Mathematics (SIAM)
We consider the problem of estimating the density of a random variable $X$ that can be sampled exactly by Monte Carlo (MC). We investigate the effectiveness of replacing MC by randomized quasi Monte Carlo (RQMC) or by stratified sampling over the unit cube, to reduce the integrated variance (IV) and the mean integrated square error (MISE) for kernel density estimators. We show theoretically and empirically that the RQMC and stratified estimators can achieve substantial reductions of the IV and the MISE, and even faster convergence rates than MC in some situations, while leaving the bias unchanged. We also show that the variance bounds obtained via a traditional Koksma-Hlawka-type inequality for RQMC are much too loose to be useful when the dimension of the problem exceeds a few units. We describe an alternative way to estimate the IV, a good bandwidth, and the MISE, under RQMC or stratification, and we show empirically that in some situations, the MISE can be reduced significantly even in high-dimensional settings.
Amal Ben Abdellah, Pierre L'Ecuyer, and Florian Puchhammer
IEEE
Array-RQMC has been proposed as a way to effectively apply randomized quasi-Monte Carlo (RQMC) when simulating a Markov chain over a large number of steps to estimate an expected cost or reward. The method can be very effective when the state of the chain has low dimension. For pricing an Asian option under an ordinary geometric Brownian motion model, for example, Array-RQMC can reduce the variance by factors in the millions. In this paper, we show how to apply this method and we study its effectiveness in case the underlying process has stochastic volatility. We show that Array-RQMC can also work very well for these models, even if it requires RQMC points in larger dimension. We examine in particular the variance-gamma, Heston, and Ornstein-Uhlenbeck stochastic volatility models, and we provide numerical results.
Florian Puchhammer
Elsevier BV
Following a result of D.~Bylik and M.T.~Lacey from 2008 it is known that there exists an absolute constant $\\eta>0$ such that the (unnormalized) $L^{\\infty}$-norm of the three-dimensional discrepancy function, i.e, the (unnormalized) star discrepancy $D^{\\ast}_N$, is bounded from below by $D_{N}^{\\ast}\\geq c (\\log N)^{1+\\eta}$, for all $N\\in\\mathbb{N}$ sufficiently large, where $c>0$ is some constant independent of $N$. This paper builds upon their methods to verify that the above result holds with $\\eta<1/(32+4\\sqrt{41})\\approx 0.017357\\ldots$
Roswitha Hofer and Florian Puchhammer
Institute of Mathematics, Polish Academy of Sciences
B. Eichinger, F. Puchhammer, and P. Yuditskii
Springer Science and Business Media LLC
We give a free parametric representation for the coefficient sequences of Jacobi matrices whose spectral measures satisfy the Killip–Simon condition with respect to two (arbitrary) disjoint intervals. This parametrization is given by means of the Jacobi flow on SMP matrices, which we introduce here.