Probabilistic programming is the natural extension of this line of development. It's an active area of research. Two such languages implemented in Scala are Figaro and Factorie. Not sure about Haskell.
(Also, weird reading a comment I made on that blog 8 years ago!)
That lead quote is highly ironic considering that MSR employs arguably the most talented collection of mathematical probabilists anywhere in the world (including academia).
Imagine a random number generator that picks 128 bits from a uniform source of 0's and 1's. Then, if I understand this correctly, if this RNG is implemented using PFP, the result will consist of 2^128 different bit-strings, each with the same probability (?)
How is this represented internally? What is the time and space complexity? Is there even a way to practically do this?
Probability monad is a really novel way to experiment with stochastic processes. I saw one in scala few weeks back. To get insight into 128-bit strings, you can transform/group the sample space such as finding the P(1st bit is 1) or get a histogram of P(# 1s in string). Otherwise, each bitstring will be almost unique. The code is really concise and execution is less than a second.
val d: Podel[Vector[Int]] = Bernoulli(p = 0.5).repeat(128).map(_.toVector)
println(d.probability(bits => bits(0) == 1))
println(d.map(bits => bits.count(_ == 1)).hist(sampleSize = 10000, optLabel = Some("Number of 1s in 128 bits")))
I can imagine a simple scheme by which this works: simply keep a "distribution" with every single variable. However, I guess this stops working correctly when there is correlation. For example, imagine that somewhere in the program, I generate a variable X with 50% probability of being 1 and 50% of being 0. Then I generate Y = -X. And finally, I compute Z = X + Y. Now, if the program does not keep track of the fact that Y and X are correlated, then Z will falsely contain nonzero values in its distribution.
Perhaps this example is what you describe? flatMap (used by the for comprehension) enables monads to combine in powerful ways. The first part shows uncorrelated X,Y such that Y ~ -X, then the second part shows correlated Y = -X such that every y = -x
{
val X: Podel[Int] = Bernoulli(p = 0.5)
val Y: Podel[Int] = X.map(-1 * _)
val Z: Podel[Int] = X + Y
val d: Podel[Int] = for {
x <- X
y <- Y
z <- Z
} yield {
z
}
d.hist(sampleSize = 10000, optLabel = Some("Uncorrelated X, Y"))
}
{
val X: Bernoulli = Bernoulli(p = 0.5)
val d: Podel[Int] = for {
x <- X
y = - 1 * x
z = x + y
} yield {
z
}
d.hist(sampleSize = 10000, optLabel = Some("y = -x"))
}
/*
So far, the weakness I find with probability monad, and simulation, is lack of precision. This makes them useless when the probability you want to find is precise such as 5.324e-11. Only formulas can give you the answer.
The naive representation that you're pointing at would indeed be very slow. Better representations examine the language of the problem as you describe it and then build efficient factor graph representations of your problem and do message passing simulation (see Factorie).
