Archive for the R Category

[very] simple rejection Monte Carlo

Posted in Books, pictures, R, University life with tags , , , , , , , , 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!

insufficient Gibbs sampling bridges as well!

Posted in Books, Kids, pictures, R, Statistics, University life with tags , , , , , , , , , , , , on March 3, 2024 by xi'an

Antoine Luciano, Robin Ryder and I posted a revised version of our insufficient Gibbs sampler on arXiv last week (along with three other revisions or new deposits of mine’s!), following comments and suggestions from referees. Thanks to this revision, we realised that the evidence based on an (insufficient) statistic was also available for approximation by a Monte Carlo estimate attached to the completed sample simulated by the insufficient sampler. Better, a bridge sampling estimator can be used in the same conditions as when the full data is available! In this new version, we thus revisited toy examples first explored in some of my ABC papers on testing (with insufficient statistics), as illustrated by both graphs on this post.

merely fiddlin

Posted in Books, Kids, Mountains, pictures, R, Running, Travel with tags , , , , , , , , , on February 25, 2024 by xi'an

For any positive, base-10 integer N, define f(N) as the number of times you have to add up its digits until you get a one-digit number. Find the smallest whole number N such that f(N) = 4. Fiddler on the Proof

Although I first gave in to the R brute force attempt, it failed to return a value for N up to 10⁹, but the solution is obvious as f(N)=p iff F(s(N))=p-1. Since f(19)=2 and f(199)=3, a number N whose sum of digits is 199 is legit, for instance

f(99999999999999999999991)=4

where there are 22 9s. And although it gets a wee bit boring, here is the associated attempt by ChatGPT³

simulating signed mixtures

Posted in Books, pictures, R, Statistics, University life with tags , , , , , , , , on February 2, 2024 by xi'an

While simulating from a mixture of standard densities is relatively straightforward, when the component densities are easily simulated, to the point that many simulation methods exploit an intermediary mixture construction to speed up the production of pseudo-random samples from more challenging distributions (see Devroye, 1986), things get surprisingly more complicated when the mixture weights can take negative values. For instance, the naïve solution consisting in first simulating from the associated mixture of positive weight components and then using an accept-reject step may prove highly inefficient since the overall probability of acceptance

{\displaystyle 1}\Big/{\displaystyle \sum_{k=1}^{P} \omega_k^+}

is the inverse of the sum of the positive weights and hence can be arbitrarily close to zero. The intuition for such inefficiency is that simulating from the positive weight components need not produce values within regions of high probability for the actual distribution

m = \sum_{k=1}^P \omega_k^+ f_k - \sum_{k=1}^N \omega_k^- g_k

since its negative weight components may remove most of the mass under the positive weight components. In other words, the negative weight components do not have a natural latent variable interpretation and the resulting mixture can be anything, as the above graph testifies.

Julien Stoehr (Paris Dauphine) and I started investigating this interesting challenge when the Master students who had been exposed to said challenge could not dent it in any meaningful way. We have now arXived a specific algorithm that proves superior to the naïve accept-reject algorithm, but also to the numerical cdf inversion (which happens to be available in this setting). Compared with the naïve version, we construct an alternative accept-reject scheme based on pairing positive and negative components as well as possible, partitioning the real line, and finding tighter upper and lower bounds on positive and negative components, respectively, towards yielding a higher acceptance rate on average. Designing a random generator of signed mixtures with enough variability and representativity proved a challenge in itself!

simulating Gumbel’s bivariate exponential distribution

Posted in Books, Kids, R, Statistics with tags , , , , , , , , , , on January 14, 2024 by xi'an

A challenge interesting enough for a sunny New Year morn, found on X validated, namely the simulation of a bivariate exponential distribution proposed by Gumbel in 1960, with density over the positive quadrant in IR²

{}^{ [(\lambda_2+rx_1)(\lambda_1+rx_2)-r]\exp[-(\lambda_1x_1+\lambda_2x_2+rx_1x_2)]}

Although there exists a direct approach based on the fact that the marginals are Exponential distributions and the conditionals signed mixtures of Gamma distributions, an accept-reject algorithm is also available for the pair, with a dominating density representing a genuine mixture of four Gammas, when omitting the X product in the exponential and the negative r in the first term. The efficiency of this accept-reject algorithm is high for r small. However, and in more direct connection with the original question, using this approach to integrate the function equal to the product of the pair, as considered in the original paper of Gumbel, is much less efficient than seeking a quasi-optimal importance function, since this importance function is yet another mixture of four Gammas that produces a much reduced variance at a cheaper cost!