5 ms·

Markov Chains vs. Simulation: Flipping a Million Little Coins

75 points by JHonaker 12y ago

neaanopri 12y ago
But by not sampling at each step, we only get the expectation value at any given time?
- JHonaker 12y ago
  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.
  - murbard2 12y ago
    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.
    - JHonaker 12y ago
      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.
      - murbard2 12y ago
        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.
    - omalleyt 12y ago
      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.)
      - murbard2 12y ago
        It's a markov chain if your state is the vector representing the size of each cell
- obastani 12y ago
  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.
termain 12y ago
Baseball is more an illustrative example than an analogy. Consider http://www.pankin.com/markov/intro.htm http://www.pankin.com/markov/intro.htm
adamtj 12y ago
"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.
- Florin_Andrei 12y ago
  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).
  - gjm11 12y ago
    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.
physPop 12y ago
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
- indeed30 12y ago
  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.
- joe_the_user 12y ago
  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/Markov_chain http://en.wikipedia.org/wiki/Continuous-time_Markov_chain http://en.wikipedia.org/wiki/Continuous-time_Markov_chain
- mturmon 12y ago
  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).
Balgair 12y ago
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.
- seslattery 12y ago
  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
- ende 12y ago
  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.

75 points by JHonaker 12y ago

neaanopri 12y ago

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

JHonaker 12y ago
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.
- murbard2 12y ago
  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.
  - JHonaker 12y ago
    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.
    - murbard2 12y ago
      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.
  - omalleyt 12y ago
    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.)
    - murbard2 12y ago
      It's a markov chain if your state is the vector representing the size of each cell
obastani 12y ago
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.

termain 12y ago

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

adamtj 12y ago

"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.

Florin_Andrei 12y ago
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).
- gjm11 12y ago
  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.

physPop 12y ago

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

indeed30 12y ago
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.
joe_the_user 12y ago
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/Markov_chain http://en.wikipedia.org/wiki/Continuous-time_Markov_chain http://en.wikipedia.org/wiki/Continuous-time_Markov_chain
mturmon 12y ago
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).

Balgair 12y ago

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.

seslattery 12y ago
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
ende 12y ago
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.