Markov Chains vs. Simulation: Flipping a Million Little Coins | Hacker News Reader

Markov Chains vs. Simulation: Flipping a Million Little Coins

thexbar.me

75 pointsJHonaker23 comments

Hi. As other posters have mentioned there is a closed form solution that comes straight from an easy differntial equation.

I'd like to make a point to help the author out, since this post is in the spirit of simple explanations: a good way of thinking of Markov models is that they are fundamentally a discreet form of differential equations, and are useful when we have finitely few states that are generally:

a) easily enumerable b) countably small enough that it doesn't make sense to just replace it with a continuous model

This is similar to the way in which martingales are discreet representations of random walks. Sadly, our education system does a bad job of discreet math, and its typically left to stats courses -- so everything has a different name and few people draw the connections, IMO!

Cheers

I'm struggling with this comment, I don't necessarily disagree with the spirit of it, but I feel that you're doing Markov chains something of a disservice by claiming they're mainly useful for "finitely few states".

You're completely ignoring immensely important topics like MCMC which (can) operate on uncountable spaces. The vast, vast majority of MCMC research and usage is with continuous random variables.

Well, Markov chains as used by programmers are generally discreet-time. Continuous-time Markov chains are continuous. Indeed, in the original model he begins with what seems to be a continuous-time Markov chain. You are correct you get to solve a linear differential equation with constant coefficients to get the distribution at any one time. Another interesting way to consider this stuff is in terms of the "logarithms" and "exponentiation" of Markov matrices.

http://en.wikipedia.org/wiki/Markov_chain http://en.wikipedia.org/wiki/Continuous-time_Markov_chain

This comment is problematic. You're confusing discrete-state-space models with discrete-time models. You can have any combination of continuous/discrete Markov process:

Discrete states, discrete time: text, etc.

Discrete states, continuous time: Poisson counting process

Continuous states, discrete time: gambling earnings

Continuous states, continuous time: Brownian motion

Of course, the term "Markov chain" is restricted to discrete time. But you have said "Markov model", which includes all processes with the Markov property (i.e., they forget old state). The model in the OP was not originally a chain, because the radioactive decay process is continuous-time, but he made it into one.

When you say that Markov chains are discretized differential equations, it tempts the reader to think "discretized in time" -- but that's called a difference equation. I think you mean it's discretized in state. I.e., there is a flow of mass from one state to another that is smooth in a diff. eq., and stochastic in the Markov chain. But that's just guesswork.

Also, I would not say a martingale is a discrete random walk. Your statement here introduces the same confusion. Like Markov processes, martingales can be continuous or discrete in both state-space and time (e.g., Brownian motion is a martingale).

I think of a martingale as a mathematical model for a fair game (expected future earnings is zero given available information).

Fair concerns, thank you for clarifying the notation, you are correct in my misuse of model/chain. Was typing fast :)

After some further reading , I also agree re: martingales. To be honest I have only ever encountered them in the common discreet-time situation.

"Simulation" in this case is just numeric integration the hard way, even though there's a closed form solution.

This is a calculus problem, not a programming problem. After you solve the real problem, the programming is trivial. If you don't solve the real problem, you can't really be sure your simulation gives you the correct answer.

Exponential decay is a lot like compound interest. Simulating this problem is like computing interest by compounding yearly instead of continuously. It's not hard for those methods to differ by a factor of two or ten or even more. You can't really know if your simulation will be affected without doing the math anyway.

Well, when computing power is plentiful, sometimes it's easier to just simulate the heck out of it. (shrug)

Back in the day, I was curious about the ratio between the volume of an N-dimensional sphere, divided by the volume of the N-dimensional cube that contains it, as a function of N.

I'm sure there might be an exact solution somewhere (it's basically the volume of the N-sphere of radius 0.5) but I was too lazy to look it up. Instead, I just did a trivial one-page Monte Carlo shoot-em-up with a lot of points. It converged quickly to a good estimate.

For the curious, as N goes up, the volume of the N-sphere is a smaller and smaller fraction of the volume of the N-cube containing it. It kinda blew my mind when I saw the result, but then it started to make sense (in higher dimensions there's more ways for space to be "wasted" in the corners).

Yup, there's an exact solution. It turns out that the volume of the n-dimensional unit ball is pi^(n/2) / (n/2)!. (The factorial function generalizes to numbers that aren't integers. For half-integers you get a rational number times sqrt(pi). E.g., if n=3 then we need (3/2)! which happens to equal 3/4.sqrt(pi), so the volume is 4/3.pi.)

So for large n we can use Stirling's approximation. To make the formulas a little shorter I'll write m=n/2. Then the volume of the unit n-ball is

V(n) = pi^m/m! ~ pi^m (e/m)^m 1/sqrt(2.pi.m)

whereas the volume of the smallest n-cube containing that ball is 2^n, so the ratio is

R(n) ~ (pi.e/2m)^m / sqrt(2.pi.m)

where, remember, m=n/2. Actually, at this point maybe it's nicer to put things back in terms of n:

R(n) ~ sqrt[(pi.e/n)^n/pi.n]

When n is large, you should see this as

R(n) ~ sqrt[(..../n)^n.....]

where the "...."s are the less important bits. So for large n this is going to get small.

Concrete not-so-large-n example: let's take n=10. Then R(10) ~= sqrt((pi.e/10)^10/10pi) which happens to be about 0.081. That's pretty accurate -- the exact answer is pi^5/3840, which is about 0.0797. So if you fit a 10-ball snugly inside a 10-cube it'll take up about 1/12 of the space.

