Archive for Non-Uniform Random Variate Generation
a minor typo, from 1986
Posted in Books, pictures, Statistics, University life with tags Bible, course, example, Luc Devroye, Monte Carlo course, multivariate Cauchy distribution, Non-Uniform Random Variate Generation, random number generation, reference, simulation, slides, typo on December 2, 2024 by xi'an
[very] simple rejection Monte Carlo
Posted in Books, pictures, R, University life with tags accept-reject algorithm, arXiv, Brussels, dominating measure, history of Monte Carlo, John von Neumann, Monte Carlo Statistical Methods, Non-Uniform Random Variate Generation, pseudo-random generator on March 29, 2024 by xi'an
“In recent years, the Rejection Monte Carlo (RMC) algorithm has emerged sporadically in literature under alternative names such as screening sampling or reject-accept sampling algorithms”
First, I was intrigued enough by a new arXival spotted in the Thalys train from Brussels to take a deeper look at it, but soon realised there was nothing of substance in the paper. Which solely recalls the fundamental of (accept-)reject algorithms, invented in the early days of computer simulation by von Neumann (even though the preprint refers to much more recent publications). Without providing the average acceptance probability as being equal to the inverse of the bounding constant [independently of the dimension of the random variable] and no mention of The Bible either… But with a standard depiction of accepted vs rejected points as uniformly dispersed on the subgraph of the proposal (as in the above taken from our very own Monte Carlo statistical Methods). Funnily enough, the most basic rejection algorithm, that is, the one based on a uniform sampling from a bounding (hyper)box is illustrated for a Normal target, although the latter has infinite support. And the paper seems to conclude on the appeal of using uniform proposals over bounding boxes, even though the increasing inefficiency against the dimension is well-known. A very simple rejection then, indeed!
why is this algorithm simulating a Normal variate?
Posted in Books, Kids, R, Statistics with tags cross validated, Devroye, Introducing Monte Carlo Methods with R, Luc Devroye, Monte Carlo Statistical Methods, Non-Uniform Random Variate Generation, normal generator, simulation on September 15, 2022 by xi'an
A backward question from X validated as to why the above is a valid Normal generator based on exponential generations. Which can be found in most textbooks (if not ours). And in The Bible, albeit as an exercise. The validation proceeds from the (standard) Exponential density dominating the (standard) Normal density and, according to Devroye, may have originated from von Neumann himself. But with a brilliant reverse engineering resolution by W. Huber on X validated. While a neat exercise, it requires on average 2.64 Uniform generations per Normal generation, against a 1/1 ratio for Box-Muller (1958) polar approach, or 1/0.86 for the Marsaglia-Bray (1964) composition-rejection method. The apex of the simulation jungle is however Marsaglia and Tsang (2000) ziggurat algorithm. At least on CPUs since, Note however that “The ziggurat algorithm gives a more efficient method for scalar processors (e.g. old CPUs), while the Box–Muller transform is superior for processors with vector units (e.g. GPUs or modern CPUs)” according to Wikipedia.
To draw a comparison between this Normal generator (that I will consider as von Neumann’s) and the Box-Müller polar generator,
#Box-Müller
bm=function(N){
a=sqrt(-2*log(runif(N/2)))
b=2*pi*runif(N/2)
return(c(a*sin(b),a*cos(b)))
}
#vonNeumann
vn=function(N){
u=-log(runif(2.64*N))
v=-2*log(runif(2.64*N))>(u-1)^2
w=(runif(2.64*N)<.5)-2
return((w*u)[v])
}
here are the relative computing times
> system.time(bm(1e8))
utilisateur système écoulé
7.015 0.649 7.674
> system.time(vn(1e8))
utilisateur système écoulé
42.483 5.713 48.222
simulating from the joint cdf
Posted in Books, Kids, pictures, R, Statistics, University life with tags Archimedean copulas, cross va, Gaussian copula, generator, inverse cdf, Monte Carlo, Non-Uniform Random Variate Generation, simulating copulas, simulation on July 13, 2022 by xi'an
An X validated question (what else?!) brought back (to me) the question of handling a bivariate cdf for simulation purposes. In the specific case of a copula when thus marginals were (well-)known…. And led me to an erroneous chain of thought, fortunately rescued by Robin Ryder! When the marginal distributions are set, the simulation setup is indeed equivalent to a joint Uniform simulation from a copula
![\mathbb P[U_1\leq u_1,U_2\leq u_2,\dots,U_d\leq u_d]=C(u_1,u_2,\dots,u_d)](https://s0.wp.com/latex.php?latex=%5Cmathbb+P%5BU_1%5Cleq+u_1%2CU_2%5Cleq+u_2%2C%5Cdots%2CU_d%5Cleq+u_d%5D%3DC%28u_1%2Cu_2%2C%5Cdots%2Cu_d%29&bg=000000&fg=B0B0B0&s=0&c=20201002)
In specific cases, as for instance the obvious example of Gaussian copulas, there exist customised simulation algorithms. Looking for more generic solutions, I turn to the Bible, where Chapter XI[an], has two entire sections XI.3.2. and XI.3.3 on the topic (even though Luc Devroye does not use the term copula there despite them being introduced in 1959 by A, Sklar, in response to a query of M. Fréchet). In addition to a study of copulas, both sections contain many specific solutions (as for instance in the [unnumbered] Table on page 585) but I found no generic simulation method. My [non-selected] answer to the question was thus to propose standard solutions such as finding one conditional since the marginals are Uniform. Which depends on the tractability of the derivatives of C(·,·).
However, being dissatisfied with this bland answer, I thought further about the problem and came up with a fallacious scheme, namely to first simulate the value p of C(U,V) by drawing a Uniform, and second simulate (U,V) conditional on C(U,V)=p. Going as far as running an R code on a simple copula, as shown above. Fallacious reasoning since (as I knew already!!!), C(U,V) is not uniformly distributed! But has instead a case-dependent distribution… As a (connected) aside, I wonder if the generator attached with Archimedean copulas has any magical feature that help with the generation of the associated copula.
A discrete Bernoulli factory
Posted in Books, Kids, Statistics with tags accept-reject algorithm, Bernoulli factory, cross validated, dice, Luc Devroye, Non-Uniform Random Variate Generation, riddle on October 18, 2021 by xi'an
A rather confusing (and now closed) question on X validated contained an interesting challenge of simulating an arbitrary discrete distribution using a single (standard) dice. It indeed made me think of the (more challenging) Bernoulli factory problem of simulating B(f(p)) using a B(p) simulator (with p unknown). I still do not see what the optimal solution is but the core challenge is to avoid simulating U(0,1) variate by exploiting the discrete nature of the target. Which may be an issue if the probabilities of the target are irrational and one is considering the cdf inversion approach. An alternative is to use an accept-reject approach, which also works for discrete distributions, by first deriving an instrumental distribution on the discrete support of the target from dice rolls, second finding the maximum of the ratio instrument to target, and third devising a discrete approach to selecting a generation with a probability taking a finite number of values. Which may prove quite costly. Finally, the least debatable approach is to turn the dice into a Uniform generator by using each draw as a digit in the base 5 representation of this Uniform variate, up to the precision desired for the resolution, and then apply the most efficient algorithm for the target distribution.
