/* 

  Linear regression in Picat.

  From  WebPPL example example/linearRegression.wppl

  Cf my Gamble model gamble_linear_regression2.rkt

  This program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Picat page: http://www.hakank.org/picat/

*/

import ppl_distributions, ppl_utils.
import util.

main => go.

/*

  var : b
  Probabilities (truncated):
  0.098272891286752: 0.0048780487804878
  0.0976680536227: 0.0048780487804878
  0.097639697054619: 0.0048780487804878
  0.097201292243307: 0.0048780487804878
  .........
  -0.095974962605678: 0.0048780487804878
  -0.096026016755556: 0.0048780487804878
  -0.096806207357495: 0.0048780487804878
  -0.097642010751878: 0.0048780487804878
  mean = 0.00498963

  var : m
  Probabilities (truncated):
  2.047109513691414: 0.0048780487804878
  2.043224314085548: 0.0048780487804878
  2.041448908522773: 0.0048780487804878
  2.040859684621964: 0.0048780487804878
  .........
  1.961568137457786: 0.0048780487804878
  1.958538106663283: 0.0048780487804878
  1.95447323808633: 0.0048780487804878
  1.954259767678693: 0.0048780487804878
  mean = 1.99932

  var : post 4
  Probabilities (truncated):
  12.559924886234032: 0.0048780487804878
  12.398210361315854: 0.0048780487804878
  11.365494224064948: 0.0048780487804878
  11.284555395389573: 0.0048780487804878
  .........
  4.612176642033882: 0.0048780487804878
  4.317635806360201: 0.0048780487804878
  4.133441850786177: 0.0048780487804878
  1.986427684394409: 0.0048780487804878
  mean = 8.02153

  var : post 5
  Probabilities (truncated):
  15.753529837842535: 0.0048780487804878
  15.239777079720653: 0.0048780487804878
  14.00784469004777: 0.0048780487804878
  13.625753270644541: 0.0048780487804878
  .........
  7.567353443015014: 0.0048780487804878
  7.324728964073918: 0.0048780487804878
  6.082954826172934: 0.0048780487804878
  5.836018244084045: 0.0048780487804878
  mean = 10.1238

  var : sigma
  Probabilities (truncated):
  3.643719813048943: 0.0048780487804878
  3.524760318527321: 0.0048780487804878
  3.466501487178647: 0.0048780487804878
  3.460750270548031: 0.0048780487804878
  .........
  0.014179505691496: 0.0048780487804878
  0.009982252014324: 0.0048780487804878
  0.00671375848899: 0.0048780487804878
  0.002868917976992: 0.0048780487804878
  mean = 0.965513


  Variable lengths:
  b = 205
  m = 205
  post 4 = 205
  post 5 = 205
  sigma = 205


*/
go ?=>
  reset_store,
  run_model(1_000_000,$model,[show_probs_trunc,mean]),
  nl,
  show_store_lengths,
  % fail,
  nl.
go => true.
  
  
model() =>
  Xs = 0..4,
  Ys = [0,2,4,6,8],

  M = normal_dist(0,1),
  B = normal_dist(0,1),
  Sigma = gamma_dist(1,1),

  ObsCrit = 0.1,

  foreach(I in 1..Xs.len)
    observe(abs(Ys[I]-(B+M*Xs[I])) <= ObsCrit)
  end,

  % predictions
  Post4=normal_dist(B + M*4,Sigma),
  Post5=normal_dist(B + M*5,Sigma),  
 
  
  if observed_ok then
    add("m",M),
    add("b",B),
    add("sigma",Sigma),
    add("post 4",Post4),
    add("post 5",Post5),    
  end.