1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
|
/*
Mixture of two Gaussians. A biased coin is thrown, if it lands heads X in mix(X)
is sampled from a Gaussian with mean 0 and variance 1. if it lands tails X is
sampled from a Gaussian with mean 5 and variance 2.
The example illustrates the use of continuous random variables and
the use of sampling, including
rejection sampling and Metropolis/Hastings. Moreover the example
illustrates the use of the predicate histogram/3 for graphing the
probability density function of continuous random variables.
*/
:- use_module(library(mcintyre)).
:- if(current_predicate(use_rendering/1)).
:- use_rendering(c3).
:- endif.
:- mc.
:- begin_lpad.
heads:0.6;tails:0.4.
% a coin is thrown. The coin is biased: with probability 0.6 it lands heads,
% with probability 0.4 it lands tails
g(X): gaussian(X,0, 1).
% X in g(X) follows a Gaussian distribution with mean 0 and variance 1
h(X): gaussian(X,5, 2).
% X in h(X) follows a Gaussian distribution with mean 5 and variance 2
mix(X) :- heads, g(X).
% if the coin lands heads, X in mix(X) is given by g(X)
mix(X) :- tails, h(X).
% if the coin lands tails, X in mix(X) is given by h(X)
:- end_lpad.
hist_uncond(Samples,NBins,Chart):-
mc_sample_arg(mix(X),Samples,X,L0),
histogram(L0,Chart,[nbins(NBins)]).
% take SAmples samples of X in mix(X) and draw a histogram with NBins bins representing
% the probability density of X
hist_rej_heads(Samples,NBins,Chart):-
mc_rejection_sample_arg(mix(X),heads,Samples,X,L0),
histogram(L0,Chart,[nbins(NBins)]).
% take Samples samples of X in mix(X) given that heads was true using
% rejection sampling and draw an
% histogram with NBins bins representing the probability density of X
hist_mh_heads(Samples,Lag,NBins,Chart):-
mc_mh_sample_arg(mix(X),heads,Samples,X,L0,[lag(Lag)]),
histogram(L0,Chart,[nbins(NBins)]).
% take Samples samples of X in mix(X) given that heads was true using
% Metropolis-Hastings and draw an
% histogram with NBins bins representing the probability density of X
hist_rej_dis(Samples,NBins,Chart):-
mc_rejection_sample_arg(mix(X),(mix(Y),Y>2),Samples,X,L0),
histogram(L0,Chart,[nbins(NBins)]).
% take Samples samples of X in mix(X) given that X>2 was true using
% rejection sampling and draw an
% histogram with NBins bins representing the probability density of X
hist_mh_dis(Samples,Lag,NBins,Chart):-
mc_mh_sample_arg(mix(X),(mix(Y),Y>2),Samples,X,L0,[lag(Lag)]),
histogram(L0,Chart,[nbins(NBins)]).
% take Samples samples of X in mix(X) given that X>2 was true using
% Metropolis-Hastings and draw an
% histogram with NBins bins representing the probability density of X
/** <examples>
?- hist_uncond(1000,40,G).
% take 1000 samples of X in mix(X) and draw a histogram with 40 bins representing
% the probability density of X
?- mc_sample_arg(mix(X),1000,X,L),histogram(L,Chart,[nbins(40)]).
% take 1000 samples of X in mix(X) and draw a histogram with 40 bins representing
% the probability density of X
?- mc_expectation(mix(X),1000,X,E).
% E=2.017964749114414
?- hist_rej_heads(1000,40,G).
% take 1000 samples of X in mix(X) given that heads was true using
% rejection sampling and draw an
% histogram with 40 bins representing the probability density of X
?- hist_mh_heads(1000,2,40,G).
% take 1000 samples of X in mix(X) given that heads was true using
% Metropolis-Hastings and draw an
% histogram with 40 bins representing the probability density of X
?- mc_mh_expectation(mix(X),heads,1000,X,E,[lag(2)]).
% E=-0.018433307290594284
?- hist_rej_dis(1000,40,G).
% take 1000 samples of X in mix(X) given that X>2 was true using
% rejection sampling and draw an
% histogram with 40 bins representing the probability density of X
?- hist_mh_dis(1000,2,40,G).
% take 1000 samples of X in mix(X) given that X>2 was true using
% Metropolis-Hastings and draw an
% histogram with 40 bins representing the probability density of X
?- mc_mh_expectation(mix(X),(mix(Y),Y>2),1000,X,E,[lag(2)]).
*/
|
