Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Start at a point $(0,0,z_0)$ and take $n$ steps of unit length in a random direction (for each step) in $\mathbb{R}^3$. Let such a walk be valid if the position of the last step, and only the last step, is under the plane $(x_n,y_n,z_n<0)$. Let a walk be almost valid if at least the last step is under the plane.

How can I sample uniformly from the set of all valid or almost valid walks without a trial and rejection method?

share|cite|improve this question
Metropolis algorithm? – mjqxxxx Oct 10 '12 at 18:59
@mjqxxxx How would I calculate an acceptance ratio? Generating the weight of a partial walk is part of the problem I think. – Hooked Oct 10 '12 at 19:03
Nice question. Are you using "step" and "point" interchangeably, or is that part of the distinction between valid and almost valid? – joriki Oct 10 '12 at 19:12
@joriki For $n$ steps, there are $n+1$ points on the walk (counting the first and last). I clarified the wording a bit to make that clear. In an almost valid walk, the points can dip below the plane before the $n$th step. In a valid walk only the last can be below the plane. – Hooked Oct 10 '12 at 19:24

The thing to note when you're uniformly sampling directions from the unit sphere is that in three dimensions that means that $z$ is distributed uniformly. So you can ignore the sideways motion and effectively have a one-dimensional random walk with steps uniformly distributed in $[-1,1]$ (after scaling the step length to $1$).

The proportion of valid walks for given $n$ as a function of $z$ is piecewise polynomial of degree at most $n$, with one polynomial for each interval between two integers. If $n$ isn't too large, you can compute these polynomials, then do a not-too-wasteful form of rejection sampling by bounding them by a constant or linear or quadratic function over $[z_0-1,z_0+1]$. It will take a bit of bookkeeping to treat the boundary at $z=0$ correctly, but it shouldn't be too difficult.

Here's an idea which I think will only work for the almost valid case, but it might work well for that.

The problem is that the piecewise polynomial distributions get increasingly complicated when you convolve them with the uniform distribution on $[-1,1]$. So we'd like to replace them with something for which the convolution with the uniform distribution is easy to compute and doesn't increase the complexity. A sum of sinusoidals fits that bill, since each sinusoidal is just multiplied by a constant by the convolution.

Start from the end with a constant density on $[z_0-n,0]$. Convolve this with the uniform distribution on $[-1,1]$ some number of times, perhaps $100$ or $1000$, as many as you can afford to handle the polynomials for. Now it strongly resembles a Gaussian. Choose a sufficiently large interval centred on the mean of the distribution at $(z_0-n)/2$ such that both the bulk of the distribution and any $z$ values you expect to reach are well inside the interval, compute a partial sum of the Fourier series for the distribution on this interval, bound the magnitude of the remaining tail and add a corresponding constant to the partial sum to arrive at a Fourier sum that approximates the distribution and bounds it from above. Since the near-Gaussian is quite smooth, you won't need too many terms to get a good approximation.

This Fourier sum is now your new target. You can very easily propagate it to $n$, even for $n\gg1000$, since the convolution is just a multiplication. Start with $(n,z_0)$, pretending that you're trying to sample from the target distribution. Then when you reach the stage of the target distribution, do a single rejection step in which you reject the generated walk if a number drawn uniformly between $0$ and the Fourier sum lies above the correct distribution. This gives you the right distribution, and then you can proceed with the correct polynomials.

I wrote the answer below when I misread the question to refer to random walks on a lattice; I'm letting it stand since it might be interesting to others in its own right.

I like the question. (I used to do lattice gauge theory simulations.)

To do this, you need to know the number $a(n,z)$ of (almost) valid walks.

If you shift $z$ by $1$ so that a walk is valid if it ends at $z=0$, which simplifies the indexing slightly, the numbers $a(n,z)$ for valid walks are given by OEIS sequence A052179. The entry gives several prescriptions for generating them. I think the most efficient one for you is probably the recurrence relation,

$$ a(n,z)=a(n-1,z-1)+4a(n-1,z)+a(n-1,z+1)\;, $$

since you only need to compute $a(n,z_0)$ once at the beginning and store the intermediate results, which will contain all values you need later in the simulation.

Once you know the values $a(n-1,z-1)$, $a(n-1,z)$ and $a(n-1,z+1)$ at your current point, all you have to do is draw a random number uniformly between $1$ and $a(n,z)$ and choose where to go according to which of the six intervals of lengths $a(n-1,z_0-1)$, $a(n-1,z_0+1)$ and four times $a(n-1,z_0)$ it falls into.

The recurrence relation is the same for almost valid walks; only the initial conditions are different.

share|cite|improve this answer
Thanks for the input joriki. You are assuming that the steps are only on a cubic lattice correct? Is it possible to apply the same treatment to a face-centered cubic or an off-lattice step? – Hooked Oct 10 '12 at 20:03
@Hooked: I see. Sorry, I misread the question -- I should have realized that this was why it said "off-lattice" in the title :-) The face-centered cubic lattice shouldn't be a problem; you'll just get a different and slightly more complicated recurrence; but continuous directions uniformly sampled from all unit vectors might pose more of a challenge. – joriki Oct 10 '12 at 20:34
I figured as much - my usual line of attack for these types of problem with they are discrete is to generate small cases and try OEIS too. There is no counterpart for continuous problems that I know of however. For context, this problem comes from a polymer walk, starting at a fixed binding site with a fixed length and ending at a membrane (the plane). I can always create the walks and reject, but the attrition rate gets high when $nL$ is about the same order as $z$. – Hooked Oct 10 '12 at 20:50
@Hooked: OK, I gave a suggestion for the unit sphere case. Whether this is practical will depend on the values of the parameters. An advantage of this approach might be that when $nL$ is close to $z$ you don't actually need most of the polynomials (whereas for $nL\gg z$ the number of coefficients to be computed is quadratic in $n$). – joriki Oct 10 '12 at 21:20
@Hooked: I added a new idea in case you need high $n$ and can't afford to handle all the polynomials. I'd be interested to know whether this might be useful. – joriki Oct 10 '12 at 23:12

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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