/* 

  Bugs book example 2.6.1 in Picat.

  Example 2.6.1 Heart transplant cost-effectiveness: risks assumed known

  Cf my Gamble model gamble_bugs_book_2_6_1.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.
% import ordset.

main => go.

/*

  var : sn
  Probabilities (truncated):
  4.716880782204347: 0.0005000000000000
  4.675344993581804: 0.0005000000000000
  3.925798726387786: 0.0005000000000000
  3.315215500801024: 0.0005000000000000
  .........
  0.000542158025204: 0.0005000000000000
  0.000245321411756: 0.0005000000000000
  0.000097872488924: 0.0005000000000000
  0.000094193058129: 0.0005000000000000
  mean = 0.504821
  HPD intervals:
  HPD interval (0.84): 0.00009419305813..0.91789005630649

  var : ot
  Probabilities:
  1: 0.8110000000000001
  0: 0.1890000000000000
  mean = 0.811
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..1.00000000000000

  var : sp
  Probabilities (truncated):
  1.797045211659614: 0.0005000000000000
  1.373157739128549: 0.0005000000000000
  1.343972065779569: 0.0005000000000000
  1.204797564147059: 0.0005000000000000
  .........
  0.000288283809216: 0.0005000000000000
  0.000269484385251: 0.0005000000000000
  0.000255497412123: 0.0005000000000000
  0.00019950625418: 0.0005000000000000
  mean = 0.192931
  HPD intervals:
  HPD interval (0.84): 0.00019950625418..0.35198453688502

  var : st
  Probabilities (truncated):
  0.0: 0.1890000000000000
  1.797045211659614: 0.0005000000000000
  1.373157739128549: 0.0005000000000000
  1.343972065779569: 0.0005000000000000
  .........
  0.000310358899269: 0.0005000000000000
  0.000299658180677: 0.0005000000000000
  0.000288283809216: 0.0005000000000000
  0.00019950625418: 0.0005000000000000
  mean = 0.153552
  HPD intervals:
  HPD interval (0.84): 0.00000000000000..0.30928916665170

  var : ic
  Probabilities (truncated):
  20000.0: 0.1890000000000000
  25391.135634978840244: 0.0005000000000000
  24119.473217385646421: 0.0005000000000000
  24031.916197338705388: 0.0005000000000000
  .........
  20000.931076697808749: 0.0005000000000000
  20000.898974542033102: 0.0005000000000000
  20000.864851427646499: 0.0005000000000000
  20000.598518762541062: 0.0005000000000000
  mean = 20460.7
  HPD intervals:
  HPD interval (0.84): 20000.00000000000000..20927.86749995510399

  var : is
  Probabilities (truncated):
  1.310351678628755: 0.0005000000000000
  1.062621777183212: 0.0005000000000000
  1.037093812041014: 0.0005000000000000
  0.95881599289762: 0.0005000000000000
  .........
  -3.119740854594099: 0.0005000000000000
  -3.834992112969625: 0.0005000000000000
  -4.423969623829041: 0.0005000000000000
  -4.593860916842535: 0.0005000000000000
  mean = -0.351269
  HPD intervals:
  HPD interval (0.84): -0.97307344859795..0.26708363404347

  var : r
  Probabilities (truncated):
  61845910.762939512729645: 0.0005000000000000
  30483331.794227954000235: 0.0005000000000000
  9877393.037721861153841: 0.0005000000000000
  9775493.481235705316067: 0.0005000000000000
  .........
  -11447167.85719901882112: 0.0005000000000000
  -14281207.398582777008414: 0.0005000000000000
  -38647528.271766297519207: 0.0005000000000000
  -204347516.037253469228745: 0.0005000000000000
  mean = -141912.0
  HPD intervals:
  HPD interval (0.84): -339470.96323667030083..262875.62472691887524


*/
go ?=>
  reset_store,
  run_model(2_000,$model,[show_probs_trunc,mean,show_hpd_intervals,hpd_intervals=[0.84]]),
  nl,
  % show_store_lengths,nl,
  % fail,
  nl.
go => true.
  
  
model() =>

  SN = exponential_dist(1/0.5), %  life without transplant (mean 2)
  OT = bern(0.8),             %  survive operation (prob 0.8)
  SP = exponential_dist(1/0.2), %  life if survive transplant (mean 5)
  ST = OT * SP,               %  total life time if choose transplant
  IC = 20000 + 3000 * ST,     %  total additional cost of transplant
  IS = ST - SN,               %  total additional survival
  R  = IC / IS,               %  individual cost per additional year

  add("sn",SN),
  add("ot",OT),
  add("sp",SP),
  add("st",ST),
  add("ic",IC),
  add("is",IS),
  add("r",R).