What if we crank the dimension up? With n=100 the ratio is about 2x10^-55. (A naive Monte Carlo simulation will be unable to distinguish this from zero.) So if you do the ball-inside-cube thing in 100 dimensions, the ball occupies less than one billion-billion-billion-billion-billion-billionth of the volume of the cube.

High-dimensional cubes are really, really spiky.

Totally true.

It's a useful heuristic approximation that

  n! ~ (n/e)^n

That is, multiplying the "n" factors from 1...n is close to multiplying (n/e) upon itself "n" times.

The actual Stirling approximation is just a sqrt(2 pi n) bigger than the above, and usually it's fine to throw that away (as above).

The point of the simulation portion was really to just explain that even though you know almost nothing about the theoretical behavior (from a probabilistic point of view) you can just apply what happens at each time point over and over to get an idea of how it behaves.

For example, using basically the same process of building from a base case to simulating a system, a few months ago, I was interested in determining roughly how many packs of cards from a certain game I would need to open to be reasonably sure I would have all card in the set. I took some probabilities that I'd gotten from other people opening hundreds of packs, and used that to simulate opening one pack. I then performed about a thousand runs of opening packs until I got every card and recorded the number of packs it took each time.

In the end I got a really nice normal distribution of the estimate of the mean number of packs it would take to get every card (hooray for the central limit theorem!).

Discrete time and continuous time exponential models can be converted from and to each other, so that they will agree anytime the discrete time model is defined.

A bit off topic: I'm currently struggling through a bio course and a lot of modern molecular biology is almost nothing but pathways for one protein to interact with another. I used to do reliability analysis and thought that the chains of proteins looked incredibly similar to Markov chains. Does anyone know of a model in biology that used markov chains explicitly? Thank you in advance.

Yes, MSMBuilder is software package that uses them explicitly in protein dynamics simulations. The rationale behind it is it allows you to take trajectories from molecular dynamics simulations, and construct a markov state model out of the conformational state space of the protein. This is extremely useful in MM simulations, because it allows for the conformational space to be sampled by many smaller trajectories, instead of one longer trajectory where it is highly unlikely to sample the full state space anyways(gets stuck in local minima). This allows for molecular dynamics simulations to be done as High Thruoughput computing instead of HPC. Its the rationale behind Folding@home.

Edited for clarification

Yes, they're all over the place and lend themselves well to nucleotide and protein sequence data. Hidden Markov Models are used in the Glimmer gene finding algorithm for example.

Baseball is more an illustrative example than an analogy. Consider http://www.pankin.com/markov/intro.htm

But by not sampling at each step, we only get the expectation value at any given time?

edit: I completely thought I was answering a different topic. Please disregard the quoted, but I will leave it for completeness.

> I'm not quite sure what you mean by this question.

> What do you mean by "not sampling at each step"? The idea isn't to characterize the entirety of the sampling function along its support, but to draw a sequence of random variables from it.

Yes but that is not what you do. You calculate the expected proportion of each atom at each step, but:

1) this does not tell you what the distribution is. Decay is a random process... what is the probability that you have more than 90% of a left at step 1000?

2) this only works because you're taking a limit over an infinite number of atoms, which takes the variance to 0. If you have a discrete system, things become hairier.

Or take the following system for instance:

There are 100 cells, each of which has the same size. At every time step, a little bit of food falls on the cells and is eaten by a single cell. A cell gets the food with a probability proportional to its weight, it then grows by one unit.

If you use expectations, everyone grows at the same rate. If you simulate, you'll see that the results follow a Pareto distribution.

Whoops, I thought I was answering a question about my rejection sampling article. My mistake. I'll address this though.

1) No it doesn't. I don't make any distributional claims for these calculations. You're right. I haven't addressed questions like P(>90%|t=1000).

2) Actually no, it doesn't take the variance to 0. For any given t, the Markov approximation is equal to the expected value of the given atom at time t. There is a variance for every atom count at every time point. I didn't address calculating variance at all. We're speaking in expected values here, the interpretation of 54.8 atoms at time t is the same as the average american household having 2.5 children.

Something I addressed in a comment when I posted this on reddit was that the Markov chain approach is the theoretical mean. If you were to actually observe a decay process like this in the real world, it would essentially be one possible realization of the simulated example.

As far as the Pareto distribution is concerned, I didn't know that actually. That's really interesting.

Ah got it, you simulate individually per atom, I thought you were describing the entire population. Indeed, you can always use that vector as the parameter to a multinomial distribution and recover the real distribution.

Regarding the Pareto distribution, random growth model are probably the reason why power laws are everywhere.

The cell example is very cool. It's not a Markov chain though. The probability that a cell grows from size 2 to size 3 (in probabilistic notation, P(2->3|2) ), for instance, would be 2/(sum of all cell sizes), and the sum of all cell sizes is growing by 1 each period (assuming the unit of food given per period is 1, etc. etc.)

It's a markov chain if your state is the vector representing the size of each cell

Yes, but you can get other quantities with this approach as well (though not in nearly as nice a form). To do this, consider the random matrix

t = Bernoulli(T) (i.e., T = E[t])

The distribution is

m_n = M_1 t_1 ... t_n,

which has expected value

M_n = E[m_n] = E[M_1 t_1 ... t_n] = M_1 T^n

Here, each of the t_i is an i.i.d. sample of t. You can similarly compute the variance:

Var[m_n] = Var[M_1 t_1 ... t_n]

I don't know a closed form for this off the top of my head, but at this point simulations are probably the best way to go.

Reply on news.ycombinator.com