1/*
    2Mixture of two Gaussians. A biased coin is thrown, if it lands heads X in mix(X)
    3is sampled from a Gaussian with mean 0 and variance 1. if it lands tails X is
    4sampled from a Gaussian with mean 5 and variance 2.
    5The example illustrates the use of continuous random variables and
    6the use of sampling, including
    7rejection sampling and Metropolis/Hastings. Moreover the example
    8illustrates the use of the predicate histogram/3 for graphing the
    9probability density function of continuous random variables.
   10*/
   11:- use_module(library(mcintyre)).   12
   13:- if(current_predicate(use_rendering/1)).   14:- use_rendering(c3).   15:- endif.   16:- mc.   17:- begin_lpad.   18
   19heads:0.6;tails:0.4. 
   20% a coin is thrown. The coin is biased: with probability 0.6 it lands heads,
   21% with probability 0.4 it lands tails
   22
   23g(X): gaussian(X,0, 1).
   24% X in g(X)  follows a Gaussian distribution with mean 0 and variance 1
   25h(X): gaussian(X,5, 2).
   26% X in h(X)  follows a Gaussian distribution with mean 5 and variance 2
   27
   28mix(X) :- heads, g(X).
   29% if the coin lands heads, X in mix(X) is given by g(X)
   30mix(X) :- tails, h(X).
   31% if the coin lands tails, X in mix(X) is given by h(X)
   32
   33:- end_lpad.   34
   35hist_uncond(Samples,NBins,Chart):-
   36  mc_sample_arg(mix(X),Samples,X,L0),
   37  histogram(L0,Chart,[nbins(NBins)]).
   38% take SAmples samples of X in mix(X) and draw a histogram with NBins bins representing 
   39% the probability density of X 
   40
   41hist_rej_heads(Samples,NBins,Chart):-
   42  mc_rejection_sample_arg(mix(X),heads,Samples,X,L0),
   43  histogram(L0,Chart,[nbins(NBins)]).
   44% take Samples samples of X in mix(X) given that heads was true using 
   45% rejection sampling and draw an
   46% histogram with NBins bins representing the probability density of X
   47
   48hist_mh_heads(Samples,Lag,NBins,Chart):-
   49  mc_mh_sample_arg(mix(X),heads,Samples,X,L0,[lag(Lag)]),
   50  histogram(L0,Chart,[nbins(NBins)]).
   51% take Samples samples of X in mix(X) given that heads was true using 
   52% Metropolis-Hastings and draw an
   53% histogram with NBins bins representing the probability density of X
   54
   55hist_rej_dis(Samples,NBins,Chart):-
   56  mc_rejection_sample_arg(mix(X),(mix(Y),Y>2),Samples,X,L0),
   57  histogram(L0,Chart,[nbins(NBins)]).
   58% take Samples samples of X in mix(X) given that X>2 was true using 
   59% rejection sampling and draw an
   60% histogram with NBins bins representing the probability density of X
   61
   62hist_mh_dis(Samples,Lag,NBins,Chart):-
   63  mc_mh_sample_arg(mix(X),(mix(Y),Y>2),Samples,X,L0,[lag(Lag)]),
   64  histogram(L0,Chart,[nbins(NBins)]).
   65% take Samples samples of X in mix(X) given that X>2 was true using 
   66% Metropolis-Hastings and draw an
   67% histogram with NBins bins representing the probability density of X