pri = NormalGamma(μ₀, κ₀, α₀, β₀)
post = posterior(pri, Normal, x)

I tried doing the same for a Bernoulli-Beta conjugate pair but didn’t work out:

julia> using Distributions, ConjugatePriors
julia> prior = Beta(10,10);
julia> poste = posterior(prior, Bernoulli, [0.2, 0.3]);
suffstats is not implemented for (Distributions.Bernoulli, Vector{Float64}).
1. error(::String)@error.jl:33
2. suffstats(::Type{Distributions.Bernoulli}, ::Vector{Float64})@genericfit.jl:5
3. posterior_canon(::Distributions.Beta{Float64}, ::Type{Distributions.Bernoulli}, ::Vector{Float64})@fallbacks.jl:3
4. posterior(::Distributions.Beta{Float64}, ::Type{Distributions.Bernoulli}, ::Vector{Float64})@fallbacks.jl:7
5. top-level scope@Local: 1[inlined]

Either the package is WIP or I miss something because I would consider Beta-Bernoulli to be very basic functionality.
Sadly there is no documentation in ConjugatePriors.jl which is a little bit frustrating.

It’s true that this package is a bit old and needs some love. IIRC @mschauer has started a PR to modernize it. As for me I hoped I would be able to contribute, but I have no time right now. In any case, you’re welcome to lend a hand

I can confirm that ReactiveMP has very strong support for accelerating computations when conjugate priors are available. Tagging @bvdmitri in case he has more to say about this.

MeasureTheory.jl also has the infrastructure to support this, but currently that’s the extent of it - such accelerations are not yet implemented.

@filchristou ReactiveMP.jl supports infrastructure to run very efficient Bayesian inference in models with conjugate pairs. A large table of rules has been implemented for many conjugate pairs (but not for all…). You can check-out this example on Beta-Bernoulli distributions: Getting Started · ReactiveMP.jl. On my laptop ReactiveMP.jl computes exact posterior with less than 1 seconds from dataset with 100_000 points:

julia> p
0.75
julia> dataset = float.(rand(rng, Bernoulli(p), n))
julia> length(dataset)
100000
julia> @btime inference(
model = Model($coin_model, length($dataset)),
data = (y = $dataset, )
)
693.036 ms (10499774 allocations: 465.60 MiB)
julia> result.posteriors[:θ][end]
Beta{Float64}(α=74987.0, β=25022.0)
julia> mean_var(result.posteriors[:θ][end])
(0.7498025177734003, 1.8757994411965801e-6)

If you are interested we have a list of examples available in the documentation and in the repository. Worth noting though that ReactiveMP supports not only conjugate pairs.