Archive for Gaussian random walk

escaping the dark side of the Moon

Posted in Books, pictures, Statistics, University life with tags , , , , , , , , , , , on March 18, 2026 by xi'an

Sub-Cauchy Sampling: Escaping the Dark Side of the Moon was recently posted on arXiv by Sebastiano Grazzi (Warwick), Sifan Liu, Gareth O. Roberts (Warwick), and Jun Yang. With an hommage to Pink Floyd’s 1973 album both Gareth and I listened to at the time. (This was for sure my first Pink Floyd album!)

This highly original work is a sequel to the stereographic projection paper by Yang, Latuszýnski and Roberts (which was itself vaguely connected to our unpublished origami sampler). As in the stereographic projection method, the Euclidean space supporting the target is turned into a spherical cap of a hyper-sphere, referred to as the complement of the dark side of the Moon (or its bright side), and defined with respect to an observer ο who was at the north pole in the original method. The proposed MCMC algorithm, the Sub-Cauchy Projection Sampler (SCS), is a random-walk-type Metropolis algorithm on the bright side and it gets its name from being uniformly ergodic for sub-Cauchy targets. An explanation for this massive achievement is that points at infinity in the Euclidean space are now mapped to the (d − 1)-dimensional boundary of the dark side rather than at the north pole of the hypersphere. Meaning that the push-forward density may remain bounded. (The random walk on the bright side involves projections for proposals ending on the dark side, while keeping the target intact.) There are several calibration parameters to the algorithm that can be tuned by variational arguments (and the goal of getting near a uniform distribution over the bright side), since optimal acceptance rates no longer apply.

Metropolis in 95 characters

Posted in pictures, R, Statistics, Travel with tags , , , , , , , , on January 2, 2020 by xi'an

Here is an R function that produces a Metropolis-Hastings sample for the univariate log-target f when the later is defined outside as another function. And when using a Gaussian random walk with scale one as proposal. (Inspired from a X validated question.)

m<-function(T,y=rnorm(1))ifelse(rep(T>1,T),
  c(y*{f({z<-m(T-1)}[1])-f(y+z[1])<rexp(1)}+z[1],z),y)

The function is definitely not optimal, crashes for values of T larger than 580 (unless one modifies the stack size), and operates the most basic version of a Metropolis-Hastings algorithm. But as a codegolf challenge (on a late plane ride), this was a fun exercise.

random wake

Posted in Books, Kids, R, Statistics with tags , , , , , , on December 27, 2017 by xi'an

Just too often on X validated, one sees questions displaying a complete ignorance of the basics that makes one purposelessly wonder what is the point of trying to implement advanced methods when missing the necessary background. And just as often, I reacted to the question by wondering out loud about this… In the current case, the question was about debugging an R code for a mixture of two exponential distributions and the random walk Metropolis algorithm that comes with it. Except that the Normal noise was replaced with a Uniform U(0,1) noise, leading to a most obviously transient Markov chain.I eventually corrected the R code, which returned a rather nicely label-switching output. And which did not necessarily help with the comprehension of the fundamentals.

normal variates in Metropolis step

Posted in Books, Kids, R, Statistics, University life with tags , , , , , , , , on November 14, 2017 by xi'an

A definitely puzzled participant on X validated, confusing the Normal variate or variable used in the random walk Metropolis-Hastings step with its Normal density… It took some cumulated efforts to point out the distinction. Especially as the originator of the question had a rather strong a priori about his or her background:

“I take issue with your assumption that advice on the Metropolis Algorithm is useless to me because of my ignorance of variates. I am currently taking an experimental course on Bayesian data inference and I’m enjoying it very much, i believe i have a relatively good understanding of the algorithm, but i was unclear about this specific.”

despite pondering the meaning of the call to rnorm(1)… I will keep this question in store to use in class when I teach Metropolis-Hastings in a couple of weeks.

puzzled by harmony [not!]

Posted in Books, Kids, Mountains, pictures, R, Running, Statistics, Travel with tags , , , , , on December 13, 2016 by xi'an

In answering yet another question on X validated about the numerical approximation of the marginal likelihood, I suggested using an harmonic mean estimate as a simple but worthless solution based on an MCMC posterior sample. This was on a toy example with a uniform prior on (0,π) and a “likelihood” equal to sin(θ) [really a toy problem!]. Simulating an MCMC chain by a random walk Metropolis-Hastings algorithm is straightforward, as is returning the harmonic mean of the sin(θ)’s.

f <- function(x){
    if ((0<x)&(x<pi)){
        return(sin(x))}else{
        return(0)}}

n = 2000 #number of iterations
sigma = 0.5
x = runif(1,0,pi) #initial x value
chain = fx = f(x)   
#generates an array of random x values from norm distribution
rands = rnorm(n,0, sigma) 
#Metropolis - Hastings algorithm
for (i in 2:n){
    can = x + rands[i]  #candidate for jump
    fcan=f(can)
    aprob = fcan/fx #acceptance probability
    if (runif(1) < aprob){
        x = can
        fx = fcan}
    chain=c(chain,fx)}
I = pi*length(chain)/sum(1/chain) #integral harmonic approximation

However, the outcome looks remarkably stable and close to the expected value 2/π, despite 1/sin(θ) having an infinite integral on (0,π). Meaning that the average of the 1/sin(θ)’s has no variance. Hence I wonder why this specific example does not lead to an unreliable output… But re-running the chain with a smaller scale σ starts producing values of sin(θ) regularly closer to zero, which leads to an estimate of I both farther away from 2 and much more variable. No miracle, in the end!