rb_prob

https://github.com/geeknees/rb_prob

It’s not ideal, but to quickly fix the ‘method missing’ issue, I’ve added the code below.

https://github.com/geeknees/rb_prob/blob/master/lib/prob.rb#L251

require 'rubygems'
require 'bundler/setup'

require 'prob'
include Probably

# Alarm example from "Artificial Intelligence - A Modern Approach" by Russel
# and Norvig Page 493 cc.
#
# Suppose you have a new fairly reliable burglar alarm at home but occasionally
# it responds to minor earthquakes. You also have two neighbors John and Mary,
# who have promised to call you at work when they hear the alarm. John always
# calls when he hears the alarm, but sometimes confuses the telephone ringing
# with the alarm and calls then, too. Mary, on the other hand, is too much in
# loud music and sometimes misses the alarm altogether.
#
# So the bayesian network will be:
#
#           B         E
#            \       /
#            _\|   |/_
#                A
#             /    \
#           |/_    _\|
#          J          M
#
#  with probabilities:
#  P(B) = 0.001
#  P(E) = 0.002
#
#  P(A| B=true, E=true)   = 0.95
#  P(A| B=true, E=false)  = 0.94
#  P(A| B=false, E=true)  = 0.29
#  P(A| B=false, E=false) = 0.001
#
#  P(J| A=true)  = 0.9
#  P(J| A=false) = 0.05
#
#  P(M| A=true)  = 0.7
#  P(M| A=false) = 0.01
#
#  where B = burglar, E = earthquake, A = alarm, J = John calls and
#  M = Mary calls
#
#  ----------------------------------------------------------------------------
#
#  Next we want to develop some 'equivalent' functions for querying that
#  network and do some benchmarks.
#

# first let's encode the probabilities from the network
# P(B)
PBurglary = choose(0.001,  :B, :notB )

# P(E)
PEarthquake = choose(0.002,  :E, :notE)

# P(A|B = b,E = e)
def p_alarm(b, e)
    pAlarmTable = {
        [:B, :E] => 0.95,
        [:B, :notE] => 0.94,
        [:notB, :E] => 0.29,
        [:notB, :notE] => 0.001
    }

    choose(pAlarmTable[[b, e]],  :A, :notA)
end

# P(J|A = a)
def p_john(a)
    choose( a == :A ? 0.9 : 0.05, :J, :notJ)
end

# P(M|A = a)
def p_mary(a)
    choose( a == :A ? 0.7 : 0.01, :M, :notM)
end

# computes the joint probability and transform result using block (if given)
# allowing to do some marginalization over one random variable by
# "leaving it out"
#
# for example:
# mk_joint_p {|b,e,a,j,m| [b,e,a]} will find P(b,e,a) = Sum(j,m) { P(b,e,a,j,m) }
#
def mk_joint_p(&blk)
    PBurglary.dep { |b|
        PEarthquake.dep {|e|
            p_alarm(b, e).dep {|a|
                p_john(a).dep { |j|
                    p_mary(a).dep {|m|
                        mk_const(if blk then blk.call([b,e,a,j,m])
                                else [b,e,a,j,m] end)
                    }
                }
            }
        }
    }
end

# compute (optionally conditional) joint probability of (free) random
# variables like mk_joint_p.
#
# To compute conditional probability set random variables to a known state.
# for example
# mk_joint_p2( {:john = :J, :mary = :M} )
# will compute
# P(B,E,A| J = true, M = true)
#
# or
# mk_joint_p2({:john = :J, :mary = :M}) {|b,e,a,j,m| b} will find
# P(B | J = true, M = true)
def mk_joint_p2( tsts = {}, &blk )
    PBurglary.dep { |b|
    condition(!tsts[:burglary] || tsts[:burglary] == b) {
        PEarthquake.dep {|e|
        condition(!tsts[:earthquake] || tsts[:earthquake] == e) {
            p_alarm(b,e).dep {|a|
            condition(!tsts[:alarm] || tsts[:alarm] == a) {
                p_john(a).dep {|j|
                condition(!tsts[:john] || tsts[:john] == j) {
                    p_mary(a).dep {|m|
                    condition(!tsts[:mary] || tsts[:mary] == m) {
                        mk_const(if blk then blk.call [b,e,a,j,m] else [b,e,a,j,m] end)
                    }}
                }}
            }}
        }}
    }}.normalize
end

# like mk_joint_p2, but using event_dep directly instead of mixing in
# condition-statements
def mk_joint_p3 (tsts = {}, &blk)
    tst_b = if_just tsts[:burglary]
    tst_e = if_just tsts[:earthquake]
    tst_a = if_just tsts[:alarm]
    tst_j = if_just tsts[:john]
    tst_m = if_just tsts[:mary]

    PBurglary.event_dep(tst_b) {|b|
        PEarthquake.event_dep(tst_e) {|e|
            p_alarm(b,e).event_dep(tst_a) {|a|
                p_john(a).event_dep(tst_j) {|j|
                    p_mary(a).event_dep(tst_m) {|m|
                        mk_const(if blk then blk.call [b,e,a,j,m] else [b,e,a,j,m] end)
                    }
                }
            }
        }
    }.normalize
