How to implement a uniform random walk on a simplex? 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. 

 A: 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.) 

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.

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.

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

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=""))

