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 Float64
s 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)