1/*
    2One-dimensional  Kalman filter. Hidden Markov model with a real
    3value as state and a real value as output. The next state is given by
    4the current state plus Gaussian noise (mean 0 and variance 2 in this example)
    5and the output is given by the current state plus Gaussian noise (mean
    60 and variance 1 in this example). 
    7This example can be considered as modeling a random walk of a single continuous 
    8state variable with noisy observations. 
    9Given that at time 0 the value 2.5 was
   10observed, what is the distribution of the state at time 1 (filtering problem)?
   11The distribution of the state is plotted in the case of having (posterior) or 
   12not having the observation (prior).
   13Liklihood weighing is used to condition the distribution on evidence on
   14a continuous random variable (evidence with probability 0).
   15CLP(R) constraints allow both sampling and weighing samples with the same
   16program.
   17From
   18Islam, Muhammad Asiful, C. R. Ramakrishnan, and I. V. Ramakrishnan. 
   19"Inference in probabilistic logic programs with continuous random variables." 
   20Theory and Practice of Logic Programming 12.4-5 (2012): 505-523.
   21http://arxiv.org/pdf/1112.2681v3.pdf
   22Russell, S. and Norvig, P. 2010. Arficial Intelligence: A Modern Approach. 
   23Third Edition, Prentice Hall, Figure 15.10 page 587
   24PRISM syntax.
   25*/
   26:- use_module(library(mcintyre)).   27:- if(current_predicate(use_rendering/1)).   28:- use_rendering(c3).   29:- endif.   30:- mc.   31:- begin_lpad.   32
   33kf_fin(N,O, T) :-
   34  msw(init,S),
   35  kf_part(0, N, S,O,_LS,T).
   36
   37
   38kf(N,O,LS) :-
   39  msw(init,S),
   40  kf_part(0, N, S,O,LS,_T).
   41
   42
   43kf_part(I, N, S, [V|RO], [S|LS], T) :-
   44  I < N, 
   45  NextI is I+1,
   46  trans(S,I,NextS),
   47  emit(NextS,I,V),
   48  kf_part(NextI, N, NextS, RO, LS, T).
   49
   50kf_part(N, N, S, [], [], S).
   51
   52trans(S,_I,NextS) :-
   53  {NextS =:= E + S},
   54  msw(trans_err,E).
   55
   56emit(NextS,_I,V) :-
   57  {V =:= NextS +X},
   58  msw(obs_err,X).
   59
   60values(init,real).
   61values(trans_err,real).
   62values(obs_err,real).
   63
   64:- set_sw(init, norm(0,1)).   65:- set_sw(trans_err,norm(0,2)).   66:- set_sw(obs_err,norm(0,1)).   67
   68% prior as in Russel and Norvig 2010, Fig 15.10
   69% transition noise as in Russel and Norvig 2010, Fig 15.10
   70% observation noise as in Russel and Norvig 2010, Fig 15.10
   71
   72:- end_lpad.   73
   74
   75
   76hist(Samples,NBins,Chart):-
   77  mc_sample_arg(kf_fin(1,_O1,Y),Samples,Y,L0),
   78  histogram(L0,Chart,[nbins(NBins)]).
   79% plot the density of the state at time 1 in case of no observation (prior)
   80% and in case of observing 2.5.
   81% Observation as in Russel and Norvig 2010, Fig 15.10