1/*
    2Factory producing widgets.
    3Consider a factory with two machines a and b. Each machine produces a widget
    4with a continuous feature. A widget is produced by machine a with probability
    50.7 and by machine b with probability b.
    6If the widget is produced by machine a, the feature is distributed as a
    7Gaussian with mean 2.0 and variance 1.0.
    8If the widget is produced by machine b, the feature is distributed as a
    9Gaussian with mean 3.0 and variance 1.0.
   10The widget then is processed by a third machine that adds a random quantity to
   11the feature distributed as a Gaussian with mean 0.5 and variance 0.1.
   12What is the distribution of the feature?
   13What is the distribution of the feature given that the widget was procuded
   14by machine a?
   15What is the distribution of the feature given that the third machine added a
   16quantity greater than 0.2?
   17What is the distribution of the feature given that the third machine added
   18a quantity of 2.0?
   19Adapted from
   20Islam, Muhammad Asiful, C. R. Ramakrishnan, and I. V. Ramakrishnan. 
   21"Inference in probabilistic logic programs with continuous random variables." 
   22Theory and Practice of Logic Programming 12.4-5 (2012): 505-523.
   23http://arxiv.org/pdf/1112.2681v3.pdf
   24*/
   25:- use_module(library(mcintyre)).   26:- use_module(library(clpr)).   27
   28:- if(current_predicate(use_rendering/1)).   29:- use_rendering(c3).   30:- endif.   31:- mc.   32:- begin_lpad.   33widget(X) :- msw(m, M),
   34msw(st(M), Z),
   35msw(pt, Y),
   36{X = Y + Z}.
   37% Ranges of RVs
   38values(m, [a,b]).
   39values(st(_), real).
   40values(pt, real).
   41% PDFs and PMFs:
   42:- set_sw(m, [0.3, 0.7]),
   43set_sw(st(a), norm(2.0, 1.0)),
   44set_sw(st(b), norm(3.0, 1.0)),
   45set_sw(pt, norm(0.5, 0.1)).   46
   47
   48:- end_lpad.   49
   50hist_uncond(Samples,NBins,Chart):-
   51  mc_sample_arg(widget(X),Samples,X,L0),
   52  histogram(L0,Chart,[nbins(NBins)]).
   53% What is the distribution of the feature?