MCMC with Gibbs Sampling and MH for straightforward stuff is
straightforward, but the subtitles of underflow, use log space or not
etc are something you guys know quite a bit more about than I do.
FWIW, a few months ago I was doing some custom Racket coded for
straightforward Gibbs (mainly) and MH in one case for customer related
data analysis, but the water gets deep quickly beyond the
straightforward, and one quickly questions the validity of their custom
sampler against considerations of a proven MCMC library. I did a brief
amount research on the current state of things with regards to WinBugs,
OpenBugs, Jags, ... After a quick overview, I sort of narrowed things
down to a relatively new arrival, STAN. https://code.google.com/p/stan/
It's a new rewrite, based on a new sampling algorithm variation
where existent libraries are getting long in the tooth. It generates
the sampling code rather then interprets BUGS DSL code and claims to
take some pains to support FFI binding (R of course) and embedding.
Just wanted to mention STAN if you haven't run across it yet, as
opposed to well know standby MCMC libs.
One day ... an FFI Racket <-> STAN would be very cool.
On Jan 4, 2013 4:47 PM, "Neil Toronto" <neil.toro...@gmail.com
<mailto:neil.toro...@gmail.com>> wrote:
I get excited about applying statistics to programming. Here's
something exciting I found today.
I'm working on a Typed Racket implementation of Chris Okasaki's
purely functional random-access lists, which are O(1) for `cons',
`first' and `rest', and basically O(log(n)) for random access. I
wanted solid randomized tests, and I figured the best way would be
bisimulation: do exactly the same, random thing to a (Listof
Integer) and an (RAList Integer), then ensure the lists are the
same. Iterate N times, each time using the altered lists for the
next bisimulation step.
There's a problem with this, though: the test runtimes are all over
the place for any fixed N, and are especially wild with large N.
Sometimes it first does a bunch of `cons' operations in a row.
Because each bisimulation step is O(n) in the list length (to check
equality), this makes the remaining steps very slow. Sometimes it
follows up with a bunch of `rest' operations, which makes the
remaining steps very fast.
More precisely, random bisimulation does a "simple random walk" over
test list lengths, where `cons' is a step from n to n+1 and `rest'
is a step from n to n-1. Unfortunately, the n's stepped on in a
simple random walk have no fixed probability distribution, and thus
no fixed average or variance. That means each bisimulation step has
no fixed average runtime or fixed variation in runtime. Wildness
ensues.
One possible solution is to generate fresh test lists from a fixed
distribution, for each bisimulation step. This is a bad solution,
though, because I want to see whether RAList instances behave
correctly *over time* as they're operated on.
The right solution is to use a Metropolis-Hastings random walk
instead of a simple random walk, to sample list lengths from a fixed
distribution. Then, the list lengths will have a fixed average,
meaning that each bisimulation step will have a fixed average
runtime, and my overall average test runtime will be O(N) instead of
wild and crazy.
First, I choose a distribution. The geometric family is good because
list lengths are nonnegative. Geometric with p = 0.05 means list
lengths should be (1-p)/p = 19 on average, but are empty with
probability 0.05 and very large with very low probability.
So I add these two lines:
(require math/distributions)
(define length-dist (geometric-dist 0.05))
In a Metropolis-Hastings random walk, you step from n to n' with
probability min(1,P(n')/P(n)); this is called "accepting" the step.
Otherwise, you "reject" the step by staying in the same spot. Adding
this new test takes only two defines and a `cond':
(define r (random 5)) ; Randomly choose an operation
(cond
[(= r 0)
;; New! Only cons with probability P(n+1)/P(n)
(define n (length lst))
(define accept (/ (pdf length-dist (+ n 1))
(pdf length-dist n)))
(cond [((random) . < . accept)
.... Accept: cons and ensure continued equality ....]
[else
.... Reject: return the same lists ....])]
[(= r 1)
.... Accept: `rest' each list and ensure continued equality
....]
.... Other random operations ....)
I left off the acceptance test for `rest' because P(n-1) > P(n) ==>
P(n-1)/P(n) > 1, so the test would always pass. (This is because the
geometric distribution's pdf is monotone. If I had chosen a Poisson
distribution, which has a bump in the middle, I'd have to do the
acceptance test for both `cons' and `rest'.)
I instrument the test loop to record list lengths, verify that their
mean is about 19, and do some density plots vs. the geometric
distribution. It all looks good: the randomized bisimulation test
does indeed operate on lists whose length is distributed
Geometric(0.05), so the overall average test runtime is O(N), and it
still tests them as they're operated on over time.
Neil ⊥
____________________
Racket Users list:
http://lists.racket-lang.org/__users
<http://lists.racket-lang.org/users>
