Probabilistic Programming

Rapid, easy model development

BAli-Phy contains a simple language for expressing probabilistic models as programs. Inference on parameters can then be performed automatically. Such languages are called probabilistic programming languages (PPL). Other well-known PPLs include BUGS, RevBayes, and Stan. The goal of the language is to allow researchers to spend their time designing models instead of designing new inference programs. The inference should take care of itself after the model is specified.

Language properties

The modeling language is a functional language, and uses Haskell syntax. Features currently implemented include:

  1. Functions work, and can be used to define random variables.
  2. Modules work, and allow code to be factored in clean manner.
  3. Packages work, and allow researchers to distribute their work separately from the BAli-Phy architecture.
  4. Optimization works, and speeds up the code written by the user by using techniques such as inlining.
  5. Composite Objects work, and can be used to define random data structures.
  6. Random control flow works, allowing if-then-else and loops that depend on random variables.

Features that are expected to be completed during 2019 include:

  • Random Trees. Random processes on trees that are not known in advance. (partially implemented)
  • Rooted Trees. Rooted trees implemented as a data structure within the language. (partially implemented)
  • JSON logging. This enables logging inferred parameters when their dimension and number is not fixed. (partially implemented)
  • Random numbers of random variables. Random variables can be conditional created, without the need for reversible-jump methods.
  • Lazy random variables. Infinite lists of random variables can be created. Random variables are only instantiated if they are accessed.
  • Type checking. Type checking will enable polymorphism and give useful error messages for program errors.

Examples

Simple linear regression

Here is a short program that performs linear regression. The program samples variables from their prior distribution using the sample function, and then conditions on the data using the observe function.

import Data.ReadFile
import Distributions
import System.Environment

xs = read_file_as_double "xs"

ys = read_file_as_double "ys"

main = do

  b <- sample $ normal 0.0 1.0

  a <- sample $ normal 0.0 1.0

  s <- sample $ exponential 1.0

  let f x = b*x + a

  sequence_ [observe y (normal mu_y s) | (x,y) <- zip xs ys, let mu_y = f x]

  return $ log_all [b %% "b", a %% "a", s %% "s"]

Interpretation:

  • let f x = b*x +a defines the prediction function f
  • b <- sample $ normal 0.0 1.0 places a prior on the slope of the line.
  • observe y (normal mu_y s) says that the data y comes from a normal distribution with mean mu_y and standard deviation s.

You can find this file here and run it as bali-phy -m LinearRegression.hs --iter=1000.

(The method of reading from files "xs" and "ys" here is kind of a hack. A high-quality interface for reading from CSV files will be easy to implement after type polymorphism is implemented.)

Random data structures and Recursion

The iid distribution returns a list of random values from another distribution. We can apply the map and sum operations to such lists to sample a sum of squares.

module Demo2 where

import Distributions

main = do

  xs <- sample $ iid 10 (normal 0.0 1.0)

  let ys = map (\x -> x*x) xs

  return $ log_all [xs %% "xs", ys %% "squares", sum ys %% "sum"]

Here (\x -> x*x) describes an un-named function that takes an argument x and returns x*x.

(Currently the number 10 of i.i.d. normal variables cannot be random. Soon we will allow random numbers of random variables, and this restriction will be relaxed.)

Random control flow I: if-then-else

The modeling language can handle graphs that change. One thing that leads to a changing graph is control-flow statements like if-then-else:

module Demo3 where

import Distributions

main = do
  i <- sample $ bernoulli 0.5
  y <- sample $ normal 0.0 1.0
  let x = if (i==1) then y else 0.0
  return $ log_all [i %% "i", x %% "x"]

Random control flow II: arrays with random subscripts

Traditional graphical modelling languages, like BUGS, allow arrays of random variables. However, they do not allow selecting a random element of these arrays. Dynamic graphs allow random subscripts. This can be used to divide observations into categories. Here different elements of ys will be exactly equal, if they belong to the same category:

module Demo4 where

import Distributions

main = do

  xs <- sample $ iid 10 (normal 0.0 1.0)

  categories <- sample $ iid 10 (categorical (replicate 10 0.1))

  let ys = [xs!!(categories!!i) | i <- [0..9]]
  return $ log_all [ys %% "ys"]

Haskell uses the !! operator to subscript a list. The C equivalent of xs!!(categories!!i)) would be something like xs[categories[i]].

Random sampling via recursive functions: Brownian Bridge

We can use recursive functions to randomly sample lists of random values. Here, we define a function random_walk that produces a list of random values starting from x0.

module Demo5 where

import Distributions

random_walk x0 n f | n < 1     = error "Random walk needs at least 1 element"
                   | n == 1    = return [x0]
                   | otherwise = do x1 <- sample $ f x0
                                    xs <- random_walk x1 (n-1) f
                                    return (x0:xs)

-- 20 element brownian bridge
main = do
  zs <- random_walk 0.0 19 (\mu -> normal mu 1.0)

  observe 1.5 $ normal (last zs) 1.0

  return $ log_all [ zs %% "zs" ]
The argument f is a function. In Haskell, we write f x instead of f(x) to apply a function. Here, f x gives the distribution of the point after x.

The observe command specifies observed data. Here we observe that the next point after element 10 of zs is 1.5. This constrains the random walk to end at 1.5, creating a Brownian Bridge.


comments and suggestions: benjamin . redelings * gmail + com