Basic Tutorial

Basic Tutorial

In this tutorial we will run through the basics of creating a model and conditioning it. This tutorial is available in notebook version.

First load Omega:

using Omega

If you tossed a coin and observed the sequqnce HHHHH, you would be a little suspicious, HHHHHHHH would make you very suspicious. Elementary probability theory tells us that for a fair coin, HHHHHHHH is just a likely outcome as HHTTHHTH. What gives? We will use Omega to model this behaviour, and see how that belief about a coin changes after observing a number of tosses.

Model the coin as a bernoulli distribution. The weight of a bernoulli determines the probability it comes up true (which represents heads). Use a beta distribution to represent our prior belief weight of the coin.

weight = β(2.0, 2.0)

A beta distribution is appropriate here because it is bounded between 0 and 1.

Draw a 10000 samples from weight using rand:

beta_samples = rand(weight, 10000)

Let's see what this distribution looks like using UnicodePlots. If you don't have it installed already install with:

] add UnicodePlots

To visualize the distribution, plot a histogram of the samples:

using UnicodePlots
   (0.0,0.1] │▇▇▇▇▇▇ 279                              │ 
   (0.1,0.2] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 727                   │ 
   (0.2,0.3] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1218       │ 
   (0.3,0.4] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1354    │ 
   (0.4,0.5] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1482 │ 
   (0.5,0.6] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1426  │ 
   (0.6,0.7] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1406   │ 
   (0.7,0.8] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 1124         │ 
   (0.8,0.9] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 702                    │ 
   (0.9,1.0] │▇▇▇▇▇▇ 282                              │ 

The distribution is symmetric around 0.5 and has support over the the interval [0, 1].

So far we have not done anything we couldn't do with Distributions.jl. A primary distinction between a package like Distribution.jl, is that Omega.jl allows you to condition probability distributions.

Create a model representing four flips of the coin. Since a coin can be heads or tales, the appropriate distribution is the bernouli distribution:

nflips = 4
coinflips_ = [bernoulli(weight, Bool) for i = 1:nflips]

Take note that weight is the random variable defined previously. bernoulli takes a type as its secoond argument; Bool indicates the result will be a Bool rather than an Int.

coinflips is a normal Julia array of Random Variables (RandVars). For reasons we will elaborate in later sections, it will be useful to have an Array-valued RandVar (instead of an Array of RandVar).

One way to do this (there are several ways discuseed later), is to use the function randarray

coinflips = randarray(coinflips_)

coinflips is a RandVar and hence we can sample from it with rand

julia> rand(coinflips)
4-element Array{Bool,1}:

Now we can condition the model. We want to find the conditional distribution over the weight of the coin given some observations.

First create some fake data

observations = [true, true, true, false]

Create a predicate that tests whether simulating from the model matches the observed data:

condition = coinflips ==ᵣ observations

condition is a random variable; we can sample from it. The function ==ᵣ (and more generally functions subscripted with ᵣ) should be read as "a realization of coinflips == observations"

We can use rand to sample from the model conditioned on condition being true:

weight_samples = rand(weight, condition, 10; alg = RejectionSample)

weight_samples is a set of 10 samples from the conditional (sometimes called posterior) distribution of weight condition on the fact that coinflips == observations.

In this case, rand takes

Plot a histogram of the weights like before:

julia> UnicodePlots.histogram(weight_samples)
   (0.1,0.2] │▇ 4                                     │ 
   (0.2,0.3] │▇▇▇ 22                                  │ 
   (0.3,0.4] │▇▇▇▇▇▇▇▇▇▇▇ 69                          │ 
   (0.4,0.5] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 147             │ 
   (0.5,0.6] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 185       │ 
   (0.6,0.7] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 226 │ 
   (0.7,0.8] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 203     │ 
   (0.8,0.9] │▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇▇ 120                 │ 
   (0.9,1.0] │▇▇▇▇ 23                                 │ 

Observe that our belief about the weight has now changed. We are more convinced the coin is biased towards heads (true).