end

# precompute joint probability to do bayesian inference using filter, map and
# query?
PJoint = mk_joint_p

puts 'P(B|M=true, J=true) :'
puts mk_joint_p3({:mary => :M, :john => :J}) {|b,e,a,j,m| b }

# puts "\njoint probability:"
# puts "=================="
# puts PJoint

# compute P(B | M=true, J=true, E=false, A=true) using all 3 different
# functions mk_joint_p, mk_joint_p2 and mk_joint_p3:
puts "\nP(B | M=true, J=true, E=false, A=true)"
puts "====================================="
puts mk_joint_p2({:mary => :M, :john => :J, :earthquake => :notE, :alarm => :A}) { |b,e,a,j,m| b }.query?(&just(:B))
puts mk_joint_p3({:mary => :M, :john => :J, :earthquake => :notE, :alarm => :A}) { |b,e,a,j,m| b }.probability(:B)
puts PJoint.filter {|b,e,a,j,m| e == :notE && j == :J && m == :M && a == :A }.query? {|b,e,a,j,m| b == :B }

# do some benchmarking:

require 'benchmark'

Benchmark.bmbm {|x|
    i = 1000
    x.report('joint probability:') {
        (1..i).each {
            mk_joint_p.filter {|b,e,a,j,m| e == :notE && j == :J && m == :M && a == :A }.query? {|b,e,a,j,m| b == :B }
        }
    }

    x.report('joint probability precomputed:') {
        (1..i).each {
            PJoint.filter {|b,e,a,j,m| e == :notE && j == :J && m == :M && a == :A}.query? {|b,e,a,j,m| b == :B}
        }
    }

    x.report('direkt:') {
        (1..i).each {
            mk_joint_p {|b,e,a,j,m|
                if e == :notE && j == :J && m == :M  && a == :A
                    [b,a]
                else
                    nil
                end
            }.query? {|b,a| b == :B}
        }
    }

    x.report('direkt with conditions:') {
        (1..i).each {
            mk_joint_p2({:mary => :M, :john => :J, :earthquake => :notE, :alarm => :A}) { |b,e,a,j,m| b }.query?(&just(:B))
        }
    }

    x.report('direkt with event condition:') {
        (1..i).each {
            mk_joint_p3({:mary => :M, :john => :J, :earthquake => :notE, :alarm => :A}) { |b,e,a,j,m| b }.query?(&just(:B))
        }
    }
}

P(B|M=true, J=true) :
B : 28.41718353643929 %
notB : 71.58281646356072 %

P(B | M=true, J=true, E=false, A=true)
=====================================
0.4847859721505931
0.4847859721505931
0.484785972150593
Rehearsal ------------------------------------------------------------------
joint probability:               0.113981   0.000574   0.114555 (  0.114670)
joint probability precomputed:   0.003608   0.000209   0.003817 (  0.003818)
direkt:                          0.048796   0.000252   0.049048 (  0.049137)
direkt with conditions:          0.015159   0.000420   0.015579 (  0.015591)
direkt with event condition:     0.014577   0.000215   0.014792 (  0.014786)
--------------------------------------------------------- total: 0.197791sec

                                     user     system      total        real
joint probability:               0.108139   0.000171   0.108310 (  0.108459)
joint probability precomputed:   0.003247   0.000047   0.003294 (  0.003271)
direkt:                          0.047736   0.000101   0.047837 (  0.047904)
direkt with conditions:          0.013683   0.000041   0.013724 (  0.013714)
direkt with event condition:     0.013739   0.000069   0.013808 (  0.013811)

[#<Benchmark::Tms:0x0000000123cacc30 @label="joint probability:", @real=0.10845900001004338, @cstime=0.0, @cutime=0.0, @stime=0.00017100000000000448, @utime=0.10813899999999999, @total=0.10830999999999999>, #<Benchmark::Tms:0x000000012337b648 @label="joint probability precomputed:", @real=0.0032709999941289425, @cstime=0.0, @cutime=0.0, @stime=4.700000000000537e-05, @utime=0.003247, @total=0.0032940000000000053>, #<Benchmark::Tms:0x00000001233bc850 @label="direkt:", @real=0.047903999977279454, @cstime=0.0, @cutime=0.0, @stime=0.00010099999999998999, @utime=0.047736, @total=0.04783699999999999>, #<Benchmark::Tms:0x0000000123caf1b0 @label="direkt with conditions:", @real=0.013714000000618398, @cstime=0.0, @cutime=0.0, @stime=4.099999999999937e-05, @utime=0.013683, @total=0.013724>, #<Benchmark::Tms:0x0000000123caf160 @label="direkt with event condition:", @real=0.013810999982524663, @cstime=0.0, @cutime=0.0, @stime=6.899999999999962e-05, @utime=0.013739000000000057, @total=0.013808000000000056>]