Modeling Dependence of Many Variables on a Single Variable
=================================
In this notebook we translate into Probabilistic Logic Programming the Figaro example from Section 6.1.2 of [A. Pfeffer, Pratical Probabilistic Programming](https://www.manning.com/books/practical-probabilistic-programming), pages 174-176.
Imagine you have found a coin of dubious origin: it might be biased or not so you don't know the probability that it'll lands heads on any given toss. The goal of this program is to estimate the bias of the coin and to predict consecutive coin tosses. This is how the Bayesian network looks like, where the root node is Bias, a real value variable that represent the probability that any given toss will be heads:
:- use_rendering(graphviz).
X = dot('digraph G {
Bias->Toss1;
Bias->Toss2;
Bias->"...";
Bias->TossN}').
Performing inference
--------------------
We show two example queries. The first is mc_lw_sample(next_toss(h),previous_tosses([h,t,h,h,h,h,h,t,h,t,h,h,h,t,h,h,h,h,h,h]),1000,P) and computes the probability P that the next toss will be heads, given the evidence expressed with the list into the predicate =previous_tosses/1=, sampling 1000 times the query. The second is mc_lw_expectation(bias(B),
previous_tosses([h,t,h,h,h,h,h,t,h,t,h,h,h,t,h,h,h,h,h,h]),
1000,B,E) and computes the average bias B of the coin given a list of previous tosses, represented by a list into the predicate =previous_tosses/1=, sampling 1000 times the query.
Code Description
----------------
The predicate =next_toss/1= has only one argument and works in two steps: evaluates the probability of a variable called =Bias= through the predicate =bias/1= which computes a beta distribution using beta(Var,Alpha,Beta) where =Var= follows a beta distribution with parameters =Alpha= and =Beta=. Then the predicate =flip/3= is used to computed the probability of the variable =T=, given a certain bias.
=previous_tosses/1= is used to assert the evidence. It computes the length of the list in input then, for each element of this list (which is a result of a previous toss, so head or tail), computes the bias using =bias/1= and the probability using =flip/3=.
:- use_module(library(mcintyre)).
:- mc.
:- begin_lpad.
flip(h,P,N):P; flip(t,P,N):1-P.
bias(X): beta(X,2,5).
next_toss(T):-
bias(Bias),
flip(T,Bias,0).
previous_tosses(Tosses):-
length(Tosses,NumTosses),
tosses(Tosses,NumTosses).
tosses([],0).
tosses([H|T],NT):-
bias(Bias),
flip(H,Bias,NT),
NT1 is NT-1,
tosses(T,NT1).
:- end_lpad.
Results
----------------------
With this two queries we can compute, given a list of previous tosses, the probability of heads on next toss
mc_lw_sample(next_toss(h),
previous_tosses([h,t,h,h,h,h,h,t,h,t,h,h,h,t,h,h,h,h,h,h]),
1000,P).
and the average bias.
mc_lw_expectation(bias(B),
previous_tosses([h,t,h,h,h,h,h,t,h,t,h,h,h,t,h,h,h,h,h,h]),
1000,B,E).