Chuch.jl aims to make it easy for anyone to perform MCMC inference in complex, and simple, probabilistic models. We aim to be:

- Fast
- Church.jl considers only the bits of the model that need to be updated during MCMC.

- Practical
- Easy to install and use.

- General
- Define new, arbitrarily complex distributions, within Church.jl.
- Allows you to define stochastic, recursive, and non-parametric distributions.
- Uses all distributions available in Distributions.jl.
- Includes a probabilistic garbage collector - so you define large models without running out of memory.

To install, use

```
Pkg.add("Church")
```

To load, use

```
using Church
```

You can define a random variable simply by calling a function,

```
a = normal()
```

Church.jl contains all the distributions in Distributions.jl - with function names that are just lowercase versions of the distribution names.

You can apply standard operators (e.g. `+`

) directly to samples,

```
b = normal() * normal()
```

To apply other functions to random variables, you must first "lift" the function. This overloads the function, so that it can deal with random variables. For instance,

```
@lift(cosh, 1)
c = cosh(normal())
```

Note that the second argument to lift is the number of arguments, which, in this case, is 1.

Combining these, we could write,

```
a = normal()
b = a * a
@lift(cosh, 1)
c = cosh(a)
d = normal(a, c)
```

To sample these random variables, we use `resample()`

, which performs a single MCMC step, and use `value(a)`

to report the value of `a`

for the current sample.
Note that `value(a)`

is ONLY provided for recording or printing the value of samples.
Any other is liable to give meaningless quantities in the best case, or cause the algorithm to no longer sample the correct distribution in the worst case.

```
for i = 1:5
#Do 5 MCMC steps.
resample(5)
#value(a) returns the value of a for the current sample.
@printf("a:% .3f, b:% .3f, c:% .3f, d:% .3f", value(a), value(b), value(c), value(d)); println()
end
#Prints:
#a: 0.978, b: 0.957, c: 1.518, d:-0.240
#a: 0.262, b: 0.069, c: 1.035, d:-1.598
#a: 0.262, b: 0.069, c: 1.035, d:-1.065
#a: 0.262, b: 0.069, c: 1.035, d:-0.458
#a:-0.182, b: 0.033, c: 1.017, d:-0.112
```

So far, we haven't done anything interesting - you could sample `a`

and `b`

in the previous sections by simply using `a = rand(Normal(0, 1))`

and `b = a*a`

.
In Church.jl, you can condition these draws on known data, using the keyword argument `condition`

,

```
normal(1, 1; condition=3)
```

In a more complete example, to sample P(a| c=10), where c ~ Normal(b, 0.1), and a ~ Normal(0, 1),

```
using Church
@lift(abs, 1)
a = normal(0, 1)
b = abs(a)
c = normal(b, 0.1; condition=3)
#Now that we're doing inference, we need to perform many sampling steps,
#for the model to converge to the correct distribution. This is known as burn-in.
resample(10^3)
println(value(a))
#Prints:
#-2.7382930822004345
```

```
using Church
#Generate some data
data = [randn(10)+6, randn(10)-6]
#The model parameters
K = 2
ms = [normal(0, 10) for i = 1:K]
vs = [gamma(2, 2) for i = 1:K]
ps = dirichlet(ones(K))
#Which mixture component does each data item belong to?
ks = [categorical(ps) for i = 1:length(data)]
for i = 1:length(data)
#Condition on the data.
normal(ms[ks[i]], vs[ks[i]]; condition=data[i])
end
resample(10^4)
@printf("m1:% .3f, m2:% .3f, v1:% .3f, v2:% .3f", value(ms[1]), value(ms[2]), value(vs[1]), value(vs[2]))
println()
map(x -> print(value(x)), ks)
println()
println((value(ps)[1], value(ps)[2]))
#Prints:
#m1:-5.575, m2: 6.024, v1: 0.947, v2: 0.940
#22222222221111111111
#(0.43393275606877524,0.5660672439312248)
```

Note that the inferred parameters are sensible, given the data.

We might want to define a mixture model with a variable number of components, for instance,

```
K = poisson(3)
ms = [normal(0, 10) for i = 1:K]
```

However, you cannot do this, because the list comprehension needs K to be an integer, not a sample. Instead, you can use a large number of components, then exploit the lazy datastructures and garbage collector in Church.jl to avoid instantiating unused mixture components. For instance,

```
using Church
using Distributions
#Generate some data
data = [randn(10)+6, randn(10)-6]
#The model parameters
ms = Mem((i::Int) -> normal(0, 10))
vs = Mem((i::Int) -> gamma(2, 2))
ps = dirichlet(10, 1.; sampler=(d,v)->Dirichlet(3*v+0.001))
#Which mixture component does each data item belong to?
ks = [categorical(ps) for i = 1:length(data)]
#Condition on the data.
for i = 1:length(data)
normal(ms[ks[i]], vs[ks[i]]; condition=data[i])
end
for i = 1:10^3
resample(10^3)
gc_church()
end
map(x -> print(value(x)), ks)
#Prints:
#89998788884444444444
```

So the model is only using 4 components.
Looking at the value of `ms`

, we see that the parameters for the other components have not been instansiated.
The other components, that have been created at some point during the sampling, have been cleaned up.

In Church.jl, a distribution is just a function, so to define a new distributioon, we just need to define a function. A very simple example, is a mixture of normal distributions with different standard deviations,

```
gsm(m::Real) = normal(m, gamma(2, 2))
```

We can also allow conditioning on new distributions, if the final call also allows conditioning. In this case,

```
gsm(m::Real; condition=nocond) = normal(m, gamma(2, 2); condition=condition)
```

Now `gsm`

can be conditioned just like any other distribution. Note that you should use nocond as the default value of condition - this is the special value indicating that the distribution is not conditioned.

In another example, we could use recursion to write down a distribution,

```
geom(p::Real) = @If(bernoulli(p), 1+geom(), 1).
```

The macro `@If`

returns `1+geom()`

if `bernoulli(p)`

is `1`

, and returns `1`

if `bernoulli(p)`

is `0`

.
It only evaluates its arguments as necassery, so we do not get infinite recursive calls to `geom`

.

Finally, you can write down random distributions.
For instance, the dirichlet distribution can be thought of as returning a vector, whose elements are positive and sum to 1, or it can be thought of as returning a categorical distribution.
The `dirichlet`

distribution returns a vector, in Church.jl.
However, we could define `fdirichlet`

, which does return a distribution,

```
fdirichlet(args...) = begin
ps = dirichlet(args...)
() -> categorical(ps)
end
#dir is a categorical distribution
dir = fdirichlet(9, 0.1)
for i = 1:10
print(value(dir()))
end
#Prints
#4444448444
```

This mechanism allows you to write down a dirichlet process, which, again, is a distribution over random distributions,

```
dp(concentration::Real, base_measure::Function) = begin
sticks = Mem(i::Int -> beta(1., concentration))
atoms = Mem(i::Int -> base_measure())
loop(i::Int) =
@If(bernoulli(sticks[i]), atoms[i], loop(i+1))
d = () -> loop(1)
end
#d is a DP
d = dp(1., normal)
for i = 1:10
println(value(d()))
end
#Prints
#-0.8241868921220118
#1.2419259288985225
#0.7481036918602126
#1.4163485859423797
#1.2419259288985225
#1.2419259288985225
#1.2419259288985225
#1.2419259288985225
#1.2419259288985225
#1.2419259288985225
```

12/27/2013

about 2 years ago

50 commits