1:- use_module(library(mcintyre)).    2:- use_module(library(matrix)).    3:- use_module(library(clpfd)).    4
    5:- if(current_predicate(use_rendering/1)).    6:- use_rendering(c3).    7:- endif.    8:- mc.    9:- begin_lpad.   10
   11% http://arxiv.org/pdf/1505.02965v2.pdf
   12
   13gp(X,Kernel,Y):-
   14  compute_cov(X,Kernel,C),
   15  gp(C,Y).
   16
   17gp(Cov,Y):gaussian(Y,Mean,Cov):-
   18  length(Cov,N),
   19  list0(N,Mean).
   20
   21compute_cov(X,Kernel,C):-
   22  length(X,N),
   23  cov(X,N,Kernel,CT,CND),
   24  transpose(CND,CNDT),
   25  matrix_sum(CT,CNDT,C).
   26
   27cov([],_,_,[],[]).
   28
   29cov([XH|XT],N,Ker,[KH|KY],[KHND|KYND]):-
   30  length(XT,LX),
   31  N1 is N-LX-1,
   32  list0(N1,KH0),
   33  maplist(call(Ker,XH),XT,KH1),
   34  call(Ker,XH,XH,KXH),
   35%  kernel(XH,XH,KXH),
   36  append([KH0,[KXH],KH1],KH),
   37  append([KH0,[0],KH1],KHND),
   38  cov(XT,N,Ker,KY,KYND).
   39
   40
   41sq_exp_p(X,XP,K):-
   42  sigma(Sigma),
   43  l(L),
   44  K is Sigma^2*exp(-((X-XP)^2)/2/(L^2)).
   45
   46
   47l(L):uniform(L,[1,2,3]).
   48
   49sigma(Sigma):uniform_dens(Sigma,-2,2).
   50
   51:- end_lpad.   52
   53%kernel(X,X,K):-!,
   54%  K is 1.27^2+0.3^2.
   55
   56sq_exp(X,XP,K):-
   57  K is exp(-((X-XP)^2)/2).
   58  %K is (1.27^2)*exp(-((X-XP)^2)/2).
   59
   60min(X,XP,K):-
   61  K is min(X,XP).
   62
   63lin(X,XP,K):-
   64  K is (X*XP+5)^2.
   65
   66ou(X,XP,K):-
   67  K is exp(-abs(X-XP)).
   68draw_fun(Kernel,C):-
   69%  X=[-1.50,-1.00,-0.75,-0.40,-0.25,0.00],
   70%  X=[-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1.00,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5],
   71  X=[-3,-2,-1,0,1,2,3],
   72  compute_cov(X,Kernel,K),
   73  mc_sample_arg_first(gp(K,Y),5,Y,L),
   74  numlist(1,5,LD),
   75  maplist(name_s,L,LD,L1),
   76  C = c3{data:_{x:x, columns:[[x|X]|L1] },
   77 % legend:_{show: false},
   78  axis:_{ x:_{ tick:_{fit:false}}}}.
   79
   80draw_fun_post(Kernel,LSR):-
   81%  X=[-1.50,-1.00,-0.75,-0.40,-0.25,0.00],
   82%  X=[-4.5,-4,-3.5,-3,-2.5,-2,-1.5,-1.00,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5],
   83  numlist(0,10,X),
   84  XT=[2,4,6,8,9],
   85  YT=[1,-0.4,-0.8,0.25,0.6],
   86  mc_lw_sample_arg(gp(X,Kernel,Y),gp(XT,Kernel,YT),10,Y,L),
   87  keysort(L,LS),
   88  reverse(LS,LSR).
   89  /*
   90  numlist(1,5,LD),
   91  maplist(name_s,L,LD,L1),
   92  C = c3{data:_{x:x, columns:[[x|X]|L1] },
   93 % legend:_{show: false},
   94  axis:_{ x:_{ tick:_{fit:false}}}}.
   95*/
   96name_s(V-_,N,[ND|V]):-
   97  atomic_concat(f,N,ND).