# Custom priors in Soss

**URL:** <https://discourse.julialang.org/t/custom-priors-in-soss/47076>\
**Category:** Probabilistic Programming\
**Created:** [September 22, 2020, 3:02pm UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076 "2020-09-22T15:02:50Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [September 22, 2020, 3:02pm UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076/1 "2020-09-22T15:02:50Z")

</div>

On Slack, Jonatan Werpers asks,

> How can I use a custom prior with a model? Basically I want to do something like:
> 
> ```julia
> my_pdf(θ) = ....
> m = @model n begin
> θ ~ my_pdf
> x ~ Bernoulli(θ) |> iid(n)
> end
> 
> ```

Generally, you’ll need three methods:

1. `Base.rand`
2. `Distributions.logpdf`
3. `Soss.xform`

For example, say we want the logistic transform of a Normal. We could do this entirely in Soss, but for now (due to Distribution.jl overhead), there may be a performance advantage to building it directly.

Starting out is pretty easy:

```julia
using Soss
import Base: rand

struct MyPrior end

Base.rand(::MyPrior) = logistic(randn())
Distributions.logpdf(MyPrior, x) = logistic(-x^2/2)

```

Note that there’s no need for the `logpdf` to be an _actual_ logpdf - normalization is irrelevant for many applications. We’ll soon be working instead in terms of _measures_, to avoid the extra computation when possible.

`Soss.xform` needs to return a transform in the sense of @Tamas_Papp’s [TransformVariables.jl](https://github.com/tpapp/TransformVariables.jl). For univariate cases this is especially simple:

```julia
Soss.xform(::MyPrior, _data=NamedTuple()) = as𝕀

```

The `_data` optional argument is just to make the dispatch work out, and it’s almost always fine to have it exactly like this.

Then we can do, e.g.,

```julia
julia> post = dynamicHMC(m(n=100), (x=x,));

julia> particles(post)
(θ = 0.245 ± 0.043,)

julia> m = @model n begin
           θ ~ MyPrior()
           x ~ Bernoulli(θ) |> iid(n)
       end;

julia> x = rand(100) .< 0.3;

julia> post = dynamicHMC(m(n=100), (x=x,));

julia> particles(post)
(θ = 0.314 ± 0.043,)

```

---

<div class="post-metadata">

**Author:** ![dilumaluthge](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dilumaluthge/32/29283_2.png) [@dilumaluthge](https://discourse.julialang.org/u/dilumaluthge)\
**Post date:** [September 22, 2020, 5:30pm UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076/2 "2020-09-22T17:30:06Z")

</div>

I would recommend implementing `Base.rand` as follows:

```julia
Base.rand(rng::Random.AbstractRNG, ::MyPrior) = logistic(randn(rng))

Base.rand(prior::MyPrior) = Base.rand(Random.GLOBAL_RNG, prior)

```

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [September 22, 2020, 5:31pm UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076/3 "2020-09-22T17:31:16Z")

</div>

Good point @dilumaluthge, that makes it much more flexible

---

<div class="post-metadata">

**Author:** ![werpers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/werpers/32/42610_2.png) [@werpers](https://discourse.julialang.org/u/werpers)\
**Post date:** [September 23, 2020, 6:58am UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076/4 "2020-09-23T06:58:42Z")

</div>

Thank you for your detailed answer!

After writing on slack I managed to get something working through a lot of trail and error. It’s not what you suggest but I think it’s working correctly.

What I came up with was the following

```julia
struct MyPrior <: Distribution{Univariate,Continuous} end
Base.minimum(::MyPrior) = 0
Base.maximum(::MyPrior) = 1
Base.rand(rng::AbstractRNG, ::MyPrior) = rejection_sampling(p,Uniform(0,1),3.5)
Distributions.pdf(::MyPrior, π) = exp(-1 + λ - b*(1-π)^3 - c*(1-π)^12)
Distributions.logpdf(::MyPrior, π::Real) = -1 + λ - b*(1-π)^3 - c*(1-π)^12

```

Since I’m pretty new to Bayesian statistics in general and completely new to PPLs I’d be very interested in hearing if there are any upsides or downsides to doing what I did compared to your solution, if it’s even correct.

From my limited knowledge of how MCMC and similar methods work I was a little bit surprised I needed to implement `rand()`. I thought the sampling methods only used the “shape” of the pdf when doing inference.

Also, what is the return `as𝕀` value in this function?

> [@cscherrer](#):
>
> ```julia
> Soss.xform(::MyPrior, _data=NamedTuple()) = as𝕀
> 
> ```

Again, thank you very much!

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [September 23, 2020, 1:56pm UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076/5 "2020-09-23T13:56:49Z")

</div>

> [@werpers](#):
>
> ```julia
> struct MyPrior <: Distribution{Univariate,Continuous} end
> Base.minimum(::MyPrior) = 0
> Base.maximum(::MyPrior) = 1
> Base.rand(rng::AbstractRNG, ::MyPrior) = rejection_sampling(p,Uniform(0,1),3.5)
> Distributions.pdf(::MyPrior, π) = exp(-1 + λ - b*(1-π)^3 - c*(1-π)^12)
> Distributions.logpdf(::MyPrior, π::Real) = -1 + λ - b*(1-π)^3 - c*(1-π)^12
> 
> ```

I guess you’re abbreviating here, right? Does `MyPrior` have `λ`, `b`, and `c` as parameters?

Otherwise, this looks good to me. You won’t need `pdf` for MCMC, but it can be nice to have around otherwise. Rather than copying the formula, I’d probably write it as

```julia
Distributions.pdf(dist::MyPrior, π) = exp(logpdf(dist, π))

```

There’s no overhead, and it will propagate any changes instead of needing to type them twice,

> [@werpers](#):
>
> From my limited knowledge of how MCMC and similar methods work I was a little bit surprised I needed to implement `rand()` . I thought the sampling methods only used the “shape” of the pdf when doing inference.

This is just a hack because Distributions.jl has no way of representing the type of the support. For example,

```julia
julia> d = MvNormal(ones(3))
ZeroMeanDiagNormal(
dim: 3
μ: 3-element FillArrays.Zeros{Float64,1,Tuple{Base.OneTo{Int64}}} = 0.0
Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
)

julia> typeof(d)
MvNormal{Float64,PDMats.PDiagMat{Float64,Array{Float64,1}},FillArrays.Zeros{Float64,1,Tuple{Base.OneTo{Int64}}}}

```

We need to know that `rand(d)` returns a Float64 vector of length 3. But there’s not a three to be found! So for now I just draw a `rand` and use that for the analysis. This will get better with [MeasureTheory.jl](https://github.com/cscherrer/MeasureTheory.jl), which we’ll be switching to when it’s mature enough.

> [@werpers](#):
>
> Also, what is the return `as𝕀` value in this function?

`as𝕀` is from TransformVariables.jl. We need to transform to a density over ℝⁿ. This will go in as one component, and all the pieces compose to get the final result. Like this:

```julia
julia> using TransformVariables
[Info: Precompiling TransformVariables [84d833dd-6860-57f9-a1a7-6da5db126cff]

julia> f = as𝕀
as𝕀

julia> f(-30.0)
9.357622968839299e-14

julia> f(20.0)
0.9999999979388463

julia> inverse(f)(0.2)
-1.3862943611198906

julia> inverse(f)(0.99)
4.595119850134589

```

---

<div class="post-metadata">

**Author:** ![werpers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/werpers/32/42610_2.png) [@werpers](https://discourse.julialang.org/u/werpers)\
**Post date:** [September 23, 2020, 8:18pm UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076/6 "2020-09-23T20:18:15Z")

</div>

Yes, I’m abbreviating, `λ` , `b` , and `c` are all constants.

The `rand` thing indeed seems like a little bit of a hack. Could you solve it by using the `length` function?

```julia
julia> d = MvNormal(ones(3))
ZeroMeanDiagNormal(
dim: 3
μ: 3-element FillArrays.Zeros{Float64,1,Tuple{Base.OneTo{Int64}}} = 0.0
Σ: [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
)

julia> length(d)
3

```

Or perhaps you need it to be more “static” in some sense?

Thanks again, I’ll be keeping an eye on these packages to learn more!

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [September 23, 2020, 8:23pm UTC](https://discourse.julialang.org/t/custom-priors-in-soss/47076/7 "2020-09-23T20:23:40Z")

</div>

> [@werpers](#):
>
> The `rand` thing indeed seems like a little bit of a hack.

Yep, it’s entirely a hack.

> [@werpers](#):
>
> Could you solve it by using the `length` function?

No, because there are plenty of values that aren’t arrays.

> [@werpers](#):
>
> Or perhaps you need it to be more “static” in some sense?

Yes! That’s the plan, in MeasureTheory.jl. Similar to how Gen goes about it, we’ll have a type parameter that specifies the return value. For “built-in” measures, we’ll probably use `SizedArray`s so we can get the shapes right without the extra overhead and annoyance 🙂