## ratio-of-uniforms [#4]

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

Possibly the last post on random number generation by Kinderman and Monahan’s (1977) ratio-of-uniform method. After fiddling with the Gamma(a,1) distribution when a<1 for a while, I indeed figured out a way to produce a bounded set with this method: considering an arbitrary cdf Φ with corresponding pdf φ, the uniform distribution on the set Λ of (u,v)’s in R⁺xX such that

0≤u≤Φοƒ[φοΦ⁻¹(u)v]

induces the distribution with density proportional to ƒ on φοΦ⁻¹(U)V. This set Λ has a boundary that is parameterised as

u=Φοƒ(x),  v=1/φοƒ(x), x∈Χ

which remains bounded in u since Φ is a cdf and in v if φ has fat enough tails. At both 0 and ∞. When ƒ is the Gamma(a,1) density this can be achieved if φ behaves like log(x)² near zero and like a inverse power at infinity. Without getting into all the gory details, closed form density φ and cdf Φ can be constructed for all a’s, as shown for a=½ by the boundaries in u and v (yellow) below

which leads to a bounded associated set Λ

At this stage, I remain uncertain of the relevance of such derivations, if only because the set A thus derived is ill-suited for uniform draws proposed on the enclosing square box. And also because a Gamma(a,1) simulation can rather simply be derived from a Gamma(a+1,1) simulation. But, who knows?!, there may be alternative usages of this representation, such as innovative slice samplers. Which means the ratio-of-uniform method may reappear on the ‘Og one of those days…

## seven permanent positions at Warwick!

Posted in pictures, Statistics, University life with tags , , , , , on November 22, 2016 by xi'an

Seven academic positions in Statistics are currently opening at the University of Warwick! They correspond to the following levels, all permanent but for the Harrison Early-Career Assistant Professorship:

• Harrison Early-Career Assistant Professor of Statistics (3 years)
• Assistant Professor of Statistics
• Senior Teaching Fellow in Statistics
• Assistant or Associate Professor of Financial Mathematics
• Assistant or Associate Professor of Statistics (two positions)
• Full Professor of Statistics

Applicants are sought with expertise in Statistics (in the wide sense, including both applied and methodological statistics, probability, probabilistic operational research and mathematical finance together with interdisciplinary topics involving one or more of these areas). Applicants for senior positions should have an excellent publication record and ability to secure research funding. Applicants for more junior positions should show exceptional promise to become leading academics. [More details.]

Informal enquires should be addressed to any of Professors Mark Steel, David Hobson, Gareth Roberts, Wilfrid Kendall, or to any other senior member of the Department. Including me. The closing date is 3 January 2017 for all positions, except for the Full Professor position which is set to 10 January 2017.

As a completely objective observer (!), I can state that this is a fantastic department, with a strong sense of community and support among the faculty, and a tremendously diverse array of research activities and topics. And thus encourage anyone interested in joining to apply for some of those positions!

## pseudo-marginal MCMC with minimal replicas

Posted in Books, Statistics, University life with tags , , , , , , on November 16, 2016 by xi'an

A week ago, Chris Sherlock, Alexandre Thiery and Anthony Lee (Warwick) arXived a short paper where they show that it is most efficient to use only one random substitute to the likelihood when considering a pseudo-marginal algorithm. Thus generalising an earlier paper of Luke Bornn and co-authors I commented earlier. A neat side result in the paper is that the acceptance probability for m replicas [in the pseudo-marginal approximation] is at most m/s the acceptance probability for s replicas when s<m. The main result is as follows:

There is a (superficial?) connection with the fact that when running Metropolis-within-Gibbs there is no advantage in doing more than one single Metropolis iteration, as the goal is to converge to the joint posterior, rather than approximating better the full conditional…

## snapshot from Warwick [jatp]

Posted in pictures, Travel, University life with tags , , , , , on November 10, 2016 by xi'an

## A tyring day…

Posted in pictures, Running, Travel with tags , , , , on November 9, 2016 by xi'an

