
Markov Chains vs. Simulation: Flipping a Million Little Coins - JHonaker
http://thexbar.me/2014/10/24/raddecay/
======
physPop
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

~~~
mturmon
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).

~~~
physPop
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.

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

~~~
mturmon
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).

------
Balgair
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
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

------
termain
Baseball is more an illustrative example than an analogy. Consider
[http://www.pankin.com/markov/intro.htm](http://www.pankin.com/markov/intro.htm)

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

~~~
JHonaker
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
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
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
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.

