Hitting Probabilities for Markov Chains

Posted on June 17, 2014 by dominicyeo

This continues my previous post on popular questions in second year exams. In the interest of keeping it under 2,500 words I’m starting a new article.

In a previous post I’ve spoken about the two types of Markov chain convergence, in particular, considering when they apply. Normally the ergodic theorem can be used to treat the case where the chain is periodic, so the transition probabilities do not converge to a stationary distribution, but do have limit points – one at zero corresponding to the off-period transitions, and one non-zero. With equal care, the case where the chain is not irreducible can also be treated.

A favourite question for examiners concerns hitting probabilities and expected hitting times of a set A. Note these are unlikely to come up simultaneously. Unless the hitting probability is 1, the expected hitting time is infinite! In both cases, we use the law of total probability to derive a family of equations satisfied by the probabilities/times. The only difference is that for hitting times, we add +1 on the right hand side, as we advance one time-step to use the law of total probability.

The case of hitting probabilities is perhaps more interesting. We have:

$h_i^A = 1,\; i\in A, \quad h_i^A=\sum_{j\in S}p_{ij}h_j^A,\; i\not\in A.$

There are two main cases of interest: where the chain is finite but has multiple closed communicating classes, and where the chain is infinite, so even though it is irreducible, a trajectory might diverge before hitting 0.

For the case of a finite non-irreducible Markov chain, this is fairly manageable, by solving backwards from states where we know the values. Although of course you could ask about the hitting probability of an open state, the most natural question is to consider the probability of ending up in a particular closed class. Then we know that the hitting probability starting from site in the closed class A is 1, and the probability starting from any site in a different closed class is 0. To find the remaining values, we can work backwards one step at a time if the set of possible transitions is sparse enough, or just solve the simultaneous equations for $\{h_i^A: i\text{ open}\}$ .

We therefore care mainly about an infinite state-space that might be transient. Typically this might be some sort of birth-and-death chain on the positive integers. In many cases, the hitting probability equations can be reduced to a quadratic recurrence relation which can be solved, normally ending up with the form

$h_i=A+B\lambda^i$ ,

where $\lambda$ might well be q/p or similar if the chain is symmetric. If the chain is bounded, typically you might know $h_0=1, h_N=0$ or similar, and so you can solve two simultaneous equations to find A and B. For the unbounded case you might often only have one condition, so you have to rely instead on the result that the hitting probabilities are the minimal solution to the family of equations. Note that you will always have $h^i_i=1$ , but with no conditions, $h^i_j\equiv 1$ is always a family of solutions.

It is not clear a priori what it means to be a minimal solution. Certainly it is not clear why one solution might be pointwise smaller than another, but in the case given above, it makes sense. Supposing that $\lambda<1$ , and A+B=1 say, then as we vary the parameters, the resulting set of ‘probabilities’ does indeed vary monotically pointwise.

Why is this true? Why should the minimum solution give the true hitting probability values? To see this, take the equations, and every time an $h_i^A$ appears on the right-hand side, substitute in using the equations. So we obtain, for $i\not\in A$ ,

$h_i^A=\sum_{j\in A}p_{ij}+\sum_{j\not\in A} p_{ij}h_j^A,$

and after a further iteration

$h_i^A=\sum_{j_1\in A}p_{ij_1}+\sum_{j_1\not\in A, j_2\in A}p_{ij_1}p_{j_1j_2}+\sum_{j_1,j_2\not\in A}p_{ij_1}p_{j_1j_2}h_{j_2}^A.$

So we see on the RHS the probability of getting from i to A in one step, and in two steps, and if keep iterating, we will get a large sum corresponding to the probability of getting from i to A in 1 or 2 or … or N steps, plus an extra term. Note that the extra term does not have to correspond to the probability of not hitting A by time N. After all, we do not yet know that $(h_{i}^A)$ as defined by the equations gives the hitting probabilities. However, we know that the probability of hitting A within N steps converges to the probability of hitting A at all, since the sequence is increasing and bounded, so if we take a limit of both sides, we get $h_i^A$ on the left, and something at least as large as the hitting probability starting from i on the right, because of the extra positive term. The result therefore follows.

It is worth looking out for related problems that look like a hitting probability calculation. There was a nice example on one of the past papers. Consider a simple symmetric random walk on the integers modulo n, arranged clockwise in a circle. Given that you start at state 0, what is the probability that your first return to state 0 involves a clockwise journey round the circle?

Because the system is finite and irreducible, it is not particularly interesting to consider the actual hitting probabilities. Also, note that if it is convenient to do so, we can immediately reduce the problem when n is even. In two steps, the chain moves from j to j+2 and j-2 with probability ¼ each, and stays at j with probability ½. So the two step chain is exactly equivalent to the lazy version of the same dynamics on n/2.

Anyway, even though the structure is different, our approach should be the same as for the hitting probability question, which is to look one step into the future. For example, to stand a chance of working, our first two moves must both be clockwise. Thereafter, we are allowed to move anticlockwise. There is nothing special about starting at 0 in defining the original probability. We could equally well ask for the probability that starting from j, the first time we hit 0 we have moved clockwise round the circle.

The only thing that is now not obvious is how to define moving clockwise round the circle, since it is not the case that all the moves have to be clockwise to have experienced a generally clockwise journey round the circle, but we definitely don’t want to get into anything complicated like winding numbers! In fact, the easiest way to make the definition is that given the hitting time of 0 is T, we demand that the chain was at state n at time T-1.

For convenience (ie to make the equations consistent) we take $h_0=0, h_n=1$ in an obvious abuse of notation, and then

$h_j=\frac12h_{j-1}+\frac12 h_{j+1},$

from which we get

$h_j=a+bj \Rightarrow h_j=\frac{j}{n}.$

Of course, once we have this in mind, we realise that we could have cut the circle at 0 (also known as n) and unfolded it to reduce the problem precisely to symmetric gambler’s ruin. In particular, the answer to the original problem is 1/2n, which is perhaps just a little surprising – maybe by thinking about the BM approximation to simple random walk, and that BM started from zero almost certainly crosses zero infinitely many times near we might have expected this probability to decay faster. But once it is unfolded into gambler’s ruin, we have the optimal stopping martingale motivation to reassure us that this indeed looks correct.

Avoiding Mistakes in Probability Exams

Posted on June 16, 2014 by dominicyeo

Over the past week, I’ve given several tutorials to second year undergraduates preparing for upcoming papers on probability and statistics. In particular, I’ve now seen a lot of solutions to a lot of past papers and specimen questions, and it’s worthwhile to consider some of the typical mistakes students can make on these questions. Of course, as with any maths exam, there’s always the possibility of a particularly subtle or involved question coming up, but if the following three common areas of difficulty can be avoided, you’re on track for doing well.

Jacobians

In a previous course, a student will learn how to calculate the pdf of a function of a random variable. Here, we move onto the more interesting and useful case of finding the (joint) density of function(s) of two or more random variables. The key thing to remember here is that manipulating pdfs is not a strange arbitrary exercise – it is just integration. It is rarely of interest to consider the value of a pdf at a single point. We can draw meaningful conclusions from a pdf or from comparison of two pdfs by integrating them.

Then the question of substituting for new random variables is precisely integration by substitution, which we are totally happy with in the one-dimensional case, and should be fairly happy with in the two-dimensional case. To get from one joint density to another, we multiply by the absolute value of the Jacobian. To ensure you get it right, it makes sense to write out the informal infinitesimal relation

$f_{U,V}(u,v) du dv = f_{X,Y}(x,y)dx dy.$

This is certainly relevant if we put integral signs in front of both sides, and explains why you obtain $f_{U,V} = \frac{d(x,y)}{d(u,v)} f_{X,Y}$ rather than the other way round. Note though that if $\frac{d(u,v)}{d(x,y)}$ is easier to calculate for some reason, then you can evaluate this and take the inverse, as your functions will almost certainly be locally bijective almost everywhere.

