Making BayesNets Bayesian

[itsdfish]

Hello,

I would like to use BayesNets.jl with Turing.jl to perform Bayesian parameter estimation. The primary requirement is that a library interacting with Turing is set up generically to support ForwardDiff.Dual types and Float64 in order to allow automatic differentiation to compute gradients. As shown below, I found that DiscreteCPD does not work because it is hardcoded to accept Float64. See this line.

A simple fix is to change that line to

const DiscreteCPD = CategoricalCPD

I tried adding constraints but it threw an error.

const DiscreteCPD = CategoricalCPD{Categorical{T,Vector{T}}} where {T<:Real}

I was wondering if you would be willing to change the line (I could also open a PR), and whether there are additional places that might need to be changed in the code.

Thanks!

MWE

using BayesNets, TikzGraphs
using ForwardDiff

c = ForwardDiff.Dual(.5,.4)
snc = ForwardDiff.Dual(.5,.4)
sc = ForwardDiff.Dual(.5,.4)

bn = BayesNet([DiscreteCPD(:c, [c,(1-c)])])
push!(bn, 
    DiscreteCPD(
        :s,  [:c], [2], 
    [Categorical([snc,(1-snc)]), 
    Categorical([sc,(1-sc)])]
    ),
)

Error

ERROR: TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 1}
Stacktrace:
  [1] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{Nothing, Float64, 1}, i1::Int64)
    @ Base ./array.jl:903
  [2] _unsafe_copyto!(dest::Vector{Float64}, doffs::Int64, src::Vector{ForwardDiff.Dual{Nothing, Float64, 1}}, soffs::Int64, n::Int64)
    @ Base ./array.jl:253
  [3] unsafe_copyto!
    @ ./array.jl:307 [inlined]
  [4] _copyto_impl!
    @ ./array.jl:331 [inlined]
  [5] copyto!
    @ ./array.jl:317 [inlined]
  [6] copyto!
    @ ./array.jl:343 [inlined]
  [7] copyto_axcheck!
    @ ./abstractarray.jl:1104 [inlined]
  [8] Vector{Float64}(x::Vector{ForwardDiff.Dual{Nothing, Float64, 1}})
    @ Base ./array.jl:563
  [9] convert
    @ ./array.jl:554 [inlined]
 [10] convert
    @ ~/.julia/packages/Distributions/1PkiH/src/univariate/discrete/discretenonparametric.jl:45 [inlined]
 [11] setindex!
    @ ./array.jl:903 [inlined]
 [12] _unsafe_copyto!(dest::Vector{Categorical{Float64, Vector{Float64}}}, doffs::Int64, src::Vector{Categorical{ForwardDiff.Dual{Nothing, Float64, 1}, Vector{ForwardDiff.Dual{Nothing, Float64, 1}}}}, soffs::Int64, n::Int64)
    @ Base ./array.jl:253
 [13] unsafe_copyto!
    @ ./array.jl:307 [inlined]
 [14] _copyto_impl!
    @ ./array.jl:331 [inlined]
 [15] copyto!
    @ ./array.jl:317 [inlined]
 [16] copyto!
    @ ./array.jl:343 [inlined]
 [17] copyto_axcheck!
    @ ./abstractarray.jl:1104 [inlined]
 [18] Vector{Categorical{Float64, Vector{Float64}}}(x::Vector{Categorical{ForwardDiff.Dual{Nothing, Float64, 1}, Vector{ForwardDiff.Dual{Nothing, Float64, 1}}}})
    @ Base ./array.jl:563
 [19] convert
    @ ./array.jl:554 [inlined]
 [20] DiscreteCPD(target::Symbol, parents::Vector{Symbol}, parental_ncategories::Vector{Int64}, distributions::Vector{Categorical{ForwardDiff.Dual{Nothing, Float64, 1}, Vector{ForwardDiff.Dual{Nothing, Float64, 1}}}})
    @ BayesNets.CPDs ~/.julia/packages/BayesNets/hZr0W/src/CPDs/categorical_cpd.jl:22
 [21] top-level scope
    @ Untitled-3:13

[mykelk]

That would be awesome! We certainly would welcome a PR!

[mykelk]

Please also add some tests to the test script as well. Thanks!

[itsdfish]

@mykelk, my apologies for the very belated reply. I never received a notification.

