# 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:

**Functions**work, and can be used to define random variables.**Modules**work, and allow code to be factored in clean manner.**Packages**work, and allow researchers to distribute their work separately from the BAli-Phy architecture.**Optimization**works, and speeds up the code written by the user by using techniques such as inlining.**Composite Objects**work, and can be used to define random data structures.**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.