It is important to take the modulus of the Jacobian, since densities cannot be negative! If this looks like a fudge, then consider the situation in one dimension. If we substitute for $x\mapsto f(x)=1-x$ , then f’ is obviously negative, BUT we also end up reversing the order of the bounds of the integral, eg [1/3, ¾] will become [2/3,1/4]. So we have a negative integrand (after multiplying by f'(x)) but bounds in the ‘wrong’ order. These two factors of -1 will obviously cancel, so it suffices just to multiply by |f'(x)| at that stage. It is harder to express in words, but a similar relation works for the Jacobian substitution.

You also need to check where the new joint density is non-zero. Suppose X, Y are supported on [0,1], then when we write $f_{X,Y}(x,y)$ we should indicate that it is 0 off this region, either by splitting into cases, or adding the indicator function $1_{\{x,y\in[0,1]\}}$ as a factor. This is even more important after substitutions, as the range of the resulting random variables might be less obvious than the originals. Eg with X,Y as above, and $U=X^2, V=X/Y$ , the resulting pdf will be non-zero only when $u\in[0,1], v\ge \sqrt{u}$ . Failing to account for this will often lead to ludicrous answers. A general rule is that you can always check that any distribution you’ve produced does actually integrate to one.

Convergence using MGFs

There are two main reasons to use MGFs and PGFs. The first is that they behave nicely when applied to (possibly random) sums of independent random variables. The independence property is crucial to allow splitting of the MGF of the sum into the product of MGFs of the summands. Of course, implicit in this argument is that MGFs determine distributions.

A key theorem of the course is that this works even in the limit, so you can use MGFs to show convergence in distribution of a family of distributions. For this, you need to show that the MGFs converge pointwise on some interval [-a,a] around 0. (Note that the moments of the distribution are given by the family of derivatives at 0, as motivation for why this condition might be necessary.) Normally for such questions, you will have been asked to define the MGF earlier in the question, and probably will have found the MGF of a particular distribution or family of distributions, which might well end up appearing as the final answer.

Sometimes such an argument might involve substituting in something unusual, like t/N, rather than t, into a known MGF. Normally a Taylor series can be used to show the final convergence result. If you have a fraction, try to cancel terms so that you only have to evaluate one Taylor series, rather than lots.

Using the Markov Property

The Markov property is initially confusing, but once we become comfortable with the statement, it is increasingly irritating to have to answer the question: “show that this process has the Markov property.” This question is irritating because in most cases we want to answer: “because it obviously does!” Which is compelling, but unlikely to be considered satisfactory in a mathematics exam. Normally we observe that the random dynamics of the next step are a function only of the present location. Looking for the word ‘independent’ in the statement of the process under discussion is a good place to start for any argument along these lines.

The most developed example of a Markov process in this course is the Poisson process. I’ve written far too much about this before, so I won’t do so again, except to say this. When we think of the Poisson process, we generally have two thoughts going through our minds, namely the equivalent definitions of IID exponential inter-arrival times, and stationary, Poisson increments (or the infinitesimal version). If we draw a sketch of a sample trajectory of this process, we can label everything up and it is clear how it all fits together. But if you are asked to give a definition of the Poisson process $(N_t)$ , it is inappropriate to talk about inter-arrival times unless you define them in terms of $N_t$ , since that is the process you are actually trying to define! It is fine to write out

$T_k:=\min\{t: N_t=k\},\quad N_t=\max\{k: Y_1+Y_2+\ldots+Y_k\le t\}$

but the relation between the two characterisations of the process is not obvious. That is why it is a theorem of the course.

We have to be particularly careful of the difference in definition when we are calculating probabilities of various events. A classic example is this. Find the distribution of $N_2$ , conditional on $T_3=1$ . It’s very tempting to come up with some massive waffle to argue that the answer is 3+Po(1). The most streamlined observation is that the problem is easy if we are conditioning instead on $N_1=3$ . We just use the independent Poisson increments definition of $(N_t)$ , with no reference to inter-arrival times required. But then the Markov property applied at time 1 says that the distribution of $(N_2)$ depends only on the value of $N_1$ , not on the process on the interval [0,1). In a sense, the condition that $T_3=1$ is giving us extra information on the behaviour of the process up to time 1, and the Markov property, which we know holds for the Poisson process, asserts precisely that the extra information doesn’t matter.

Popoviciu’s Inequality

Posted on May 29, 2014 by dominicyeo

I’ve just returned to the UK after an excellent stay at the University of British Columbia. More about that will follow in some posts which are being queued. Anyway, I flew back in time to attend the last day of the camp held at Oundle School to select the UK team for this year’s International Mathematical Olympiad, to be held in Cape Town in early July. I chose to give a short session on inequalities, which is a topic I did not enjoy as a student and do not enjoy now, but perhaps that makes it a particularly suitable choice?

We began with a discussion of convexity. Extremely occasionally in olympiads, and merely slightly occasionally in real life, an inequality arises which can be proved by showing that a given function is convex in all its arguments, hence its maximum must be attained at a boundary value in each variable.

In general though, our main experience of convexity will be through the medium of Jensen’s inequality. A worthwhile check is to consider one form of the statement of Jensen’s inequality, with two arguments. We are always given a convex function f defined on an interval I=[a,b], and $x,y\in I$ , and weights $\alpha,\beta$ which sum to 1. Then

$\alpha f(x)+\beta f(y)\ge f(\alpha x+\beta y).$

How do we prove this? Well, in fact this is the natural definition of convexity for a function. There had initially been vague murmurings that convexity should be defined as a property of the second derivative of the function. But this is somewhat unsatisfactory, as the function $f(x)=|x|$ is certainly convex, but the second derivative does not exist at x=0. One could argue that the second derivative may not be finite at x=0, but is nonetheless positive by defining it as a limit which happens to be infinite in this case. However, I feel it is uncontroversial to take the case of Jensen given above as the definition of convexity. It is after all a geometric property, so why raise objections to a geometric definition?

The general statement of Jensen’s inequality, with the natural definitions, is

$\sum_{i} \alpha_i f(x_i)\ge f(\sum_{i}\alpha_ix_i).$

This is sometimes called Weighted Jensen in the olympiad community, with ‘ordinary’ Jensen following when the weights are all 1/n. In a probabilistic context, we write

$\mathbb{E}[f(X)]\ge f(\mathbb{E}X),$

for X any random variable supported on the domain of f. Naturally, X can be continuous as well as discrete, giving an integral version of the discretely weighted statement.

Comparing ‘ordinary’ Jensen and ‘weighted’ Jensen, we see an example of the situation where the more general result is easier to prove. As is often the case in these situations, this arises because the more general conditions allow more ‘elbow room’ to perform an inductive argument. A stronger statement means that assuming the induction hypothesis is more useful! Anyway, I won’t digress too far onto the proof of discrete ‘weighted’ Jensen as it is a worthwhile exercise for olympiad students.

What I wanted to discuss principally was an inequality due to Tiberiu Popoviciu:

$\frac13[f(x)+f(y)+f(z)]+f(\frac{x+y+z}{3})\ge \frac23[f(\frac{x+y}{2})+f(\frac{y+z}{2})+f(\frac{z+x}{2})].$

We might offer the following highly vague intuition. Jensen asserts that for sums of the form $\sum f(x_i)$ , you get larger sums if the points are more spread out. The effect of taking the mean is immediately to bring all the points as close together as possible. But Popoviciu says that this effect is so pronounced that even with only half the weight on the outer points (and the rest as close together as possible), it still dominates a system with the points twice as close together.

So how to prove it? I mentioned that there is, unsurprisingly, a weighted version of this result, which was supposed to act as a hint to avoid getting too hung up about midpoints. One can draw nice diagrams with a triangle of points (x,f(x)), (y,f(y)), (z,f(z)) and draw midpoints, medians and centroids, but the consensus seemed to be that this didn’t help much.

I had tried breaking up the LHS into three symmetric portions, and using weighted Jensen to obtain terms on the RHS, but this also didn’t yield much, so I warned the students against this approach unless they had a specific reason to suppose it might succeed.

Fortunately, several of the students decided to ignore this advice, and though most fell into a similar problem I had experienced, Joe found that by actively avoiding symmetry, a decomposition into two cases of Jensen could be made. First we assume WLOG that $x\le y \le z$ , and so by standard Jensen, we have

$\frac13[f(x)+f(y)]\ge \frac23 f(\frac{x+y}{2}).$

It remains to show

$\frac13 f(z)+f(\frac{x+y+z}{3})\ge \frac23[f(\frac{x+z}{2})+f(\frac{y+z}{2})].$

If we multiply by ¾, then we have an expression on each side that looks like the LHS of Weighted Jensen. At this point, it is worth getting geometric again. One way to visualise Jensen is that for a convex function, a chord between two points on the function lies above the function. (For standard Jensen with two variables, in particular the midpoint lies above the function.) But indeed, suppose we have values $x_1<x_2<y_2<y_1$ , then the chord between $f(x_1),f(y_1)$ lies strictly above the chord between $f(x_2),f(y_2)$ . Making precisely such a comparison gives the result required above. If you want to be more formal about it, you could consider replacing the values of f between $x_2,y_2$ with a straight line, then applying Jensen to this function. Linearity allows us to move the weighting in and out of the brackets on the right hand side, whenever the mean lies in this straight line interval.

Hopefully the diagram above helps. Note that we can compare the heights of the blue points (with the same abscissa), but obviously not the red points!

In any case, I was sceptical about whether this method would work for the weighted version of Popoviciu’s inequality

$\alpha f(x)+\beta f(y) + \gamma f(z)+f(\alpha x+\beta y+\gamma z)\ge$

$(\alpha+\beta)f(\frac{\alpha x+\beta y}{\alpha+\beta})+(\beta+\gamma)f(\frac{\beta y + \gamma z}{\beta+\gamma})+(\gamma+\alpha)f(\frac{\gamma z+\alpha x}{\gamma+\alpha}).$

It turns out though, that it works absolutely fine. I would be interested to see a solution to the original statement making use of the medians and centroid, as then by considering general Cevians the more general inequality might follow.

That’s all great, but my main aim had been to introduce one trick which somewhat trivialises the problem. Note that in the original statement of Popoviciu, we have a convex function, but we only evaluate it at seven points. So for given x,y,z, it makes no difference if we replace the function f with a piece-wise linear function going through the correct seven points. This means that if we can prove that the inequality for any convex piece-wise linear function with at most eight linear parts then we are done.

(There’s a subtlety here. Note that we will prove the inequality for all such functions and all x,y,z, but we will only use this result when x,y,z and their means are the points where the function changes gradient.)

So we consider the function

$g_a(x)=\begin{cases}0& x\le 0\\ ax & x\ge 0\end{cases}$

for some positive value of a. It is not much effort to check that this satisfies Popoviciu. It is also easy to check that the constant function, and the linear function g(x)=bx also satisfy the inequality. We now prove that we can write the piece-wise linear function as a sum of functions which satisfy the inequality, and hence the piece-wise linear function satisfies the inequality.

Suppose we have a convex piecewise linear function h(x) where $x_1<x_2<\ldots<x_n$ are the points where the derivative changes. We write

$a_i=h'(x_i+)-h'(x_i'),\quad a_0=h'(x_1-),$

for the change in gradient of h around point $x_i$ . Crucially, because h is convex, we have $a_i\ge 0$ . Then we can write h as

$h(x)=C+ a_0x+g_{a_1}(x-x_1)+\ldots+g_{a_{n}}(x-x_n),$

for a suitable choice of the constant C. This result comes according to [1] as an intermediate step in a short paper of Hardy, Littlewood and Polya, which I can’t currently find online. Note that inequalities are preserved under addition (but not under subtraction) so it follows that h satisfies Popoviciu, and so the original function f satisfies it too for the values of x,y,z chosen. These were arbitrary (but were used to construct h), and hence f satisfies the inequality for all x,y,z.

Some further generalisations can be found in [1]. With more variables, there are more interesting combinatorial aspects that must be checked, regarding the order of the various weighted means.

[1] – D. Grinberg – Generalizations of Popoviciu’s Inequality. arXiv

A Month in Vancouver

Posted on May 22, 2014 by dominicyeo

I’ve spent the past ~3 weeks at the University of British Columbia in Vancouver, working on some projects related to forest fires and frozen percolation with Balazs Rath. In advance, I’d imagined that this might be a prolific time for new articles on this blog, but in the end it was more like the exact opposite. After spending all of every day playing with technical details on our models (which do not lend themselves naturally to informal writing…) if I felt like doing maths in the evenings, it was more important to keep working on the questions in hand.

Anyway, I had a great time mathematically and otherwise, so I feel a shameless tourism post is justified to report what I enjoyed most about Vancouver.

Asian Food

It sounds ridiculous to return from a trip and reference this as the highlight, but the most noticeable long-term result of a month in Vancouver is that I will probably never be able to order sushi in the UK again. The typical house roll cost about £3 and contained so much seafood that it was essentially impossible to manipulate. Or maybe that was just my (lack of) chopstick skills? Anyway, as I’ve observed in other parts of Canada, groceries bizarrely expensive relative to eating out, at least relative to Europe. By my third week I had basically given up and found a different East Asian cuisine every other night. I ate Japanese, Chinese, Korean, Thai, Malaysian, Taiwanese, Vietnamese and others that didn’t feel constrained by a single national identity. All were uniformly excellent, unfailingly welcoming, and made commendable efforts to include smoked salmon in every single dish.

Outdoor Pursuits

For a country where a donut shop is an integral part of the national identity, the total absence of fat people was immediately noticeable. In their place was a city where everyone seemed to go running, hiking, skiing, quad-biking, skidoo-ing and so on at the slightest opportunity. I haven’t seen that much lycra since breakfast in college on the morning of May Bumps…

For a major world city, the wilderness was tantalisingly close, and remarkably accessible. The public transit system spread far into the mountains and up the coast, and offered you any 90 minutes worth of travel for $2. A breezy ferry trip across the sound to Bowen Island offered the chance to hike up Mt Gardner, an unrelenting and totally solitary trek up through the pine forest. I was conscious that if I slipped somewhere, there wouldn’t be much left to find once the bears had got to me, but the views over the inlet from the summit made it all worthwhile. A trip to Mt Seymour on Victoria Day, then down past the vertigo-inducing Lynn Suspension bridge and a trail back towards civilisation was equally successful. On both occasions, felt I’d very much earned that donut.

UBC

Of course, the main reason for the trip was to do worthwhile mathematics, and I certainly haven’t experienced a nicer university campus to do exactly that. A five minute walk and 500 steps away from the Wysteria-covered maths department is the 7km long Wreck Beach, a sunny haven for Vancouver’s hipsters and nudists, protected from direct exposure to the Pacific by Vancouver Island. Of more interest to me were the 30km or so of trails through the forest in the Pacific Spirit Regional Park, directly across the road from where I was staying in Wesbrook Village at the south end of the campus.

Most importantly, I’m extremely grateful for the invitation to work at UBC for these past few weeks. It’s been productive, and an excellent opportunity to meet some new people and learn about some interesting aspects of probability. I very much hope an opportunity will arise to go back soon!

Coupling from the Past

Posted on April 30, 2014 by dominicyeo

In a long series of previous posts I have talked about mixing times for Markov chains. We consider how long it takes for the distribution of a particular Markov chain to approach equilibrium. We are particularly interested in the asymptotics when some parameter of the model grows, such as the size of the state space, grows to infinity.

But why are we interested in the underlying problem? The idea of Markov Chain Monte Carlo methods is to sample from an intractable distribution by instead sampling from a Markov chain which approximates the distribution well at large times. A distribution might be intractable because it is computationally demanding to work out the normalising constant, or it might be distributed uniformly on a complicated combinatorial set. If, however, the distribution is the equilibrium distribution of some Markov chain, then we know how to at least sample from a distribution which is close to the one we want. But we need to know how long to run the process. We will typically tolerate some small error in approximating the distribution (whether we measure this in terms of total variation distance or some other metric doesn’t really matter at this heuristic level), but we need to know how it scale. If we double the size of the system, do we need to double the number of iterations of the chain, or square it. This is really important if we are going to use this for large real-world models with finite computing power!

Sometimes though, an approximation is not enough. If we want an exact sample from the equilibrium distribution, Markov chains typically will not help us as it is only in very artificial examples that the distribution after some finite time is actually the equilibrium distribution. One thing that we might use is a stationary time, which is a stopping time T, for which $X_T\stackrel{d}{=}\pi$ . Note that there is one trivial way to do this. We can sample Y from distribution $\pi$ before starting the process, then stop X at the first time T for which $X_T=Y$ . But this is no help really, as we need to have Y in the first place!

So we are really interested in less trivial stationary times. Perhaps the best example is the top-to-random shuffle. Here we are given a pack of labelled cards, WLOG initially in descending order at each step we move the top card in the pile to a randomly-chosen location in the pile (which includes back onto the top). Then it turns out that the first time we move the card originally at the bottom from the top to somewhere is a strong stationary time. This is fairly natural, as by this time, every card has been involved in at least one randomising event.

Anyway, so this gives a somewhat artificial way to sample from the uniform distribution on a pack of cards. This strong stationary time is almost surely finite, with distribution given by the coupon collector problem, for which the expectation grows as $n\log n$ , where n is the number of cards.

The problem with this method is that it is not easy in general to come up with a non-contrived stationary time such as this one. The idea of coupling from the past, discussed by some previous authors but introduced in this context by Propp and Wilson in the mid ’90s, is another method to achieve perfect sampling from the equilibrium distribution of a Markov chain. The idea here is to work backwards rather than forwards. The rest of this post, which discusses this idea, is based on the talk given at the Junior Probability Seminar by Irene, and on the chapter in the Levin, Peres, Wilmer book.

The key to the construction is a coupling of the transitions of a Markov chain. In the setting of a simple random walk, we have by construction a coupling of the transitions. It doesn’t matter which state we are at: we toss a coin to decide whether to move up or down, and we can do this without reference to our current position. Levin, Peres and WIlmer call this a random mapping representation in general, and it is yet another concept that is less scary than its definition might suggest.

Given a transition matrix P on state space S, such a representation is a function

$\phi: S\times[0,1]\rightarrow S,\text{ s.t. }\mathbb{P}(\phi(i,U)=j)=p_{ij},$

where U is a U(0,1) random variable independent of choice of i. In particular, once we have the random value of u, we can consider $\phi(i,u)$ as i varies, to obtain a random map $S\rightarrow S$ . Crucially, this map is not necessarily a bijection.

Note first that there are many possibilities for constructing the representation $\phi$ . For some chains, and some representations, in particular random walks on vertex-transitive graphs (such as SRW – only for now we are restricting attention to finite state spaces) it is possible to choose $\phi$ so that it always gives a bijection, but it is also always possible to choose it so that there is some probability it doesn’t give a bijection.

Let $U_1,U_2,\ldots$ be an IID sequence of U[0,1] random variables, and write $\phi_i$ for the random map induced by $U_i$ . Then consider the sequence of iterated maps:

$\phi_1, \phi_1\circ \phi_2, \ldots, \phi_1\circ\ldots\circ\phi_n,$

and let T be the (random) smallest time such that the image of $\phi_1\circ\ldots\circ \phi_T$ is a single state. Ie, as we go backwards in time through the maps $\phi_i$ , we are gradually losing various states, corresponding to the maps not being bijections. Since the state space is finite, and the probability of not being a bijection is positive, it can be shown that T is almost surely finite. The claim then is that

$Y=\text{Im}(\phi_1\circ\ldots\circ \phi_T)$

is distributed as the equilibrium distribution of the chain. We finish by proving this.

Proof: Since the algorithm terminates after finite time almost surely, given any $\epsilon>0$ , we can choose N such that the probability the algorithm stops in at most N steps is greater than $1-\epsilon$ .

Now run the Markov chain from time -N, started in the equilibrium distribution, with the transition from time -t to -(t-1) given by the random mapping driven by $U_t$ . Thus at time 0, the distribution of the chain is still the equilibrium distribution. But if we condition on the event that $T\le N$ , then $X_0=\phi_1\circ \ldots \circ\phi_n(X_{-N})=Y$ regardless of the initial value. So $\mathbb{P}(X_0\ne Y)<\epsilon$ , and hence the result follows, since $\epsilon>0$ was arbitrary.

What makes this easier than strong stationary times is that we don’t have to be clever to come up with the stopping time. It is however still important to know how long on average it takes to run the algorithm. At the end of her talk, Irene showed how to adapt this algorithm to deal with Probabilistic Cellular Automata. Roughly speaking, these are a sequence of infinite strings of 0s and 1s. The value of some element is determined randomly as a function of the values in the row underneath, say the element directly underneath and the two either side. In that setting, if you start with a finite subsequence and couple from the past by looking down to lower rows, each time you drop down a row you consider one further element, so in fact the coupling from the past algorithm has to eliminate possibilities fast enough to make up for this, if we want to terminate almost surely in finite time.

Here’s a link to the paper which discusses this in fuller detail.

Generating random samples pt. 3 – Metropolis-Hastings algorithm (jyyuan.wordpress.com)
Harmonic Analysis-day3 (lifeasgrad.wordpress.com)
Chernoff-Hoeffding Bounds for Markov Chains (mybiasedcoin.blogspot.com)
Reversing Markov Chains (eventuallyalmosteverywhere.wordpress.com)

Sticky Brownian Motion

Posted on April 24, 2014 by dominicyeo

This follows on pretty much directly from the previous post about reflected Brownian motion. Recall that this is a process defined on the non-negative reals which looks like Brownian motion away from 0. We consider whether RBM is the only such process, and how any alternative might be constructed as a limit of discrete-time Markov processes.

One of the alternatives is called Sticky Brownian motion. This process spends more time at 0 than reflected Brownian motion. In fact it spends some positive proportion of time at 0. My main aim here is to explain why some intuitive ideas I had about how this might arise are wrong.

The first thought was to ensure that each visit to 0 last some positive measure of time. This could be achieved by staying at 0 for an Exp(1) duration, every time the process visited it. It doesn’t seem unreasonable that this might appear as the limit of a standard SRW adjusted so that on each visit to 0 the walker waits for a time given by independent geometric distributions. These distributions are memoryless, so that is fine, but by Blumenthal’s 0-1 Law, starting from 0 a Brownian motion hits zero infinitely many times before any small time t. So in fact the process described above would be identically zero as before it gets anywhere it would have to spend some amount of time at 0 given by an infinite sum of Exp(1) RVs.

We will return later to the question of why the proposed discrete-time model will still converge to reflected BM rather than anything more exotic. First though, we should discount the possibility of any intermediate level of stickiness, where the set of times spent at 0 still has measure zero, but the local time at 0 grows faster than for standard reflected BM. We can define the local time at 0 through a limit

$L_t=\lim_{\epsilon\downarrow 0}\frac{1}{2\epsilon}\text{Leb}(\{0\le s \le t: |B_t|<\epsilon\})$

of the measure of time spent very near 0, rescaled appropriately. So if the measure of the times when the process is at 0 is zero, then the local time is determined by behaviour near zero rather than by behaviour at zero. More precisely, on the interval $[-\epsilon,\epsilon]$ , the process behaves like Brownian motion, except on a set of measure zero, so the local time process should look the same as that of BM itself. Note I don’t claim this as a formal proof, but I hope it is a helpful heuristic for why you can’t alter the local time process without altering the whole process.

At this stage, it seems sensible to define Sticky Brownian motion. For motivation, note that we are looking for a process which spends a positive measure of time at 0. So let’s track this as a process, say $C_t$ . Then the set of times when C is increasing is sparse, as it coincides with the process being 0, but we know we cannot wait around at 0 for some interval of time without losing the Markov property. So C shares properties with the local time of a reflected BM. The only difference is that the measure of times when C is increasing is positive here, but zero for the local time.

So it makes sense to construct the extra time spent at zero from the local time of a standard reflected BM. The heuristic is that we slow down the process whenever it is at 0, so that local time becomes real time. We can also control the factor by which this slowing-down happens, so define

$\sigma(s)=\rho L(s)+s,$

where L is the local time process of an underlying reflected BM, and $\rho>0$ is a constant. So $\sigma$ is a map giving a random time-change. Unsurprisingly, we now define Sticky BM as the reflected BM with respect to this time-change. To do this formally, it is easiest to define a family of stopping times $\{\tau_t\}$ , such that $\sigma(\tau_t)=t, \tau_{\sigma(s)}=s$ , then if X is the reflected BM, define $Y_t=X_{\tau_t}$ for the sticky BM.

It is worth thinking about what the generator of this process should be. In particular, why should it be different to reflected BM? The key observation is that the drift of the underlying reflected BM is essentially infinite at 0. By slowing down the process at 0, this drift becomes finite. So the martingale associated with sticky BM is precisely a time-changed version of the martingale associated with the underlying reflected BM, but this time-change is precisely what is required to give a generator. We get:

$(\mathcal{L}f)(x)=\begin{cases}\frac12f''(x)&\quad x>0\\ \rho^{-1}f'(0) &\quad x=0.\end{cases}$

Now that we have the generator, it starts to become apparent how sticky BM might appear as a limit of discrete-time walks. The process must look like mean-zero, unit-variance RW everywhere except near 0, where the limiting drift should be $\rho^{-1}$ . Note that when considering the limiting drift near zero, we are taking a joint limit in x and h. The order of this matters. As explained at the end of the previous article, we only need to worry about the limiting drift along sequences of $x,h$ such that $a_h(x)\rightarrow 0$ . If no such sequences exist, or the limiting drift along any of these is infinite, then we actually have a reflected boundary condition.

This highlights one confusing matter about convergence of reflected processes. The boundary of the discrete-time process should converge to the boundary of the reflected process, but we also have to consider where reflective behaviour happens. Can we get sticky BM with reflection only at the boundary in the discrete-time processes? The answer turns out to be no. At the start of this article, I proposed a model of SRW with geometric waiting times whenever the origin was visiting. What is the limit of this?

The trick is to consider how long the discrete process spends near 0, after rescaling. It will spend a multiple 1/p more time at 0 itself, where p is the parameter of the geometric distribution, but no more time than expected at any point $x\in(0,\epsilon)$ . But time spent in $(0,\epsilon)$ dominates time spent at 0 before this adjustment, so must also dominate it after the adjustment, so in the limit, the proportion of time spent at 0 is unchanged, and so in particular it cannot be positive.

Because of all of this, in practice it seems that most random walks we might be interested in converge (if they converge to a process at all) to a reflected SDE/diffusion etc, rather than one with sticky boundary conditions. I feel I’ve been talking a lot about Markov processes converging, so perhaps next, or at least soon, I’ll write some more technical things about exactly what conditions and methods are required to prove this.

REFERENCES

S. Varadhan – Chapter 16 from a Lecture Course at NYU can be found here.

Markovian Excursions (eventuallyalmosteverywhere.wordpress.com)
Brownian motion and stochastic calculus e-book downloads (okroxon.wordpress.com)
254A, Notes 3b: Brownian motion and Dyson Brownian motion (terrytao.wordpress.com)

Reflected Brownian Motion

Posted on April 21, 2014 by dominicyeo

A standard Brownian motion is space-homogeneous, meaning that the behaviour of $B_{T+t}-B_T$ does not depend on the value of $B_T$ . By Donsker’s Theorem, such a Brownian motion is also the limit in a process space of any homogeneous random walk with zero-drift and constant variance, after suitable rescaling.

In many applications, however, we are interested in real-valued continuous-time Markov processes that are defined not on the whole of the real line, but on the half-line $\mathbb{R}_{\ge 0}$ . So as BM is the fundamental real-valued continuous-time Markov process, we should ask how we might adjust it so that it stays non-negative. In particular, we want to clarify uniqueness, or at least be sure we have found all the sensible ways to make this adjustment, and also to consider how Donsker’s Theorem might work in this setting.

We should consider what properties we want this non-negative BM to have. Obviously, it should be non-negative, but it is also reasonable to demand that it looks exactly like BM everywhere except near 0. But since BM has a scale-invariance property, it is essentially meaningful to say ‘near 0’, so we instead demand that it looks exactly like BM everywhere except at 0. Apart from this, the only properties we want are that it is Markov and has continuous sample paths.

A starting point is so-called reflected Brownian motion, defined by $X_t:=|B_t|$ . This is very natural and very convenient for analysis, but there are some problems. Firstly, this has the property that it looks like Brownian motion everywhere except 0 only because BM is space-homogeneous but also symmetric, in the sense that $B_t\stackrel{d}{=}-B_t$ . This will be untrue for essentially any other process, so as a general method for how to keep stochastic processes positive, this will be useless. My second objection is a bit more subtle. If we consider this as an SDE, we get

$dX_t=\text{sign}(B_t)dB_t.$

This is a perfectly reasonable SDE but it is undesirable, because we have a function of B as coefficient on the RHS. Ideally, increments of X would be a function of X, and the increments of B, rather than the values of B. That is, we would expect $X_{t+\delta t}-X_t$ to depend on $X_t$ and on $(B_{t+s}-B_t, 0\le s\le \delta t)$ , but not on $B_t$ itself, as that means we have to keep track of extra information while constructing X.

So we need an alternative method. One idea might be to add some non-negative process to the BM so that the sum stays non-negative. If this process is deterministic and finite, there there is some positive probability that the sum will eventually be negative, so this won’t do. We are looking therefore so a process which depends on the BM. Obviously we could take $\max(-B_t,0)$ , but this sum would then spend macroscopic intervals of time at 0, and these intervals would have the Raleigh distribution (for Brownian excursions) rather than the exponential distribution, hence the process given by the sum would not be memoryless and Markov.

The natural alternative is to look for an increasing process $A_t$ , and then it makes sense to talk about the minimal increasing process that has the desired property. A moment’s thought suggests that $A_t=-min_{s\le t}B_t$ satisfies this. So we have the decomposition

$B_t=-A_t+S_t,$

where $S_t$ is the height of B above its running minimum. So S is an ideal alternative definition of reflecting BM. In particular, when B is away from its minimum, $dB_t=dS_t$ , so this has the property that it evolves exactly as the driving Brownian motion.

What we have done is to decompose a general continuous process into the sum of a decreasing continuous process and a non-negative process. This is known as the Skorohod problem, and was the subject of much interest, even in the deterministic case. Note that process A has the property that it is locally constant almost everywhere, and is continuous, yet non-constant. Unsurprisingly, since A only changes when the underlying BM is 0, A is continuous with respect to the local time process at 0. In fact, A is the local time process of the underlying Brownian motion, by comparison with the construction by direct reflection.

One alternative approach is to look instead at the generator of the process. Recall that the generator of a process is an operator on some space of functions, with $\mathcal{L}f$ giving the infinitissimal drift of $f(X_t)$ . In the case of Brownian motion, the generator $(\mathcal{L}f)(x)=\frac12 f''(x)$ for bounded smooth functions f. This is equivalent to saying that

$f(X_t)-f(X_0)-\int_0^t \frac12 f''(X_s)ds$ (*)

is a martingale. This must hold also for reflected Brownian motion, whenever x is greater than 0. Alternatively, if the function f is zero in a small neighbourhood of 0, it should have the same generator with respect to reflected BM. Indeed, for a general smooth bounded function f, we can still consider the expression (*) with respect to reflected BM. We know this expression behaves as a martingale except when X is zero. If f'(0)>0, and T is some hitting time of 0, then $f(X_{T+\delta T})-f(X_T)\ge 0$ , hence the expression (*) is a submartingale. So if we restrict attention to functions with f'(0)=0, the generator remains the same. Indeed, by patching together all such intervals, it can be argued that even if f'(0) is not zero,

$f(X_t)-f(X_0)-\int_0^t \frac12 f''(X_s)ds - f'(0)A_t$

is a martingale, where A is the local time process at zero.

I was aware when I started reading about this that there was another family of processes called ‘Sticky Brownian Motion’ that shared properties with Reflected BM, in that it behaves like standard BM away from zero, but is also constrained to the non-negative reals. I think this will get too long if I also talk about that here, so that can be postponed, and for now we consider reflected BM as a limit of reflected (or other) random walks, bearing in mind that there is at least one other candidate to be the limit.

Unsurprisingly, if we have a family of random walks constrained to the non-negative reals, that are zero-drift unit-variance away from 0, then if they converge as processes, the limit is Brownian away from zero, and non-negative. Note that “away from 0” means after rescaling. So the key aspect is behaviour near zero.

What is the drift of reflected BM at 0? We might suspect it is infinite because of the form of the generator, but we can calculate it directly. Given $X_0=0$ , we have:

$\frac{\mathbb{E}X_t}{t}=\frac{\mathbb{E}|B_t|}{t}=\frac{\sqrt{t}\mathbb{E}|B_1|}{t},$

so letting $t\rightarrow 0$ , we see indeed that the drift is infinite at 0.

For convergence of discrete processes, we really need the generators to converge. Typically we index the discrete-time processes by the time unit h, which tends to 0, and $b_h(x),a_h(x)$ are the rescaled drift and square-drift from x. We assume that we don’t see macroscopic jumps in the limit. For the case of simple random walk reflected at 0, it doesn’t matter exactly how we construct the joint limit in h and x, as the drift is uniform on x>0, but in general this does matter. I don’t want to discuss sticky BM right now, so it’s probably easiest to be vague and say that the discrete Markov processes converge to reflected BM so long they don’t spend more time than expected near 0 in the limit, as the title ‘sticky’ might suggest.

The two ways in which this can happen is if the volatility term $a_h(x)$ is too small, in which case the process looks almost deterministic near 0, or if the drift doesn’t increase fast enough. And indeed, this leads to two conditions. The first is straightforward, if $a_h(x)$ is bounded below, in the sense that $\liminf_{h,x\rightarrow 0} a_h(x)\ge C>0$ , then we have convergence to reflected BM. Alternatively, the only danger can arise down those subsequences where $a_h(x)\rightarrow 0$ , so if we have that $b_h(x)\rightarrow +\infty$ whenever $h,x,a_h(x)\rightarrow 0$ , then this convergence also holds.

Next time I’ll discuss what sticky BM means, what it doesn’t mean, why it isn’t easy to double the local time, and how to obtain sticky BM as a limit of discrete random walks in a similar way to the above.

REFERENCES

S. Varadhan – Chapter 16 from a Lecture Course at NYU can be found here.

Lamperti Walks

Posted on April 18, 2014 by dominicyeo

The theory of simple random walks on the integer lattice is a classical topic in probability theory. Polya proved in the 1920s that such a SRW on $\mathbb{Z}^d$ is recurrent only for d=1 or 2. The argument is essentially combinatorial. We count the number of possible paths from 0 back to itself and show that this grows fast enough that even with the probabilistic penalty of having a particular long path we will still repeatedly see this event happening. In larger dimensions there is essentially ‘more space’ at large distances, at least comparatively, so a typical walk is more likely to escape into this space.

As Kakutani (of the product martingale theorem) said, and was subsequently quoted as the dedication on every undergraduate pdf about random walks: “A drunk man will find his way home, whereas a drunk bird may get lost forever.”

But transience in some sense a long-distance property. We can fiddle with the transition rates near zero and, so long as we don’t make anything deterministic this shouldn’t affect transience properties. Obviously if we have a (space-)homogeneous nearest-neighbour random walk on the integers with non-zero drift the process will be transient: it drifts towards positive infinity if the drift is positive. But can we have a random walk with non-zero drift, but where the drift tends to zero at large distances fast enough, and the process is still recurrent? What is the correct scaling for the decay of the drift to see interesting effects?

The answers to these questions is seen in the so-called Lamperti random walks, which were a recurring theme of the meeting on Aspects of Random Walks held in Durham this week. Thanks to the organisers for putting on such an excellent meeting. I hadn’t known much about this topic before, so thought it might be worth writing a short note.

As explained above, we consider time-homogeneous random walks. It will turn out that the exact distributions of the increments is not hugely important. Most of the properties we might care about will be determined only by the first two moments, which we define as:

$\mu_1(x)=\mathbb{E}[X_{t+1}-X_t | X_t=x],$

$\mu_2=\mathbb{E}[(X_{t+1}-X_t)^2 | X_t=x].$

Note that because the drift will be asymptotically zero, the second term is asymptotically equal to the variance of the increment. It will also turn out that the correct scaling for $\mu_1$ to see a phase transition is $\mu_1(x)\sim \frac{c}{x}$ .

We begin by seeing how this works in the simplest possible example, from Harris (1952). Let’s restrict attention to a random walk on the non-negative integers, and impose the further condition that increments are +1 or -1. In the notation of a birth-and-death process from a first course on Markov chains, we can set:

$p_j:=\mathbb{P}(X_{t+1}=j+1| X_t=j), \quad q_j=1-p_j.$

We will set $p_j=\frac12 + \frac{c}{2j}$ . Then a condition for transience is that

$1+\frac{q_1}{p_1}+\frac{q_1q_2}{p_1p_2}+\ldots <\infty.$

In our special case:

$\frac{q_1\ldots q_r}{p_1\ldots p_r}\approx\frac{(r-2c)(r-1-2c)(r-2-2c)\ldots}{r!}\approx \frac{1}{r^{2c}}.$

So we can deduce that this sum converges if c>1/2, giving transience. A similar, but slightly more complicated calculation specifies the two regimes of recurrence. If -1/2<=c<=1/2 then the chain is null-recurrent, meaning that the expected time to return to any given state is infinite. If c<-1/2, then it is positive recurrent.

In general, we assume $\mu_1(x)\sim \frac{c}{x}$ and $\mu_2(x)\approx s^2$ . In the case above, obviously $s^2=1$ . The general result is that under mild assumptions on the increment distributions, for instance a $(2+\epsilon)$ -moment, if we define $r=-\frac{2c}{s^2}$ , then the RW is transient if r<-1, positive-recurrent if r>1, and null-recurrent otherwise. This is the main result of Lamperti.

To explain why we have parameterised exactly like this, it makes sense to talk about the more general proof methods, as obviously the direct Markov chain calculation won’t work in general. The motivating idea is that we can deal well with the situation where the drift is zero, so let’s transform the random walk so that the drift becomes zero. A function of a Markov chain that is more stable (in some sense) that the original MC, for analysis at least, is sometimes called a Lyapunov function. Here, the sensible thing is to consider $Y_t=X_t^\gamma$ , for some exponent $\gamma>0$ .

So long as our distributions are fairly well-behaved (eg a finite $2+\epsilon$ -moment), we can calculate the drift of Y as

$\mathbb{E}[Y_{t+1}-Y_t| X_t=x]=\frac{\gamma}{2}x^{\gamma-2}(2c+(1-\gamma)s^2) +o(x^{\gamma-2}).$

In particular, taking $\gamma=1+r$ results in a random walk that is ‘almost’ a martingale. Note that the original RW was almost a martingale, in the sense that the drift is asymptotically zero, but now it is zero to second order as well.

To draw any rigorous conclusions, we need to be careful about exactly how precise this approximation is, but we won’t worry about that now. In particular, we need to know whether we can take this approximation over the optional stopping theorem, as this allows us to say:

$\mathbb{P}(X\text{ hits }x\text{ before 0})=\mathbb{P}(Y\text{ hits }x^\gamma\text{ before 0})\sim x^{-\gamma}.$

This is particularly useful for working out the expected excursion time away from 0, which precisely leads to the condition for null-recurrence.

In his talk, Ostap Hryniv showed that this Lyapunov function analysis can be taken much further, to derive much more precise results about excursions, maxima and ergodicity. Results of Menshikov and Popov from the 90s further specify the asymptotics for the invariant distribution, if it exists, in terms of r.

One cautionary remark I should make is that earlier I implied that once we know the drift of such a random walk is zero, we have recurrence. This is true on $\mathbb{Z}$ with very mild restrictions, but is not necessarily true in higher dimensions. For example, consider the random walk on $\mathbb{R}^2$ , where conditional on $X_t$ , the increment is $X_{t+1}-X_t$ is of length 1 and perpendicular to the vector $X_t$ . The two possible directions are equally likely. The drift is therefore 0 everything, and the second moment is also well-behaved, but note that $||X_t||^2=t^2$ , just by considering Pythagoras. So in higher dimensions, we have to be a bit more careful, and put restrictions on the covariance structure of the increment distributions.

As a final comment, note that from Lamperti’s result, we can re-derive Polya’s result about SRW in higher dimensions. If we have $X_t$ an SRW on $\mathbb{Z}^d$ , then consider $Y_t=||X_t||$ . By considering a couple of examples in two-dimensions, it is clear that this is not Markov. But the methods we considered above for the Lamperti walks were really martingale methods rather than Markov chain methods. And indeed this process Y has asymptotically zero drift with the right scaling. Here,

$c=\frac{1}{2}(1-\frac{1}{d}),\quad s^2=\frac{1}{d},$

and so r=d-1, leading to exactly the result we know to be true, that the SRW is transient precisely in three dimensions and higher.

REFERENCES

Harris – First Passage and Recurrence Distributions (1952)

The slides from Ostap Hryniv’s talk, on which this was based, can be found here.

Probability that a uniformly random sequence is already sorted

Local to Global in Point Set Combinatorics

Posted on April 8, 2014 by dominicyeo

I’ve spent the last few days in Cambridge, helping to run the first camp to select the UK team for the International Mathematical Olympiad. We’ve had two 4.5 hour exams, replicating what six of the students will face in the final competition, which this year is held in Cape Town, and also lots of problem-solving sessions, introducing the students to different areas of mathematics relevant to the olympiad. For various reasons, I ended up choosing to discuss combinatorial problems involving points and lines with both sets of students. This is an area which is very popular with problem setters (witness for example Question Two at last year’s IMO in Colombia) but which, because of a lack of any real ‘theory’, is often not treated at these camps.

Anyway, 90 minutes is not a lot of time to solve a large number of problems, so I was interested in the notion of using the time to try to teach the students how to have the right sort of ideas to deal with these. Obviously if there were an easy way to teach people to have good ideas, mathematically or otherwise, life would be a whole lot more straightforward. So, given we can’t do that, the next best thing is to have lots of ideas of any kind, and to get practice at choosing which ones seem best reliably and quickly.

We have two problems to consider that are very similar. I’ll state the first one:

Q1: a strip of width w is the set of all points which lie on or between two parallel lines that are a distance w apart. Let S be a set of n>2 points such that any three can be covered by a strip of width 1. Prove that S can be covered by a strip of width 2.

This sort of thing comes up all the time in olympiad combinatorics and beyond. We have some condition that is local, in the sense that it only affects triples of points, but we want to show that the same condition works for all the points together.

The first step is to be clear about what we are looking for. We are looking for a strip, but we have also been given lots of strips, namely the strips covering each triple of points. So the natural thing to do is to construct the strip we are looking for from the strips we have.

The next observation offered was that we might want to combine two strips of width 1 to get a strip of width 2. This is certainly an idea, but a quick check makes you realise it isn’t a very promising idea. Why? Well unless two strips are parallel, it takes an infinite-width strip to cover their union. Ok, so although that wasn’t useful in itself, it now encourages us to focus on constructing the strip we are looking for from a single one of the strips we are given.

It’s not yet clear exactly how we are going to construct the desired strip from the strip we are given. We should be alive to the possibility that 2 might be a weak bound, and actually $1+\epsilon$ would suffice. Or we might have tried some examples. One thing we might have is a strong example of where a strip of width 2 is necessary, and where we can’t add any points which don’t fit in this strip without breaking the first condition.

Maybe this doesn’t help so much. So perhaps we should think not about how to construct the big strip from the little strip, but instead, how to choose the little strip from the large set we have at our disposal. It is at this stage that I feel it is really important to use all the power at our disposal. In many combinatorics problems we have dealing with a discrete set like [n], where the roles of the labels are interchangeable, and so it can sometimes to be difficult to pick out a particular element. For example, it doesn’t matter for induction whether we take [n-1] then add {n} or start with [2,n] and add {1}. When the construction is geometric, we should use the geometry to make life easier in the construction. For example, if we want to induct, we could add the left-most point in some co-ordinate axes, or some point on the convex hull.

In our particular example, we are still trying to decide which strip to work with. If we consider the triple of points which is in some sense hardest to cover with a strip, this might be a sensible thing to work with. Now we need to decide what it means to be hard to cover a triangle with a strip. We might observe that if the triangle has an altitude of height h, we can certainly cover the triangle with a strip of width h. Further playing around convinces us that the minimum width of a strip required to cover a triangle is given by the shortest altitude. So our intuition that it’s hard to cover a triangle if it is large has been made a bit more concrete. We might also observe that we don’t have much flexibility for a strip if one of the lengths of the triangle is very long. For fixed strip width, as the length of the triangle gets large, the range of angles the strip could sit at (if it is possible at all) gets small.

This motivates considering the least convenient triangle as the one with maximal width and maximal height of altitude. We might try both orders, but certainly considering the triangle ABC such that d(A,B) is maximal among the set of points, and d(AB, C) is maximal once A,B are determined, is useful. The fact that we can cover ABC with one of the strips tells us from the above that $d(AB,C)\le 1$ , and suddenly we are almost done, as we now have a strong clue for what the strip we are looking for should be. If we take all the points x such that $d(AB,x)\le 1$ , we have a strip of width 2, and by construction of C we must have included all points of the set.

The key ingredient here was using what we already had, realising it probably wasn’t going to be helpful to attempt to use *all* of what we already had, then coming up with a sensible way to choose which pre-existing strip to develop.

With this in mind, the next question seems a lot easier:

Q2: Given n points in the plane such that any three lie within some circle of radius 1, prove that *all* the points lie within some circle of radius 1.

Circulating round the room, there were lots of pairs of compasses out and various diagrams of interlocking circles snaking around pieces of A4. But after discussing Q1, suddenly all seems a lot more clear. Similarly to the previous case, we have a lot of circles, and we are looking for a circle. Conveniently, the circle we are looking for is even the same size as the circles we have.

So we need to choose which circle we already have is likely to be most useful. Suppose we have three points very close together. Then we know very little about the circle containing those points – it could be almost anything. By contrast, if we have three points which form a triangle with circumradius equal to 1, then there is in fact no choice about what the covering circle could be. So this is most likely to be useful. So we conjecture that we should consider the triangle with largest circumradius. In fact, this is a problem: if the triangle is obtuse, then the circumradius might actually be greater than 1, but the points still satisfy the condition. In that case, it makes sense to consider the circle with diameter given by the longest side of the triangle.

With that clarification, we now have a collection of circles with radii all $\le 1$ , which cover the triangles. We conjecture that the circle with largest radius might cover all the points. This turns out to be correct, and not too difficult to prove. Both these problems are in some sense extremal. Only here the tricky part was deciding what to extremise. But all were agreed that with some meta-thinking, we were able to decide fairly accurately which ideas were leading towards progress, and following them led fairly rapidly to this extremal property.

Mixing of the Noisy Voter Model

Posted on March 23, 2014 by dominicyeo

I’ve been spending a bit of time recently thinking about some interacting particle systems and mixing time problems, so it was interesting to find a short paper due to Harishchandra Ramadas on arXiv a couple of weeks ago combining these two areas. I enjoyed reading it, so I thought I would summarise its arguments here, in preparation for perhaps talking about it at this week’s junior probability seminar.

The voter model is a graph-based interacting particle system. The underlying motivation is to describe the spread of opinions through a population. Each vertex represents an agent who has an opinion which might vary over time. At exponential rates, an agent will replace their opinion with the opinion of one of their neighbours, chosen uniformly at random. Note that, obviously, this might not actually involve a change of opinion.

One problem with this model is that we might be interested in convergence to equilibrium, but this might be hard to specify, as the state where everything is 0 and the state where everything is 1 are absorbing, and similarly any binomial distribution on the set of vertices is invariant. It is also not entirely realistic to assume that future updates to the system are entirely based on the current configuration. In any actual application, we should allow for some, possibly very small, influence from the outside, whether due to random errors in the update rule, or some external factors which haven’t been explicitly accounted for.

So we consider instead the noisy voter model where, in addition to the dynamics described above, each vertex changes state at rate $\delta$ , independently of the configuration or the history. [Note that we are now considering the set of opinions {0,1} ]

We are interested in the mixing time of this system. It is not immediately clear that this should have an invariant distribution, but it seems reasonable since the noise, which converts 0 to 1 and 1 to 0 at equal rates if the proportion of 0s and 1s is equal, moves the system towards something like Bin(V(G),1/2), although the correlations need to take account of the local geometry.

Anyway, we assume $\delta>0$ is fixed, and n tends to infinity. Ramadas proves that the dynamics have the property of optimal temporal mixing, which states that given copies of the process X, Y starting from different initial conditions

$||X_t-Y_t||\le C n e^{-ct},$

for suitable choice of constants. Note that it doesn’t matter whether X and Y are coupled or independent. Total variation distance is a function of distribution, not of probability space. I don’t have a great intuition for exactly what this means, except to make the obvious comment that the distance from uniformity grows at most uniformly as the size of the system grows. Note that this immediately implies that the mixing time is O(log n). The rest of the paper proves that the noisy voter model has this property.

The method is to consider the relation:

$||X-Y||\le \mathbb{P}(X\ne Y),$

which says that the total variation distance gives the minimum of the RHS probability over couplings of X and Y. Anyway, as discussed in the previous post on duality in such interacting particle systems, the key step is to consider the process as a deterministic walk through a random environment. Ie, we apply the same moves to configurations started from different initial conditions.

The first step is to break down the noise mechanism to make it easier to couple everything. Instead of saying that a vertex changes state at rate $\delta$ , we instead have two noise processes. One sets a vertex’s state to 0 at rate $\delta$ , and another sets a vertex’s state to 1 at rate $\delta$ independently. Since at any configuration, for each vertex, precisely one of these processes has no effect.

So for the coupling, we run the same noise processes, and we choose the neighbours from which to update our opinions the same in both processes. The only difference is the initial configuration. With these two processes, X and Y defined on the same probability space, Ramadas defines $Z_t(v)=\mathbf{1}\{X_t(v)\neq Y_t(v)\}$ , which counts how many vertices are different between the two processes. Almost surely, eventually all the Zs will be 0, and this state is absorbing because of how we’ve defined this coupling. We are interested to find out how long it takes for this event, which we could think of like a strong stationary time, to happen.

To control this, we track the evolution in the probabilities of $Z_t(v)=0$ jointly. We note that after a noise event at vertex v, we definitely have $Z_t(v)=0$ . After an opinion-change event, whether $Z_t(v)$ changes depends on whether $Z_t(u)$ is 1, where u is the uniformly-chosen vertex from which v inherits its opinion in both processes. We obtain the following ODEs:

$\frac{d}{dt}\mathbb{P}(Z_t(v)=1)=\frac{1}{d(v)}\left(\sum_{u:u\sim v}\mathbb{P}(Z_t(v)=0,Z_t(u)=1)\right)$

$-\frac{1}{d(v)}\left(\sum_{u:u\sim v}\mathbb{P}(Z_t(v)=1,Z_t(u)=0)\right)-2\delta \mathbb{P}(Z_t(v)=1).$

We interpret this sum as a negative bit plus two terms which we would hope would cancel each other out. It seems impractical to attempt to formalise this heuristic. However, we can instead choose to consider at time t to consider the ODE for whichever vertex v maximises $Z_t(v)$ . It’s then clear that the RHS is less than $-2\delta \mathbb{P}(Z_t(v)=1)$ .

Right-differentiabililty of $Z_t(v_{max})$ follows from differentiability of $Z_t(v)$ . The number of vertices is finite, and almost surely the vertex attaining this maximum can’t change too often. In particular, from the ODE, we must have $Z_t(v_{max})$ bounded above by $Ce^{-2\delta t}$ . From which it follows that

$\mathbb{P}(X\ne Y)\le nCe^{-2\delta t },$

and we’ve shown this optimal temporal mixing property, and hence the mixing time has order log n.

What I’ve discussed so far is taken entirely from Section 3 of Ramadas’s paper. My first thought was whether this bound is actually correct in scale, as we have used only the size of the graph rather than any finer information about its graph structure. The author discusses that it is not possible to use a general result due to Hayes and Sinclair on this topic in this setting, as it is not reversible.

I certainly can’t answer this question, but I will make a comment on the extreme examples. If the graph is the empty graph on n vertices, then the voter model dynamics are null, so we essentially have just a coupon collector problem going on. It doesn’t lose any generality to start from entirely 0s, then at time t, the distribution of the system will be Bin(n,p), with p a function of time, but not n. Now note that

$||\mathrm{Bin}(n,p)-\mathrm{Bin}(n,q)||_{TV}=\Theta(\sqrt{n}|p-q|),$

so in this example, we need to choose a time so that $p=\frac12+\Theta(\sqrt{n})$ . This is really a statement about the central limit theorem. It’s clear that Bin(n,1/2) should genuinely be the equilibrium distribution here. The evolution of p(t) is exponential, and hence we get the O(log n) scaling emerging. Note that the presence of $\sqrt{n}$ rather than any other power of n makes no difference, though it would be relevant if we actually wanted to know the constant.

At the other end of the connectivity spectrum, we consider the dynamics on the complete graph. Here, let’s count the number of 1s, and call it W(t) say. Then W goes from k to k+1 at rate $\frac{(n-k)\delta}{2}+\frac{k(n-k)}{n}$ , and from k to k-1 at rate $\frac{k\delta}{2}+\frac{k(n-k)}{n}$ . In particular, note in expectation, the contributions from the voter model dynamics cancel. This fits with our intuition that it is solely the noise driving the system towards equilibrium.

If we take an ODE approximation, we get:

$\frac{dX}{dt}=\frac{n}{2}-\delta x,$

from which it follows that

$t(x)=\alpha-\frac{1}{\delta}\log(C-\delta x),$

which gives precisely the correct scaling, when we take x=n/2.

I would be very interested to think about or hear about arguments for a general lower bound across all graphs.

REFERENCES

Hayes, Sinclair (2007) – A general lower bound for mixing of single-site dynamics on graphs (arXiv)

Ramadas (2014) – Mixing of the Noisy Voter Model (arXiv link at top of article)

Eventually Almost Everywhere

A blog about probability and olympiads by Dominic Yeo

Hitting Probabilities for Markov Chains

Avoiding Mistakes in Probability Exams

Popoviciu’s Inequality

A Month in Vancouver

Coupling from the Past

Sticky Brownian Motion

Reflected Brownian Motion

Lamperti Walks

Local to Global in Point Set Combinatorics

Mixing of the Noisy Voter Model

Related articles

Related articles

Related articles

Related articles

Related articles

Related articles

Related articles

Related articles