1/*
    2Random arithmetic function from http://forestdb.org/models/arithmetic.html
    3The model generatively defines a random arithmetic function.
    4The problem is to predict the value returned by the function given one or
    5two couples of input-output, i.e., to compute a conditional probability.
    6Translated from the Church functional probabilistic programming language.
    7Sampling is necessary as queries have an infinite number of explanations.
    8Both rejection sampling and Metropolis/Hastings can be applied.
    9*/
   10:- use_module(library(mcintyre)).   11
   12:- if(current_predicate(use_rendering/1)).   13:- use_rendering(c3).   14:- endif.   15
   16:- mc.   17
   18:- begin_lpad.   19
   20eval(X,Y):-
   21  random_fn(X,0,F),
   22  Y is F.
   23
   24op(L,+):0.5;op(L,-):0.5.
   25
   26
   27
   28random_fn(X,L,F):-
   29  comb(L),
   30  random_fn(X,l(L),F1),
   31  random_fn(X,r(L),F2),
   32  op(L,Op),
   33  F=..[Op,F1,F2].
   34
   35random_fn(X,L,F):-
   36  \+ comb(L),
   37  base_random_fn(X,L,F).
   38
   39comb(_):0.3.
   40
   41base_random_fn(X,L,X):-
   42  identity(L).
   43
   44base_random_fn(_X,L,C):-
   45  \+ identity(L),
   46  random_const(L,C).
   47
   48identity(_):0.5.
   49
   50random_const(_,C):discrete(C,[0:0.1,1:0.1,2:0.1,3:0.1,4:0.1,
   51  5:0.1,6:0.1,7:0.1,8:0.1,9:0.1]).
   52
   53
   54
   55
   56
   57:- end_lpad.

?- mc_mh_sample(eval(2,4),eval(1,3),100,P, [mix(100),lag(3),successes(T),failures(F)]). % perform Metropolis Hastings sampling of eval(2,Y) given that % eval(1,3) is true (100 samples, 100 mixing samplesi, lag 3) % expected result % T = 17, % F = 83, % P = 0.17.

?- mc_mh_sample(eval(2,4),eval(1,3),100,P, [lag(3),successes(T),failures(F)]). % perform Metropolis Hastings sampling of eval(2,Y) given that % eval(1,3) is true % expected result % T = 17, % F = 83, % P = 0.17.

?- mc_mh_sample(eval(2,4),(eval(0,2),eval(1,3)),100,P, [lag(3),successes(T),failures(F)]). % perform Metropolis Hastings sampling of eval(2,Y) given that % eval(0,2) and eval(1,3) are true % expected result % T = 100, % F = 0, % P = 1.

?- mc_gibbs_sample(eval(2,4),eval(1,3),100,P, [mix(100),successes(T),failures(F)]). % perform Metropolis Hastings sampling of eval(2,Y) given that % eval(1,3) is true (100 samples, 100 mixing samples, lag 3) % expected result % T = 17, % F = 83, % P = 0.17.

?- mc_gibbs_sample(eval(2,4),eval(1,3),500,P,[]). % expected result % P = 0.17.

?- mc_gibbs_sample(eval(2,4),eval(1,3),100,P, [successes(T),failures(F),mix(100)]). % perform Gibbs sampling of eval(2,Y) given that % eval(1,3) is true % expected result % T = 17, % F = 83, % P = 0.17.

?- mc_gibbs_sample(eval(2,4),(eval(0,2),eval(1,3)),100,P, []). % perform Gibbs sampling of eval(2,Y) given that % eval(0,2) and eval(1,3) are true % expected result % P = 1.

?- mc_rejection_sample(eval(2,4),eval(1,3),100,P, [successes(T),failures(F)]). % perform rejection sampling of eval(2,4) given that eval(1,3) is true % expected result % T = 10, % F = 90, % P = 0.1.

?- mc_mh_sample_arg(eval(2,Y),(eval(0,2),eval(1,3)),100,Y,V, [mix(100),lag(3)]). % sample arg Y of eval(2,Y) given that % eval(0,2) and eval(1,3) are true % Sample using Metropolis Hastings % expected result % V = [[4]-100]. ?- mc_mh_sample_arg(eval(2,Y),(eval(0,2),eval(1,3)),100,Y,V, [mix(100),lag(3)]),argbar(V,C).

?- mc_mh_sample_arg(eval(2,Y),eval(1,3),100,Y,V, [mix(100),lag(3)]). % sample arg Y of eval(2,Y) given that % eval(1,3) is true % Sample using Metropolis Hastings % expected result % V = [[3]-52, [6]-20, [5]-16, [4]-12] ?- mc_mh_sample_arg(eval(2,Y),eval(1,3),100,Y,V, [mix(100),lag(3)]),argbar(V,C).

?- mc_rejection_sample_arg(eval(2,Y),eval(1,3),100,Y,V). % sample argument Y of eval(2,Y) given that % eval(1,3) is true % Sample using rejection sampling % expected result % V = [[3]-79, [4]-8, [6]-8, [2]-5]. ?- mc_rejection_sample_arg(eval(2,Y),eval(1,3),100,Y,V),argbar(V,C).

?- mc_expectation(eval(2,Y),100,Y,E). % what is the expected value of Y in eval(2,Y)? % expected result % E = 3.48

?- mc_mh_expectation(eval(2,Y),eval(1,3),100,Y,E, [mix(100),lag(3)]). % what is the expected value of Y in eval(2,Y) given that eval(1,3) is true? % expected result % E = 3.52 ?- mc_rejection_expectation(eval(2,Y),eval(1,3),100,Y,E). % what is the expected value of Y in eval(2,Y) given that eval(1,3) is true? % expected result % E = 3.06 */