# Random arithmetic functions In this example we want to show how to perform conditional inference in an approximate way using sampling. In particular, we will show how to use rejection sampling and Metropolis-Hastings. This example generatively defines a random arithmetic function. The problem is to predict the value returned by the function given one or two couples of input-output, i.e., to compute a conditional probability. This program is translated from the example http://forestdb.org/models/arithmetic.html in the Church functional probabilistic programming language. Sampling is necessary as queries have an infinite number of explanations. ## Full program The full program of this example is
:- use_module(library(mcintyre)). :- if(current_predicate(use_rendering/1)). :- use_rendering(c3). :- endif. :- mc. :- begin_lpad. eval(X,Y):- random_fn(X,0,F), Y is F. op(+):0.5;op(-):0.5. random_fn(X,L,F):- comb(L), random_fn(X,l(L),F1), random_fn(X,r(L),F2), op(Op), F=..[Op,F1,F2]. random_fn(X,L,F):- \+ comb(L), base_random_fn(X,L,F). comb(_):0.3. base_random_fn(X,L,X):- identity(L). base_random_fn(_X,L,C):- \+ identity(L), random_const(L,C). identity(_):0.5. random_const(L,0):0.1;random_const(L,1):0.1;random_const(L,2):0.1; random_const(L,3):0.1;random_const(L,4):0.1;random_const(L,5):0.1; random_const(L,6):0.1;random_const(L,7):0.1;random_const(L,8):0.1; random_const(L,9):0.1. :- end_lpad.
We know that the random function return 3 for input 1 and we want to compute the probability that it returns 4 for input 2. We thus need to ask a conditional query and sampling is necessary as queries have an infinite number of explanations. The simplest approach is to use rejection sampling: you first query the evidence and, if the query is successful, query the goal in the same sample, otherwise the sample is discarded. This can be done with == mc_rejection_sample(:Query:atom,:Evidence:atom,+Samples:int, -Successes:int,-Failures:int,-Probability:float). == that takes =Samples= samples of =Query= given that =Evidence= is true. An example of use of the above predicate is
mc_rejection_sample(eval(2,4),eval(1,3),1000,T,F,P).
that perform rejection sampling of eval(2,4) given that eval(1,3) is true. Differently from exact inference, in approximate inference the evidence can be a conjunction of atoms, so if you also know that eval(0,2) is true, you can use
mc_rejection_sample(eval(2,4),(eval(0,2),eval(1,3)),1000,T,F,P).
and, as you can see, the query with more evidence is now almost certainly true. In Metropolis-Hastings MCMC, a Markov chain is produced using the algorithm of [Nampally, Ramakrishnan, 2014]: after a sample, a number of sampled probabilistic choices are deleted and the others are retained for the next sample. The sample is accepted with a probability of min{1,N0/N1} where N0 is the number of choices sampled in the previous sample and N1 is the number of choices sampled in the current sample. Metropolis-Hastings is usually much faster than rejection sampling because less samples are discarded. To use Metropolis-Hastings, the following predicate is available == mc_mh_sample(:Query:atom,:Evidence:atom,+Samples:int,+Lag:int, -Successes:int,-Failures:int,-Probability:float). == where =Lag= is the number of sampled choices to forget before taking a new sample. For example
mc_mh_sample(eval(2,4),eval(1,3),10000,1,T,F,P).
takes 10000 accepted samples and returns in =T= the number of samples where eval(2,4) is true, in =F= the number of samples where eval(2,4) is false and in =P= the estimated probability (=T/10000=). You can also compute conditional expectations with == mc_mh_expectation(:Query:atom,:Evidence:atom,+N:int,+Lag:int,?Arg:var,-Exp:float). == as in
mc_mh_expectation(eval(2,Y),eval(1,3),1000,1,Y,E).
that computes the expectation of argument =Y= of eval(2,Y) given that eval(1,3) is true by taking 1000 samples using Metropolis-Hastings MCMC.
-- Complete example: [arithm.pl](example/inference/arithm.pl) -- - References: - Nampally, Arun, and C. R. Ramakrishnan. _Adaptive MCMC-Based Inference in Probabilistic Logic Programs_. arXiv preprint arXiv:1403.6036 (2014). - http://forestdb.org/models/arithmetic.html
-- [Back to Tutorial](tutorial/tutorial.swinb)