4
$\begingroup$

I am looking for an uniform random walk algorithm on a simplex for MCMC purposes. Hence the process should on average spent the same time in any given area. I want this to be my proposal algorithm. The Dirichlet distribution has all sorts of problems at proposing samples near the boundaries. This is well documented.

The only algorithm I could find was https://en.wikipedia.org/wiki/User%3aSkinnerd/Simplex_Point_Picking. It looks like the algorithm will "struggle" to escape when near the boundary of the simplex i.e for values close to zero xnew = xold*random_number. Am I correct in that this algorithm is not a random walk on a simplex? or am I missing something? I implemented the algorithm and as I thought too many samples are at the boundary of the simplex.

The MCMC samples from the linked algorithm.

$\endgroup$
6
  • 2
    $\begingroup$ I wouldn't expect a random walk on a simplex to be uniform; what's the basis on which you say that there's too many near the boundary? How do you know what the behavior of a random walk on the simplex should look like near the boundary? $\endgroup$
    – Glen_b
    Sep 11 '17 at 11:10
  • $\begingroup$ Your random walk must be uniform, right? If so, I think you should add this point in your question. $\endgroup$ Sep 11 '17 at 11:10
  • $\begingroup$ The link describes how to perform a uniform random walk. I should have made it clearer that I wanted a uniform random walk. Edited. $\endgroup$
    – Travis
    Sep 11 '17 at 11:40
  • 1
    $\begingroup$ ^ the algorithm described in the linked Wikipedia user's page is erroneous, as explained in my answer to the aforementioned question $\endgroup$ Sep 11 '17 at 15:52
  • $\begingroup$ @JuhoKokkala I did try your suggestion using the correct Metropolis Hastings algorithm, but then the samples are even less uniform. $\endgroup$
    – Travis
    Sep 12 '17 at 6:36
6
$\begingroup$

I understand your simplex to be the specific subset

$$\Delta^d = \{(x_1,x_2,\ldots, x_d)\mid 0 \le x_1 \le x_2 \le \cdots \le x_d \le 1\} \subset \mathbb{R}^d.$$

(If not, a linear map onto any other $d$-dimensional simplex will have constant Jacobean and therefore will convert any uniform distribution on this simplex to a uniform distribution on the target.)

Start anywhere at $\mathbf{x}^{(0)}\in\mathbb{R}^d$. Generate independent, identically distributed increments $\mathbf{h}^{(i)} \in \mathbb{R}^d$, $i=1, 2, \ldots, n$, according to your random walk. Compute their cumulative sums

$$\mathbf{x}^{(i)} = \mathbf{x}^{(i-1)} + \mathbf{h}^{(i)}$$

for $i=1, 2, \ldots, n$.

Consider the reflections through all hyperplanes that are parallel to the coordinate planes and situated at integral distances from them. These generate a reflection group $H$ that includes translations by vectors all of whose coordinates are integral multiples of $2$. A fundamental domain for $H$ is the unit cube

$$I^d = \{(x_1, \ldots, x_d)\mid 0 \le x_i \lt 1, i=1, 2, \ldots, d\}.$$

Replacing each $\mathbf{x}^{(i)}$ by its unique representative $\mathbf{y}^{(i)}$ in $I^d$ produces a random walk on $I^d$ that asymptotically is uniform. (Indeed, when the standard deviation of the coordinate increments is not too small compared to $1$, the uniformity is very rapidly achieved: see method (8) at https://stats.stackexchange.com/a/117711/919.)

Figure showing scatterplot matrix of y for d=4

Here is a scatterplot matrix of the components of $\mathbb{y}$ for $d=4$, $n=10^4$. Increments were drawn from a Normal distribution with standard deviation $d$.

This random walk projects to a uniform walk on the simplex $\Delta^d$ simply by sorting each vector $\mathbb{y}^{(i)}$. Sorting can be viewed as including $H$ in a larger reflection group $G$ whose generators also contain all permutations of the coordinates.

Figure showing scatterplot matrix of z

The projections onto planes do not look uniform--they shouldn't be.

Please note that the individual coordinates cannot be uniform or even identically distributed and that their increments will not be identically distributed either.

Histograms of z coordinates

These are histograms of individual coordinates ($z_1$ through $z_4$ are shown left to right).

Histograms of increments of z

These are histograms of the increments of the individual coordinates. Note how the middle histograms differ from the outer ones.

For the purpose of generating proposals in an MCMC algorithm, you may generate the original walk $\mathbf{x}^{(i)}$ sequentially: the new proposal is its reduction modulo $G$.


The description of the group $G$ in terms of two sets of reflections allows realizations of these walks to be efficiently generated: beginning with a walk in $\mathbb{R}^d$, reduce the coordinates modulo $H$ and then reduce the results by sorting the coordinates. To reduce coordinates modulo $H$, exploit the fact that $H$ contains all the translations by multiples of $2$: reduce all coordinates modulo $2$, and then apply the reflection

$$u \to 2-u$$

for any coordinate value $1\le u \lt 2$. Finally, sort the coordinates to obtain the corresponding point in the simplex.

Here is working R code to illustrate.

#
# Uniform random walk on the simplex 0 <= x[1] <= ... <= x[d] <= 1,
#
d <- 4
n <- 1e4          # Number of points
sigma <- d        # SD of lengths Normal-component increments
start <- rep(0, d)
#
# Generate the walk.
#
dx <- matrix(rnorm(n*d, sd=sigma/sqrt(d)), nrow=n) # n X d matrix of increments
x <- apply(dx, 2, cumsum)                          # n X d matrix of raw values
y <- x %% 2
y <- apply(y, 2, function(u) ifelse(u >= 1, 2-u, u))
z <- t(apply(y, 1, sort))
#
# Display the walk `z` in various ways,
# first with scatterplot matrices of `y` and `z`.
#
colnames(y) <- paste("y", 1:d, sep=".")
colnames(z) <- paste("z", 1:d, sep=".")
pairs(y, cex=1/4, col="#00000020")
pairs(z, cex=1/4, col="#00000020")
#
# Histograms of `z` and its increments `dz`.
#
dz <- apply(z, 2, diff)
par(mfrow=c(1, d))
invisible(apply(z, 2, hist, freq=FALSE, breaks=30, xlab="z", main=""))
invisible(apply(dz, 2, hist, freq=FALSE, breaks=30, xlab="dz", main=""))
$\endgroup$
6
  • $\begingroup$ Thanks. I need $\sum_{i=1}^{d}x_i =1$; it is not clear to me that this is a constraint in your algorithm? $\endgroup$
    – Travis
    Sep 12 '17 at 8:29
  • $\begingroup$ I appreciate your effort. I will implement this and have a look. $\endgroup$
    – Travis
    Sep 12 '17 at 9:48
  • $\begingroup$ @whuber: brilliant! $\endgroup$
    – Xi'an
    Sep 12 '17 at 11:26
  • $\begingroup$ Travis, my initial remarks address this constraint (assuming your $x_i$ are nonnegative). The transformation from the simplex I use to the simplex you want is to take the successive differences of the coordinates. This is a one-to-one area-preserving correspondence between the two. $\endgroup$
    – whuber
    Sep 12 '17 at 14:04
  • $\begingroup$ Can you please show (by connecting the dots) the random walk? $\endgroup$
    – Travis
    Sep 13 '17 at 15:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.