A strange thing happened. The code did not work with automatic differentiation when I tried it several years ago. The example above still doesn't work, and I was unable to get it to work after several changes to the source code. However, it runs with Turing.jl now. I am quite surprised and don't have an explanation. For future reference, here is an example:

Code

using BayesNets
using DataFrames
using Distributions 
using DynamicPPL: @addlogprob!
using Turing

import Distributions: logpdf
import Distributions: loglikelihood

struct BayesNetDist{T <: Real} <: DiscreteMultivariateDistribution
    c::T
    sc̄::T
    sc::T
    rc̄::T
    rc::T
    s̄r̄w::T
    s̄rw::T
    sr̄w::T
    srw::T
end

function BayesNetDist(; c, sc̄, sc, rc̄, rc, s̄r̄w, s̄rw, sr̄w, srw)
    return BayesNetDist(c, sc̄, sc, rc̄, rc, s̄r̄w, s̄rw, sr̄w, srw)
end

function generate_model(d::BayesNetDist)
    (; c, sc̄, sc, rc̄, rc, s̄r̄w, s̄rw, sr̄w, srw) = d
    return generate_model(; c, sc̄, sc, rc̄, rc, s̄r̄w, s̄rw, sr̄w, srw)
end

function logpdf(d::BayesNetDist, data::DataFrame)
    bn = generate_model(d)
    return logpdf(bn, data)
end

loglikelihood(d::BayesNetDist, data::DataFrame) =
    logpdf(d::BayesNetDist, data::DataFrame)

function generate_model(d::BayesNetDist)
    (; c, sc̄, sc, rc̄, rc, s̄r̄w, s̄rw, sr̄w, srw) = d
    bn = BayesNet([DiscreteCPD(:c, [c, (1-c)])])

    # Add conditional probability table for Cloudy
    push!(bn,
        CategoricalCPD(
            :s, [:c], [2],
            [Categorical([sc, (1-sc)]),
                Categorical([sc̄, (1-sc̄)])]
        ),
        # Add conditional probability table for Rain
        CategoricalCPD(
            :r, [:c], [2],
            [Categorical([rc, (1-rc)]),
                Categorical([rc̄, (1-rc̄)])]
        ),
        # Add conditional probability table for Wet
        CategoricalCPD(
            :w, [:s, :r], [2, 2],
            [Categorical([srw, (1-srw)]),
                Categorical([sr̄w, (1-sr̄w)]),
                Categorical([s̄rw, (1-s̄rw)]),
                Categorical([s̄r̄w, (1-s̄r̄w)])]
        )
    )
    return bn
end

@model function model(data)
    # define priors
    c ~ Beta(1, 1)
    sc̄ ~ Beta(1, 1)
    sc ~ Beta(1, 1)
    rc̄ ~ Beta(1, 1)
    rc ~ Beta(1, 1)
    s̄r̄w ~ Beta(1, 1)
    sr̄w ~ Beta(1, 1)
    s̄rw ~ Beta(1, 1)
    srw ~ Beta(1, 1)

    # likelihood of data
    @addlogprob! logpdf(BayesNetDist(; c, sc̄, sc, rc̄, rc, s̄r̄w, s̄rw, sr̄w, srw), data)
end

true_parms = (
    c = 0.5,
    sc = 0.10,
    sc̄ = 0.50,
    rc̄ = 0.20,
    rc = 0.80,
    srw = 0.99,
    sr̄w = 0.90,
    s̄rw = 0.90,
    s̄r̄w = 0.0
)
bn = generate_model(BayesNetDist(; true_parms...))
data = rand(bn, 200)

chains = sample(model(data), NUTS(1000, 0.65), MCMCThreads(), 1000, 4);

Version Info
Julia 1.11.7

  [ba4760a4] BayesNets v3.5.2
  [a93c6f00] DataFrames v1.8.0
  [31c24e10] Distributions v0.25.121
  [366bfd00] DynamicPPL v0.37.5
  [fce5fe82] Turing v0.40.3

[mykelk]

Wow! That's amazing! Thanks for sharing your example!

[itsdfish]

No problem! Its the famous sprinkler/rain/wet pavement example.

I'll close this issue, but if you think the example is useful for the documentation, I'd be happy to submit a PR.