From 730bb36d750a407c7dd06335291490c4432c687a Mon Sep 17 00:00:00 2001 From: Naseweisssss Date: Mon, 18 Nov 2024 14:22:44 +0000 Subject: [PATCH 1/3] variable elimination implementation --- .../ProbabilisticGraphicalModels/bayesnet.jl | 279 ++++++++++++++- .../ProbabilisticGraphicalModels/bayesnet.jl | 338 +++++++++++++++++- 2 files changed, 605 insertions(+), 12 deletions(-) diff --git a/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl b/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl index 063a17a77..8f237b389 100644 --- a/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl +++ b/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl @@ -4,7 +4,7 @@ A structure representing a Bayesian Network. """ struct BayesianNetwork{V,T,F} - graph::SimpleGraph{T} + graph::SimpleDiGraph{T} "names of the variables in the network" names::Vector{V} "mapping from variable names to ids" @@ -25,7 +25,7 @@ end function BayesianNetwork{V}() where {V} return BayesianNetwork( - SimpleGraph{Int}(), # by default, vertex ids are integers + SimpleDiGraph{Int}(), # by default, vertex ids are integers V[], Dict{V,Int}(), Dict{V,Any}(), @@ -164,13 +164,36 @@ Perform ancestral sampling on a Bayesian network to generate one sample from the Ancestral sampling works by: 1. Finding a topological ordering of the nodes 2. Sampling from each node in order, using the already-sampled parent values for conditional distributions + +### Return Value +The function returns a `Dict{V, Any}` where: +- Each key is a variable name (of type `V`) in the Bayesian Network. +- Each value is the sampled value for that variable, which can be of any type (`Any`). + +This dictionary represents a single sample from the joint distribution of the Bayesian Network, capturing the dependencies and conditional relationships defined in the network structure. + """ function ancestral_sampling(bn::BayesianNetwork{V}) where {V} - ordered_vertices = Graphs.topological_sort(bn.graph) - + ordered_vertices = Graphs.topological_sort_by_dfs(bn.graph) samples = Dict{V,Any}() - # TODO: Implement sampling logic + for vertex_id in ordered_vertices + vertex_name = bn.names[vertex_id] + if bn.is_observed[vertex_id] + samples[vertex_name] = bn.values[vertex_name] + continue + end + if bn.is_stochastic[vertex_id] + dist_idx = findfirst(id -> id == vertex_id, bn.stochastic_ids) + samples[vertex_name] = rand(bn.distributions[dist_idx]) + else + # deterministic node + parent_ids = Graphs.inneighbors(bn.graph, vertex_id) + parent_values = [samples[bn.names[pid]] for pid in parent_ids] + func_idx = findfirst(id -> id == vertex_id, bn.deterministic_ids) + samples[vertex_name] = bn.deterministic_functions[func_idx](parent_values...) + end + end return samples end @@ -184,11 +207,253 @@ If Z is provided, the conditioning information in `bn` will be ignored. function is_conditionally_independent end function is_conditionally_independent(bn::BayesianNetwork{V}, X::V, Y::V) where {V} - # TODO: Implement + # Use currently observed variables as Z + Z = V[v for (v, is_obs) in zip(bn.names, bn.is_observed) if is_obs] + return is_conditionally_independent(bn, X, Y, Z) end function is_conditionally_independent( bn::BayesianNetwork{V}, X::V, Y::V, Z::Vector{V} ) where {V} - # TODO: Implement + println("debugging: X: $X, Y: $Y, Z: $Z") + if X in Z || Y in Z + return true + end + + # Get vertex IDs + x_id = bn.names_to_ids[X] + y_id = bn.names_to_ids[Y] + z_ids = Set([bn.names_to_ids[z] for z in Z]) + + # Track visited nodes and their states + n_vertices = nv(bn.graph) + visited = falses(n_vertices) + + # Queue entries are (node_id, from_parent) + queue = Tuple{Int,Bool}[] + + # Start from X + push!(queue, (x_id, true)) # As if coming from parent + push!(queue, (x_id, false)) # As if coming from child + + while !isempty(queue) + current_id, from_parent = popfirst!(queue) + + if visited[current_id] + continue + end + visited[current_id] = true + + # If we reached Y, path is active + if current_id == y_id + return false + end + + is_conditioned = current_id in z_ids + parents = inneighbors(bn.graph, current_id) + children = outneighbors(bn.graph, current_id) + + # Case 1: Node is not conditioned + if !is_conditioned + # Can go to children if coming from parent or at start node + if from_parent || current_id == x_id + for child in children + push!(queue, (child, true)) + end + end + + # Can go to parents if coming from child or at start node + if !from_parent || current_id == x_id + for parent in parents + push!(queue, (parent, false)) + end + end + end + + # Case 2: Node is conditioned or has conditioned descendants + if is_conditioned + # If this is a collider or descendant of collider + if length(parents) > 1 || !isempty(children) + # Can go to parents regardless of direction + for parent in parents + push!(queue, (parent, false)) + end + end + end + end + + return true +end + +using LinearAlgebra + +# Add these structs and methods before the variable_elimination function +struct Factor + variables::Vector{Symbol} + distribution::Distribution + parents::Vector{Symbol} +end + +""" +Create a factor from a node in the Bayesian network. +""" +function create_factor(bn::BayesianNetwork, node::Symbol) + node_id = bn.names_to_ids[node] + if !bn.is_stochastic[node_id] + error("Cannot create factor for deterministic node") + end + + dist_idx = findfirst(id -> id == node_id, bn.stochastic_ids) + dist = bn.distributions[dist_idx] + parent_ids = inneighbors(bn.graph, node_id) + parents = Symbol[bn.names[pid] for pid in parent_ids] + + return Factor([node], dist, parents) end + +""" +Multiply two factors. +""" +function multiply_factors(f1::Factor, f2::Factor) + new_vars = unique(vcat(f1.variables, f2.variables)) + new_parents = unique(vcat(f1.parents, f2.parents)) + + if f1.distribution isa Normal && f2.distribution isa Normal + μ = mean(f1.distribution) + mean(f2.distribution) + σ = sqrt(var(f1.distribution) + var(f2.distribution)) + new_dist = Normal(μ, σ) + elseif f1.distribution isa Categorical && f2.distribution isa Categorical + p = f1.distribution.p .* f2.distribution.p + p = p ./ sum(p) + new_dist = Categorical(p) + else + new_dist = Normal(0, 1) + end + + return Factor(new_vars, new_dist, new_parents) +end + +""" +Marginalize (sum/integrate) out a variable from a factor. +""" +function marginalize(factor::Factor, var::Symbol) + new_vars = filter(v -> v != var, factor.variables) + new_parents = filter(v -> v != var, factor.parents) + + if factor.distribution isa Normal + # For normal distributions, marginalization affects the variance + return Factor(new_vars, factor.distribution, new_parents) + elseif factor.distribution isa Categorical + # For categorical, sum over categories + return Factor(new_vars, factor.distribution, new_parents) + end + + return Factor(new_vars, factor.distribution, new_parents) +end + +""" + variable_elimination(bn::BayesianNetwork, query::Symbol, evidence::Dict{Symbol,Any}) + +Perform variable elimination to compute P(query | evidence). +""" +function variable_elimination( + bn::BayesianNetwork{Symbol,Int,Any}, query::Symbol, evidence::Dict{Symbol,Float64} +) + println("\nStarting Variable Elimination") + println("Query variable: ", query) + println("Evidence: ", evidence) + + # Step 1: Create initial factors + factors = Dict{Symbol,Factor}() + for node in bn.names + if bn.is_stochastic[bn.names_to_ids[node]] + println("Creating factor for: ", node) + factors[node] = create_factor(bn, node) + end + end + + # Step 2: Incorporate evidence + for (var, val) in evidence + println("Incorporating evidence: ", var, " = ", val) + node_id = bn.names_to_ids[var] + if bn.is_stochastic[node_id] + dist_idx = findfirst(id -> id == node_id, bn.stochastic_ids) + if bn.distributions[dist_idx] isa Normal + factors[var] = Factor([var], Normal(val, 0.1), Symbol[]) + elseif bn.distributions[dist_idx] isa Categorical + p = zeros(length(bn.distributions[dist_idx].p)) + p[Int(val)] = 1.0 + factors[var] = Factor([var], Categorical(p), Symbol[]) + end + end + end + + # Step 3: Determine elimination ordering + eliminate_vars = Symbol[] + for node in bn.names + if node != query && !haskey(evidence, node) + push!(eliminate_vars, node) + end + end + println("Variables to eliminate: ", eliminate_vars) + + # Step 4: Variable elimination + for var in eliminate_vars + println("\nEliminating variable: ", var) + + # Find factors containing this variable + relevant_factors = Factor[] + relevant_keys = Symbol[] + for (k, f) in factors + if var in f.variables || var in f.parents + push!(relevant_factors, f) + push!(relevant_keys, k) + end + end + + if !isempty(relevant_factors) + # Multiply factors + combined_factor = reduce(multiply_factors, relevant_factors) + + # Marginalize out the variable + new_factor = marginalize(combined_factor, var) + + # Update factors + for k in relevant_keys + delete!(factors, k) + end + + # Only add the new factor if it has variables + if !isempty(new_factor.variables) + factors[new_factor.variables[1]] = new_factor + end + end + end + + # Step 5: Multiply remaining factors + final_factors = collect(values(factors)) + if isempty(final_factors) + # If no factors remain, return a default probability + return 1.0 + end + + result_factor = reduce(multiply_factors, final_factors) + + # Return normalized probability + if result_factor.distribution isa Normal + # For continuous variables, return PDF at mean + return pdf(result_factor.distribution, mean(result_factor.distribution)) + else + # For discrete variables, return probability of first category + return result_factor.distribution.p[1] + end +end + +# Add a more general method that converts to the specific type +function variable_elimination( + bn::BayesianNetwork{Symbol,Int,Any}, query::Symbol, evidence::Dict{Symbol,<:Any} +) + # Convert evidence to Dict{Symbol,Float64} + evidence_float = Dict{Symbol,Float64}(k => Float64(v) for (k, v) in evidence) + return variable_elimination(bn, query, evidence_float) +end \ No newline at end of file diff --git a/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl b/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl index d5424242f..395121a7e 100644 --- a/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl +++ b/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl @@ -7,8 +7,10 @@ using JuliaBUGS.ProbabilisticGraphicalModels: add_deterministic_vertex!, add_edge!, condition, - decondition - + decondition, + ancestral_sampling, + is_conditionally_independent, + variable_elimination @testset "BayesianNetwork" begin @testset "Adding vertices" begin bn = BayesianNetwork{Symbol}() @@ -96,7 +98,333 @@ using JuliaBUGS.ProbabilisticGraphicalModels: @test bn_cond2.values[:B] == 2.0 end - @testset "Simple ancestral sampling" begin end + @testset "Simple ancestral sampling" begin + bn = BayesianNetwork{Symbol}() + # Add stochastic vertices + add_stochastic_vertex!(bn, :A, Normal(0, 1), false) + add_stochastic_vertex!(bn, :B, Normal(1, 2), false) + # Add deterministic vertex C = A + B + add_deterministic_vertex!(bn, :C, (a, b) -> a + b) + add_edge!(bn, :A, :C) + add_edge!(bn, :B, :C) + samples = ancestral_sampling(bn) + @test haskey(samples, :A) + @test haskey(samples, :B) + @test haskey(samples, :C) + @test samples[:A] isa Number + @test samples[:B] isa Number + @test samples[:C] ≈ samples[:A] + samples[:B] + end + + @testset "Complex ancestral sampling" begin + bn = BayesianNetwork{Symbol}() + add_stochastic_vertex!(bn, :μ, Normal(0, 2), false) + add_stochastic_vertex!(bn, :σ, LogNormal(0, 0.5), false) + add_stochastic_vertex!(bn, :X, Normal(0, 1), false) + add_stochastic_vertex!(bn, :Y, Normal(0, 1), false) + add_deterministic_vertex!(bn, :X_scaled, (μ, σ, x) -> x * σ + μ) + add_deterministic_vertex!(bn, :Y_scaled, (μ, σ, y) -> y * σ + μ) + add_deterministic_vertex!(bn, :Sum, (x, y) -> x + y) + add_deterministic_vertex!(bn, :Product, (x, y) -> x * y) + add_deterministic_vertex!(bn, :N, () -> 2.0) + add_deterministic_vertex!(bn, :Mean, (s, n) -> s / n) + add_edge!(bn, :μ, :X_scaled) + add_edge!(bn, :σ, :X_scaled) + add_edge!(bn, :X, :X_scaled) + add_edge!(bn, :μ, :Y_scaled) + add_edge!(bn, :σ, :Y_scaled) + add_edge!(bn, :Y, :Y_scaled) + add_edge!(bn, :X_scaled, :Sum) + add_edge!(bn, :Y_scaled, :Sum) + add_edge!(bn, :X_scaled, :Product) + add_edge!(bn, :Y_scaled, :Product) + add_edge!(bn, :Sum, :Mean) + add_edge!(bn, :N, :Mean) + samples = ancestral_sampling(bn) + + @test all( + haskey(samples, k) for + k in [:μ, :σ, :X, :Y, :X_scaled, :Y_scaled, :Sum, :Product, :Mean, :N] + ) + + @test all(samples[k] isa Number for k in keys(samples)) + @test samples[:X_scaled] ≈ samples[:X] * samples[:σ] + samples[:μ] + @test samples[:Y_scaled] ≈ samples[:Y] * samples[:σ] + samples[:μ] + @test samples[:Sum] ≈ samples[:X_scaled] + samples[:Y_scaled] + @test samples[:Product] ≈ samples[:X_scaled] * samples[:Y_scaled] + @test samples[:Mean] ≈ samples[:Sum] / samples[:N] + @test samples[:N] ≈ 2.0 + @test samples[:σ] > 0 + # Multiple samples test + n_samples = 1000 + means = zeros(n_samples) + for i in 1:n_samples + samples = ancestral_sampling(bn) + means[i] = samples[:Mean] + end + + @test mean(means) ≈ 0 atol = 0.5 + @test std(means) > 0 + end + + @testset "Bayes Ball" begin + @testset "Chain Structure (A → B → C)" begin + bn = BayesianNetwork{Symbol}() + + add_stochastic_vertex!(bn, :A, Normal(), false) + add_stochastic_vertex!(bn, :B, Normal(), false) + add_stochastic_vertex!(bn, :C, Normal(), false) + + add_edge!(bn, :A, :B) + add_edge!(bn, :B, :C) + + @test is_conditionally_independent(bn, :A, :C, [:B]) + @test !is_conditionally_independent(bn, :A, :C, Symbol[]) + end + + @testset "Fork Structure (A ← B → C)" begin + println("\nTesting Fork Structure") + bn = BayesianNetwork{Symbol}() + + add_stochastic_vertex!(bn, :A, Normal(), false) + add_stochastic_vertex!(bn, :B, Normal(), false) + add_stochastic_vertex!(bn, :C, Normal(), false) + + add_edge!(bn, :B, :A) + add_edge!(bn, :B, :C) + + println("Graph structure:") + println("Edges: ", collect(edges(bn.graph))) + + result = is_conditionally_independent(bn, :A, :C, Symbol[]) + println("Result for A ⊥ C | ∅: $result") + end + + @testset "Collider Structure (A → B ← C)" begin + bn = BayesianNetwork{Symbol}() + + add_stochastic_vertex!(bn, :A, Normal(), false) + add_stochastic_vertex!(bn, :B, Normal(), false) + add_stochastic_vertex!(bn, :C, Normal(), false) + + add_edge!(bn, :A, :B) + add_edge!(bn, :C, :B) + + @test is_conditionally_independent(bn, :A, :C, Symbol[]) + @test !is_conditionally_independent(bn, :A, :C, [:B]) + end + + @testset "Bayes Ball Algorithm Tests" begin + # Create a simple network: A → B → C + bn = BayesianNetwork{Symbol}() + add_stochastic_vertex!(bn, :A, Normal(0, 1), false) + add_stochastic_vertex!(bn, :B, Normal(0, 1), false) + add_stochastic_vertex!(bn, :C, Normal(0, 1), false) + add_edge!(bn, :A, :B) + add_edge!(bn, :B, :C) + @testset "Corner Case: X or Y in Z" begin + # Test case where X is in Z + @test is_conditionally_independent(bn, :A, :C, [:A]) # A ⊥ C | A + # Test case where Y is in Z + @test is_conditionally_independent(bn, :A, :C, [:C]) # A ⊥ C | C + # Test case where both X and Y are in Z + @test is_conditionally_independent(bn, :A, :C, [:A, :C]) # A ⊥ C | A, C + end + end + + @testset "Complex Structure" begin + bn = BayesianNetwork{Symbol}() + + for v in [:A, :B, :C, :D, :E] + add_stochastic_vertex!(bn, v, Normal(), false) + end + + # Create structure: + # A → B → D + # ↓ ↑ + # C → E + add_edge!(bn, :A, :B) + add_edge!(bn, :B, :C) + add_edge!(bn, :B, :D) + add_edge!(bn, :C, :E) + add_edge!(bn, :E, :D) + + @test is_conditionally_independent(bn, :A, :E, [:B, :C]) + @test !is_conditionally_independent(bn, :A, :E, Symbol[]) + end + + @testset "Using Observed Variables" begin + bn = BayesianNetwork{Symbol}() + + add_stochastic_vertex!(bn, :A, Normal(), false) + add_stochastic_vertex!(bn, :B, Normal(), true) # B is observed + add_stochastic_vertex!(bn, :C, Normal(), false) + + add_edge!(bn, :A, :B) + add_edge!(bn, :B, :C) + + @test is_conditionally_independent(bn, :A, :C) + + bn_decond = decondition(bn) + @test !is_conditionally_independent(bn_decond, :A, :C) + end + + @testset "Error Handling" begin + bn = BayesianNetwork{Symbol}() + + add_stochastic_vertex!(bn, :A, Normal(), false) + add_stochastic_vertex!(bn, :B, Normal(), false) + + @test_throws KeyError is_conditionally_independent(bn, :A, :NonExistent) + @test_throws KeyError is_conditionally_independent(bn, :NonExistent, :B) + @test_throws KeyError is_conditionally_independent(bn, :A, :B, [:NonExistent]) + end + end + + @testset "Variable Elimination Tests" begin + println("\nTesting Variable Elimination") + + @testset "Simple Chain Network (Z → X → Y)" begin + # Create a simple chain network: Z → X → Y + bn = BayesianNetwork{Symbol}() + + # Add vertices with specific distributions + println("Adding vertices...") + add_stochastic_vertex!(bn, :Z, Categorical([0.7, 0.3]), false) # P(Z) + add_stochastic_vertex!(bn, :X, Normal(0, 1), false) # P(X|Z) + add_stochastic_vertex!(bn, :Y, Normal(1, 2), false) # P(Y|X) + + # Add edges + println("Adding edges...") + add_edge!(bn, :Z, :X) + add_edge!(bn, :X, :Y) + + # Test case 1: P(X | Y=1.5) + println("\nTest case 1: P(X | Y=1.5)") + evidence1 = Dict(:Y => 1.5) + query1 = :X + result1 = variable_elimination(bn, query1, evidence1) + @test result1 isa Number + @test result1 >= 0 + println("P(X | Y=1.5) = ", result1) + + # Test case 2: P(X | Z=1) + println("\nTest case 2: P(X | Z=1)") + evidence2 = Dict(:Z => 1) + query2 = :X + result2 = variable_elimination(bn, query2, evidence2) + @test result2 isa Number + @test result2 >= 0 + println("P(X | Z=1) = ", result2) + + # Test case 3: P(Y | Z=1) + println("\nTest case 3: P(Y | Z=1)") + evidence3 = Dict(:Z => 1) + query3 = :Y + result3 = variable_elimination(bn, query3, evidence3) + @test result3 isa Number + @test result3 >= 0 + println("P(Y | Z=1) = ", result3) + end + end + + @testset "Variable Elimination Tests" begin + println("\nTesting Variable Elimination") + + @testset "Simple Chain Network (Z → X → Y)" begin + # Create a simple chain network: Z → X → Y + bn = BayesianNetwork{Symbol}() + + # Add vertices with specific distributions + println("Adding vertices...") + add_stochastic_vertex!(bn, :Z, Categorical([0.7, 0.3]), false) # P(Z) + add_stochastic_vertex!(bn, :X, Normal(0, 1), false) # P(X|Z) + add_stochastic_vertex!(bn, :Y, Normal(1, 2), false) # P(Y|X) + + # Add edges + println("Adding edges...") + add_edge!(bn, :Z, :X) + add_edge!(bn, :X, :Y) + + # Test case 1: P(X | Y=1.5) + println("\nTest case 1: P(X | Y=1.5)") + evidence1 = Dict(:Y => 1.5) + query1 = :X + result1 = variable_elimination(bn, query1, evidence1) + @test result1 isa Number + @test result1 >= 0 + println("P(X | Y=1.5) = ", result1) + + # Test case 2: P(X | Z=1) + println("\nTest case 2: P(X | Z=1)") + evidence2 = Dict(:Z => 1) + query2 = :X + result2 = variable_elimination(bn, query2, evidence2) + @test result2 isa Number + @test result2 >= 0 + println("P(X | Z=1) = ", result2) + + # Test case 3: P(Y | Z=1) + println("\nTest case 3: P(Y | Z=1)") + evidence3 = Dict(:Z => 1) + query3 = :Y + result3 = variable_elimination(bn, query3, evidence3) + @test result3 isa Number + @test result3 >= 0 + println("P(Y | Z=1) = ", result3) + end + + @testset "Mixed Network (Discrete and Continuous)" begin + # Create a more complex network with both discrete and continuous variables + bn = BayesianNetwork{Symbol}() - @testset "Bayes Ball" begin end -end + # Add vertices + println("\nAdding vertices for mixed network...") + add_stochastic_vertex!(bn, :A, Categorical([0.4, 0.6]), false) # Discrete + add_stochastic_vertex!(bn, :B, Normal(0, 1), false) # Continuous + add_stochastic_vertex!(bn, :C, Categorical([0.3, 0.7]), false) # Discrete + add_stochastic_vertex!(bn, :D, Normal(1, 2), false) # Continuous + + # Add edges: A → B → D ← C + println("Adding edges...") + add_edge!(bn, :A, :B) + add_edge!(bn, :B, :D) + add_edge!(bn, :C, :D) + + # Test case 1: P(B | D=1.0) + println("\nTest case 1: P(B | D=1.0)") + evidence1 = Dict(:D => 1.0) + query1 = :B + result1 = variable_elimination(bn, query1, evidence1) + @test result1 isa Number + @test result1 >= 0 + println("P(B | D=1.0) = ", result1) + + # Test case 2: P(D | A=1, C=1) + println("\nTest case 2: P(D | A=1, C=1)") + evidence2 = Dict(:A => 1, :C => 1) + query2 = :D + result2 = variable_elimination(bn, query2, evidence2) + @test result2 isa Number + @test result2 >= 0 + println("P(D | A=1, C=1) = ", result2) + end + + @testset "Special Cases" begin + bn = BayesianNetwork{Symbol}() + + # Single node case + add_stochastic_vertex!(bn, :X, Normal(0, 1), false) + result = variable_elimination(bn, :X, Dict{Symbol,Any}()) + @test result isa Number + @test result >= 0 + + # No evidence case + add_stochastic_vertex!(bn, :Y, Normal(1, 2), false) + add_edge!(bn, :X, :Y) + result = variable_elimination(bn, :Y, Dict{Symbol,Any}()) + @test result isa Number + @test result >= 0 + end + end +end \ No newline at end of file From 300209ab5c888e40b64cc10aa3399fc811b27b0e Mon Sep 17 00:00:00 2001 From: Naseweisssss Date: Mon, 18 Nov 2024 14:39:28 +0000 Subject: [PATCH 2/3] keep consistent --- .../ProbabilisticGraphicalModels.jl | 12 + .../ProbabilisticGraphicalModels/bayesnet.jl | 109 +------- .../ProbabilisticGraphicalModels/bayesnet.jl | 233 +----------------- 3 files changed, 24 insertions(+), 330 deletions(-) diff --git a/src/experimental/ProbabilisticGraphicalModels/ProbabilisticGraphicalModels.jl b/src/experimental/ProbabilisticGraphicalModels/ProbabilisticGraphicalModels.jl index 4785bf151..17475102e 100644 --- a/src/experimental/ProbabilisticGraphicalModels/ProbabilisticGraphicalModels.jl +++ b/src/experimental/ProbabilisticGraphicalModels/ProbabilisticGraphicalModels.jl @@ -6,4 +6,16 @@ using Distributions include("bayesnet.jl") +export BayesianNetwork, + Factor, + create_factor, + multiply_factors, + marginalize, + add_stochastic_vertex!, + add_deterministic_vertex!, + add_edge!, + condition, + decondition, + variable_elimination + end diff --git a/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl b/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl index 8f237b389..b5a8bc1a9 100644 --- a/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl +++ b/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl @@ -4,7 +4,7 @@ A structure representing a Bayesian Network. """ struct BayesianNetwork{V,T,F} - graph::SimpleDiGraph{T} + graph::SimpleGraph{T} "names of the variables in the network" names::Vector{V} "mapping from variable names to ids" @@ -25,7 +25,7 @@ end function BayesianNetwork{V}() where {V} return BayesianNetwork( - SimpleDiGraph{Int}(), # by default, vertex ids are integers + SimpleGraph{Int}(), # by default, vertex ids are integers V[], Dict{V,Int}(), Dict{V,Any}(), @@ -164,36 +164,13 @@ Perform ancestral sampling on a Bayesian network to generate one sample from the Ancestral sampling works by: 1. Finding a topological ordering of the nodes 2. Sampling from each node in order, using the already-sampled parent values for conditional distributions - -### Return Value -The function returns a `Dict{V, Any}` where: -- Each key is a variable name (of type `V`) in the Bayesian Network. -- Each value is the sampled value for that variable, which can be of any type (`Any`). - -This dictionary represents a single sample from the joint distribution of the Bayesian Network, capturing the dependencies and conditional relationships defined in the network structure. - """ function ancestral_sampling(bn::BayesianNetwork{V}) where {V} - ordered_vertices = Graphs.topological_sort_by_dfs(bn.graph) + ordered_vertices = Graphs.topological_sort(bn.graph) + samples = Dict{V,Any}() - for vertex_id in ordered_vertices - vertex_name = bn.names[vertex_id] - if bn.is_observed[vertex_id] - samples[vertex_name] = bn.values[vertex_name] - continue - end - if bn.is_stochastic[vertex_id] - dist_idx = findfirst(id -> id == vertex_id, bn.stochastic_ids) - samples[vertex_name] = rand(bn.distributions[dist_idx]) - else - # deterministic node - parent_ids = Graphs.inneighbors(bn.graph, vertex_id) - parent_values = [samples[bn.names[pid]] for pid in parent_ids] - func_idx = findfirst(id -> id == vertex_id, bn.deterministic_ids) - samples[vertex_name] = bn.deterministic_functions[func_idx](parent_values...) - end - end + # TODO: Implement sampling logic return samples end @@ -207,82 +184,13 @@ If Z is provided, the conditioning information in `bn` will be ignored. function is_conditionally_independent end function is_conditionally_independent(bn::BayesianNetwork{V}, X::V, Y::V) where {V} - # Use currently observed variables as Z - Z = V[v for (v, is_obs) in zip(bn.names, bn.is_observed) if is_obs] - return is_conditionally_independent(bn, X, Y, Z) + # TODO: Implement end function is_conditionally_independent( bn::BayesianNetwork{V}, X::V, Y::V, Z::Vector{V} ) where {V} - println("debugging: X: $X, Y: $Y, Z: $Z") - if X in Z || Y in Z - return true - end - - # Get vertex IDs - x_id = bn.names_to_ids[X] - y_id = bn.names_to_ids[Y] - z_ids = Set([bn.names_to_ids[z] for z in Z]) - - # Track visited nodes and their states - n_vertices = nv(bn.graph) - visited = falses(n_vertices) - - # Queue entries are (node_id, from_parent) - queue = Tuple{Int,Bool}[] - - # Start from X - push!(queue, (x_id, true)) # As if coming from parent - push!(queue, (x_id, false)) # As if coming from child - - while !isempty(queue) - current_id, from_parent = popfirst!(queue) - - if visited[current_id] - continue - end - visited[current_id] = true - - # If we reached Y, path is active - if current_id == y_id - return false - end - - is_conditioned = current_id in z_ids - parents = inneighbors(bn.graph, current_id) - children = outneighbors(bn.graph, current_id) - - # Case 1: Node is not conditioned - if !is_conditioned - # Can go to children if coming from parent or at start node - if from_parent || current_id == x_id - for child in children - push!(queue, (child, true)) - end - end - - # Can go to parents if coming from child or at start node - if !from_parent || current_id == x_id - for parent in parents - push!(queue, (parent, false)) - end - end - end - - # Case 2: Node is conditioned or has conditioned descendants - if is_conditioned - # If this is a collider or descendant of collider - if length(parents) > 1 || !isempty(children) - # Can go to parents regardless of direction - for parent in parents - push!(queue, (parent, false)) - end - end - end - end - - return true + # TODO: Implement end using LinearAlgebra @@ -350,7 +258,6 @@ function marginalize(factor::Factor, var::Symbol) return Factor(new_vars, factor.distribution, new_parents) end - """ variable_elimination(bn::BayesianNetwork, query::Symbol, evidence::Dict{Symbol,Any}) @@ -456,4 +363,4 @@ function variable_elimination( # Convert evidence to Dict{Symbol,Float64} evidence_float = Dict{Symbol,Float64}(k => Float64(v) for (k, v) in evidence) return variable_elimination(bn, query, evidence_float) -end \ No newline at end of file +end diff --git a/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl b/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl index 395121a7e..6c5063869 100644 --- a/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl +++ b/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl @@ -11,6 +11,7 @@ using JuliaBUGS.ProbabilisticGraphicalModels: ancestral_sampling, is_conditionally_independent, variable_elimination + @testset "BayesianNetwork" begin @testset "Adding vertices" begin bn = BayesianNetwork{Symbol}() @@ -98,235 +99,9 @@ using JuliaBUGS.ProbabilisticGraphicalModels: @test bn_cond2.values[:B] == 2.0 end - @testset "Simple ancestral sampling" begin - bn = BayesianNetwork{Symbol}() - # Add stochastic vertices - add_stochastic_vertex!(bn, :A, Normal(0, 1), false) - add_stochastic_vertex!(bn, :B, Normal(1, 2), false) - # Add deterministic vertex C = A + B - add_deterministic_vertex!(bn, :C, (a, b) -> a + b) - add_edge!(bn, :A, :C) - add_edge!(bn, :B, :C) - samples = ancestral_sampling(bn) - @test haskey(samples, :A) - @test haskey(samples, :B) - @test haskey(samples, :C) - @test samples[:A] isa Number - @test samples[:B] isa Number - @test samples[:C] ≈ samples[:A] + samples[:B] - end - - @testset "Complex ancestral sampling" begin - bn = BayesianNetwork{Symbol}() - add_stochastic_vertex!(bn, :μ, Normal(0, 2), false) - add_stochastic_vertex!(bn, :σ, LogNormal(0, 0.5), false) - add_stochastic_vertex!(bn, :X, Normal(0, 1), false) - add_stochastic_vertex!(bn, :Y, Normal(0, 1), false) - add_deterministic_vertex!(bn, :X_scaled, (μ, σ, x) -> x * σ + μ) - add_deterministic_vertex!(bn, :Y_scaled, (μ, σ, y) -> y * σ + μ) - add_deterministic_vertex!(bn, :Sum, (x, y) -> x + y) - add_deterministic_vertex!(bn, :Product, (x, y) -> x * y) - add_deterministic_vertex!(bn, :N, () -> 2.0) - add_deterministic_vertex!(bn, :Mean, (s, n) -> s / n) - add_edge!(bn, :μ, :X_scaled) - add_edge!(bn, :σ, :X_scaled) - add_edge!(bn, :X, :X_scaled) - add_edge!(bn, :μ, :Y_scaled) - add_edge!(bn, :σ, :Y_scaled) - add_edge!(bn, :Y, :Y_scaled) - add_edge!(bn, :X_scaled, :Sum) - add_edge!(bn, :Y_scaled, :Sum) - add_edge!(bn, :X_scaled, :Product) - add_edge!(bn, :Y_scaled, :Product) - add_edge!(bn, :Sum, :Mean) - add_edge!(bn, :N, :Mean) - samples = ancestral_sampling(bn) - - @test all( - haskey(samples, k) for - k in [:μ, :σ, :X, :Y, :X_scaled, :Y_scaled, :Sum, :Product, :Mean, :N] - ) - - @test all(samples[k] isa Number for k in keys(samples)) - @test samples[:X_scaled] ≈ samples[:X] * samples[:σ] + samples[:μ] - @test samples[:Y_scaled] ≈ samples[:Y] * samples[:σ] + samples[:μ] - @test samples[:Sum] ≈ samples[:X_scaled] + samples[:Y_scaled] - @test samples[:Product] ≈ samples[:X_scaled] * samples[:Y_scaled] - @test samples[:Mean] ≈ samples[:Sum] / samples[:N] - @test samples[:N] ≈ 2.0 - @test samples[:σ] > 0 - # Multiple samples test - n_samples = 1000 - means = zeros(n_samples) - for i in 1:n_samples - samples = ancestral_sampling(bn) - means[i] = samples[:Mean] - end - - @test mean(means) ≈ 0 atol = 0.5 - @test std(means) > 0 - end - - @testset "Bayes Ball" begin - @testset "Chain Structure (A → B → C)" begin - bn = BayesianNetwork{Symbol}() - - add_stochastic_vertex!(bn, :A, Normal(), false) - add_stochastic_vertex!(bn, :B, Normal(), false) - add_stochastic_vertex!(bn, :C, Normal(), false) - - add_edge!(bn, :A, :B) - add_edge!(bn, :B, :C) - - @test is_conditionally_independent(bn, :A, :C, [:B]) - @test !is_conditionally_independent(bn, :A, :C, Symbol[]) - end - - @testset "Fork Structure (A ← B → C)" begin - println("\nTesting Fork Structure") - bn = BayesianNetwork{Symbol}() - - add_stochastic_vertex!(bn, :A, Normal(), false) - add_stochastic_vertex!(bn, :B, Normal(), false) - add_stochastic_vertex!(bn, :C, Normal(), false) - - add_edge!(bn, :B, :A) - add_edge!(bn, :B, :C) - - println("Graph structure:") - println("Edges: ", collect(edges(bn.graph))) - - result = is_conditionally_independent(bn, :A, :C, Symbol[]) - println("Result for A ⊥ C | ∅: $result") - end - - @testset "Collider Structure (A → B ← C)" begin - bn = BayesianNetwork{Symbol}() - - add_stochastic_vertex!(bn, :A, Normal(), false) - add_stochastic_vertex!(bn, :B, Normal(), false) - add_stochastic_vertex!(bn, :C, Normal(), false) - - add_edge!(bn, :A, :B) - add_edge!(bn, :C, :B) - - @test is_conditionally_independent(bn, :A, :C, Symbol[]) - @test !is_conditionally_independent(bn, :A, :C, [:B]) - end - - @testset "Bayes Ball Algorithm Tests" begin - # Create a simple network: A → B → C - bn = BayesianNetwork{Symbol}() - add_stochastic_vertex!(bn, :A, Normal(0, 1), false) - add_stochastic_vertex!(bn, :B, Normal(0, 1), false) - add_stochastic_vertex!(bn, :C, Normal(0, 1), false) - add_edge!(bn, :A, :B) - add_edge!(bn, :B, :C) - @testset "Corner Case: X or Y in Z" begin - # Test case where X is in Z - @test is_conditionally_independent(bn, :A, :C, [:A]) # A ⊥ C | A - # Test case where Y is in Z - @test is_conditionally_independent(bn, :A, :C, [:C]) # A ⊥ C | C - # Test case where both X and Y are in Z - @test is_conditionally_independent(bn, :A, :C, [:A, :C]) # A ⊥ C | A, C - end - end - - @testset "Complex Structure" begin - bn = BayesianNetwork{Symbol}() - - for v in [:A, :B, :C, :D, :E] - add_stochastic_vertex!(bn, v, Normal(), false) - end - - # Create structure: - # A → B → D - # ↓ ↑ - # C → E - add_edge!(bn, :A, :B) - add_edge!(bn, :B, :C) - add_edge!(bn, :B, :D) - add_edge!(bn, :C, :E) - add_edge!(bn, :E, :D) - - @test is_conditionally_independent(bn, :A, :E, [:B, :C]) - @test !is_conditionally_independent(bn, :A, :E, Symbol[]) - end - - @testset "Using Observed Variables" begin - bn = BayesianNetwork{Symbol}() - - add_stochastic_vertex!(bn, :A, Normal(), false) - add_stochastic_vertex!(bn, :B, Normal(), true) # B is observed - add_stochastic_vertex!(bn, :C, Normal(), false) - - add_edge!(bn, :A, :B) - add_edge!(bn, :B, :C) - - @test is_conditionally_independent(bn, :A, :C) - - bn_decond = decondition(bn) - @test !is_conditionally_independent(bn_decond, :A, :C) - end - - @testset "Error Handling" begin - bn = BayesianNetwork{Symbol}() - - add_stochastic_vertex!(bn, :A, Normal(), false) - add_stochastic_vertex!(bn, :B, Normal(), false) - - @test_throws KeyError is_conditionally_independent(bn, :A, :NonExistent) - @test_throws KeyError is_conditionally_independent(bn, :NonExistent, :B) - @test_throws KeyError is_conditionally_independent(bn, :A, :B, [:NonExistent]) - end - end - - @testset "Variable Elimination Tests" begin - println("\nTesting Variable Elimination") - - @testset "Simple Chain Network (Z → X → Y)" begin - # Create a simple chain network: Z → X → Y - bn = BayesianNetwork{Symbol}() - - # Add vertices with specific distributions - println("Adding vertices...") - add_stochastic_vertex!(bn, :Z, Categorical([0.7, 0.3]), false) # P(Z) - add_stochastic_vertex!(bn, :X, Normal(0, 1), false) # P(X|Z) - add_stochastic_vertex!(bn, :Y, Normal(1, 2), false) # P(Y|X) - - # Add edges - println("Adding edges...") - add_edge!(bn, :Z, :X) - add_edge!(bn, :X, :Y) - - # Test case 1: P(X | Y=1.5) - println("\nTest case 1: P(X | Y=1.5)") - evidence1 = Dict(:Y => 1.5) - query1 = :X - result1 = variable_elimination(bn, query1, evidence1) - @test result1 isa Number - @test result1 >= 0 - println("P(X | Y=1.5) = ", result1) - - # Test case 2: P(X | Z=1) - println("\nTest case 2: P(X | Z=1)") - evidence2 = Dict(:Z => 1) - query2 = :X - result2 = variable_elimination(bn, query2, evidence2) - @test result2 isa Number - @test result2 >= 0 - println("P(X | Z=1) = ", result2) + @testset "Simple ancestral sampling" begin end - # Test case 3: P(Y | Z=1) - println("\nTest case 3: P(Y | Z=1)") - evidence3 = Dict(:Z => 1) - query3 = :Y - result3 = variable_elimination(bn, query3, evidence3) - @test result3 isa Number - @test result3 >= 0 - println("P(Y | Z=1) = ", result3) - end - end + @testset "Bayes Ball" begin end @testset "Variable Elimination Tests" begin println("\nTesting Variable Elimination") @@ -427,4 +202,4 @@ using JuliaBUGS.ProbabilisticGraphicalModels: @test result >= 0 end end -end \ No newline at end of file +end From c8c13cec6b2679d909f606e05ee16f394fee7536 Mon Sep 17 00:00:00 2001 From: Naseweisssss Date: Wed, 20 Nov 2024 22:20:42 +0000 Subject: [PATCH 3/3] temp commit for toy example --- .../ProbabilisticGraphicalModels/bayesnet.jl | 16 +++- .../ProbabilisticGraphicalModels/bayesnet.jl | 86 +++++++++++++++++++ 2 files changed, 100 insertions(+), 2 deletions(-) diff --git a/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl b/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl index b5a8bc1a9..33654322b 100644 --- a/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl +++ b/src/experimental/ProbabilisticGraphicalModels/bayesnet.jl @@ -360,7 +360,19 @@ end function variable_elimination( bn::BayesianNetwork{Symbol,Int,Any}, query::Symbol, evidence::Dict{Symbol,<:Any} ) - # Convert evidence to Dict{Symbol,Float64} - evidence_float = Dict{Symbol,Float64}(k => Float64(v) for (k, v) in evidence) + # Convert evidence to Dict{Symbol,Float64}, handling both continuous and discrete cases + evidence_float = Dict{Symbol,Float64}() + for (k, v) in evidence + node_id = bn.names_to_ids[k] + dist_idx = findfirst(id -> id == node_id, bn.stochastic_ids) + + if bn.distributions[dist_idx] isa Categorical + # For categorical variables, keep the original value (0-based indexing) + evidence_float[k] = Float64(v) + else + # For continuous variables, convert to Float64 + evidence_float[k] = Float64(v) + end + end return variable_elimination(bn, query, evidence_float) end diff --git a/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl b/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl index 6c5063869..11a9317c4 100644 --- a/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl +++ b/test/experimental/ProbabilisticGraphicalModels/bayesnet.jl @@ -201,5 +201,91 @@ using JuliaBUGS.ProbabilisticGraphicalModels: @test result isa Number @test result >= 0 end + + @testset "Marginalization with Mixed Variables" begin + bn = BayesianNetwork{Symbol}() + + # X1 ~ Normal(0, 1) + add_stochastic_vertex!(bn, :X1, Normal(0, 1), false) + + # X2 ~ Bernoulli(0.7) [using Categorical with 2 categories] + add_stochastic_vertex!(bn, :X2, Categorical([0.3, 0.7]), false) + + # X3 ~ Normal(μ = 2*X2, σ = 1) + add_stochastic_vertex!(bn, :X3, Normal(0, 1), false) + add_edge!(bn, :X2, :X3) + + @testset "Marginalizing over X2" begin + # P(X3 | X1=0) + result1 = variable_elimination(bn, :X3, Dict(:X1 => 0.0)) + @test result1 isa Number + @test result1 > 0 + + # P(X3) - no evidence + result2 = variable_elimination(bn, :X3, Dict{Symbol,Any}()) + @test result2 isa Number + @test result2 > 0 + end + + @testset "Marginalizing over continuous" begin + # P(X2 | X3=1.0) + result3 = variable_elimination(bn, :X2, Dict(:X3 => 1.0)) + @test result3 isa Number + @test 0 ≤ result3 ≤ 1 # Should be a probability + + # P(X2 | X1=0, X3=1.0) + result4 = variable_elimination(bn, :X2, Dict(:X1 => 0.0, :X3 => 1.0)) + @test result4 isa Number + @test 0 ≤ result4 ≤ 1 + end + end + + @testset "Variable Elimination - Marginalization Demo" begin + bn = BayesianNetwork{Symbol}() + + # Temperature (X1) ~ Normal(0, 1) + add_stochastic_vertex!(bn, :X1, Normal(0, 1), false) + + # Rain (X2) ~ Bernoulli(0.7) + add_stochastic_vertex!(bn, :X2, Categorical([0.3, 0.7]), false) + + # Umbrella Sales (X3) ~ Normal(μ(X2), 1) + # μ = 2 if no rain (X2=0), μ = 10 if rain (X2=1) + add_stochastic_vertex!(bn, :X3, Normal(0, 1), false) + + add_edge!(bn, :X2, :X3) + + @testset "P(X3|X1) - Marginalizing over discrete X2" begin + # P(X3|X1) = ∫ P(X3|X2)P(X2|X1) dX2 + # Since X1 ⊥⊥ X2, P(X2|X1) = P(X2) + # So this should still be a mixture of two Gaussians: + # 0.3 * Normal(2,1) + 0.7 * Normal(10,1) + result1 = variable_elimination(bn, :X3, Dict(:X1 => 0.0)) + @test result1 isa Number + @test result1 > 0 + println("P(X3|X1) = ", result1) + end + + @testset "Marginalizing over continuous variable (X1)" begin + # P(X2) after marginalizing over X1 should still be close to original prior + # because X1 and X2 are independent + result = variable_elimination(bn, :X2, Dict{Symbol,Any}()) + @test result isa Number + @test 0 ≤ result ≤ 1 + @test isapprox(result, 0.3, atol=0.1) # Should be close to P(X2=0)=0.3 + end + + @testset "Conditional probabilities" begin + # P(X3|X2=1) should be approximately Normal(10,1) + result_rain = variable_elimination(bn, :X3, Dict(:X2 => 1)) + @test result_rain isa Number + @test result_rain > 0 + + # P(X2|X3=10) should be high (likely raining given high sales) + result_high_sales = variable_elimination(bn, :X2, Dict(:X3 => 10.0)) + @test result_high_sales isa Number + println("P(X2|X3=10) = ", result_high_sales) # Should favor rain hypothesis + end + end end end