An exciting and explosive Election day in Warwick, if unconnected with the US elections, since my bike front tube first lost its valve to the pump [if not to the trump!] and then the replacement tube exploded one hour later, presumably a combination of high pressure and hot temperatures. [Disclaimer: Any similarity with current events and overblown egos is purely based on hot air. Any hope of seeing the situation deflating nicely is alas fizzling out really fast…]

## ratio-of-uniforms [#3]

Posted in Books, pictures, R, Statistics with tags , , , , , on November 4, 2016 by xi'an

Being still puzzled (!) by the ratio-of-uniform approach, mostly failing to catch its relevance for either standard distributions in a era when computing a cosine or an exponential is negligible, or non-standard distributions for which computing bounds and boundaries is out-of-reach, I kept searching for solutions that would include unbounded densities and still produce compact boxes, as this seems essential for accept-reject simulation if not for slice sampling. And after exploring some dead-ends (in tune with running in Venezia!), I came upon the case of the generalised logistic transform

$h(\omega)=\omega^a/(1+\omega^a)$

which ensures that the [ratio-of-almost-uniform] set I defined in my slides last week

$\mathfrak{H}=\left\{(u,v);\ 0\le u\le h(f(v/g(u))\right\}$

is bounded in u. Since the transform g is the derivative of the inverse of h (!),

$g(y)=a^{-1}y^{(1-a)/a}/(1-y)^{(1-3a)/a}$

the parametrisation of the boundary of H is

$u(x)=f(x)^a/(1+f(x)^a)\ v(x)=a^{-1}xf(x)^{(a-1)/a}(1+f(x)^a)^2$

which means it remains bounded if (a) a≤1 [to ensure boundedness at infinity] and (b) the limit of v(x) at zero [where I assume the asymptote stands] is bounded. Meaning

$\lim_{x\to 0} xf(x)^{2a+1/a-1}<\infty$

Working a wee bit more on the problem led me to realise that resorting to an arbitrary cdf Φ instead of the logistic cdf could solve the problem for most distributions, including all Gammas! Indeed, the boundary of H is now

$u(x)=\Phi(f(x))^a\ v(x)=a^{-1}xf(x)^{(a-1)/a}/\varphi(f(x))$

which means it remains bounded if φ has very heavy tails, like 1/x². To handle the explosion when x=0. And an asymptote itself at zero, to handle the limit at infinity when f(x) goes to zero.

## ratio-of-uniforms

Posted in Books, pictures, R, Statistics with tags , , , , on October 24, 2016 by xi'an

One approach to random number generation that had always intrigued me is Kinderman and Monahan’s (1977) ratio-of-uniform method. The method is based on the result that the uniform distribution on the set A of (u,v)’s in R⁺xX such that

0≤u²≤ƒ(v/u)

induces the distribution with density proportional to ƒ on V/U. Hence the name. The proof is straightforward and the result can be seen as a consequence of the fundamental lemma of simulation, namely that simulating from the uniform distribution on the set B of (w,x)’s in R⁺xX such that

0≤w≤ƒ(x)

induces the marginal distribution with density proportional to ƒ on X. There is no mathematical issue with this result, but I have difficulties with picturing the construction of efficient random number generators based on this principle.

I thus took the opportunity of the second season of [the Warwick reading group on] Non-uniform random variate generation to look anew at this approach. (Note that the book is freely available on Luc Devroye’s website.) The first thing I considered is the shape of the set A. Which has nothing intuitive about it! Luc then mentions (p.195) that the boundary of A is given by

u(x)=√ƒ(x),v(x)=x√ƒ(x)

which then leads to bounding both ƒ and x→x²ƒ(x) to create a box around A and an accept-reject strategy, but I have trouble with this result without making further assumptions about ƒ… Using a two component normal mixture as a benchmark, I found bounds on u(.) and v(.) and simulated a large number of points within the box to end up with the above graph that indeed the accepted (u,v)’s were within this boundary. And the same holds with a more ambitious mixture: