Informally,

Omega comes with a small number of built-in primitive probability distributions. One example is the standard uniform:

using Omega
xs = StdUniform{Float64}()

xs represents a collection of random variables, which we call a random variable class. To get the nth random variable in a random variable class we use the nth function

x1 = nth(xs, 1)

A more convenient way to write this is using ~

x1 = 1 ~ xs

Random variables, such as x1 represent quantities whose true value is uncertain. There are many ways in English to express

All of the probability distributions in the Distributions.jl package are random variable classes.

using  Distributions
x1 = @~ Uniform(0.0, 1.0)

sample. To construct another random variable x2, we do the same.

x2 = @~ Uniform(0.0, 1.0)

x1 and x2 are identically distributed and independent (i.i.d.).

julia> randsample((x1, x2))
(0.5602978842341093, 0.9274576159629635)

Contrast this with:

julia> randsample((x1, x1))
(0.057271529749001626, 0.057271529749001626)

Composition

A single primitive random variable by itself isn't very interesting or useful. In order to model complex things in the world we will typically need to compose together random variables into more complex ones.

There are two ways to compose random variables: the statistical style, which can be less verbose, and more intuitive, but has some limitations, and the explicit style, which is more general.

Statistical style is convenient because it allows us to treat a T-valued random variable as if it is a value of type T. For instance @~ Uniform(0.0, 1.0) is Float64-valued random variable, and hence using the statistical style, we can add, multiply, divide them as if they were values of type Float64.

To add two random variables we use the dot-form of the + function:

x3 = x1 .+ x2
randsample(@joint(x1, x2, x3))
(x1 = 0.3940869955915858, x2 = 0.5380135911080237, x3 = 0.9321005866996095)

x3 is a random variable.

We can apply inequalities in the same way:

p = x3 .> 1.0
julia> randsample(@joint(x3, p))
(x3 = 1.2748335833273177, p = true)

Random variable families which take parameters of type T can in the same way take T-valued random variables.

μ = @~ Normal(0, 1)
n = @~ Normal.(μ, 1)

Suppose you write your own function which take Float64s as input:

myfunc(x::Float64, y::Float64) = (x * y)^2

We can't apply myfunc to random variables; it will cause a method error

julia> myfunc(x1, x2)
ERROR: MethodError: no method matching myfunc...

However this is easily remedied with dot notation

myfunc.(x1, x2)

Explicit Style

The above style is convenient but has a few limitations. For example, we can't use it with Julia's conditional branching constructs like if

if p
  3
else
  4
end
ERROR: TypeError: non-boolean (OmegaCore.Var.Pw{Tuple{OmegaCore.Var.Pw{Tuple{Member{Uniform{Float64}, Int64}, Member{Uniform{Float64}, Int64}}, typeof(+)}, Float64}, typeof(>)}) used in boolean context

To create random variables in the explicit style, create a normal Julia function that takes as input ω

For instance, to define a bernoulli distribution in explicit style:

x(ω::Ω) = @~ StdNormal{Float64}()(ω) > 0.5

x is just a normal julia function.

All of the primitive distributions can be used in explicit style by passing the ω object as the first parameter (type constraints are added just to show that the return values are not random variables but elements. But don't add them to your own code! It will prevent automatic differentiation based inference procedures from working):

function x(ω)
  if (@~ Bernoulli(0.5))(ω)::Bool
    (@~ Normal(0, 1))(ω)::Float64
  else (@~ Bernoulli(0.5))(ω)::Bool
    (@~ Beta(2.0, 2.0))(ω)::Float64
  end
end

Statistical style and functional style can be combined naturally. For example:

x(ω) = (@~ Bernoulli(ω)) > 0.5 ? StdUniform(ω)^2 : sqrt.(StdUniform(ω))
y = Normal(0.0, 1.0)
z = x .+ y

Random Variable Families

Often we want to parameterize a random variable. To do this we create functions which return random variables as a function of some parameters. For example, we can make a Uniform distribution family (without using Distributions.jl) by defining a function which maps a and b to a random variable

"Uniform distribution between `a` and `b`"

unif(a, b) = StdUniform{Float64}() .* (b - a) + b  


# x is uniformly distributed between 10 and 20
x = unif(10, 20)

And hence if we wanted to create a method that created independent uniformly distributed random variables, we could do it like so:

unif2(a,b) =~ rng -> rand(rng) * (b - a) + a

# x is distributed between 30 and 40 (and independent of y)
x = unif2(30, 40)

# y is distributed between 30 and 40 (and independent of x)
y = unif2(30, 40)