/* 

  Medical test in Picat.

  https://www.math.hmc.edu/funfacts/ffiles/30002.6.shtml
  """
  Suppose that you are worried that you might have a rare disease. You decide to get tested,
  and suppose that the testing methods for this disease are correct 99 percent of the time
  (in other words, if you have the disease, it shows that you do with 99 percent probability,
  and if you don't have the disease, it shows that you do not with 99 percent probability).
  Suppose this disease is actually quite rare, occurring randomly in the general population
  in only one of every 10,000 people.

  If your test results come back positive, what are your chances that you actually have the disease?

  Do you think it is approximately: (a) .99, (b) .90, (c) .10, or (d) .01?

  Surprisingly, the answer is (d), less than 1 percent chance that you have the disease! 
  """


  Exact probabilities via my Gamble model:
  var : disease
  #f: 0.9901960784313726
  #t: 0.009803921568627458
  mean: 0.009803921568627458

  This is a port of my Gamble model gamble_medical_test.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.

main => go.

/* 
  * Observe test == true
    var : disease
    Probabilities:
    false: 0.9909909909909910
    true: 0.0090090090090090
    mean = [false = 0.990991,true = 0.00900901]

    var : test
    Probabilities:
    true: 1.0000000000000000
    mean = [true = 1.0]

  * Observe test == false
    var : disease
    Probabilities:
    false: 1.0000000000000000
    mean = [false = 1.0]

    var : test
    Probabilities:
    false: 1.0000000000000000
    mean = [false = 1.0]

  * Observe disease == true
    var : disease
    Probabilities:
    true: 1.0000000000000000
    mean = [true = 1.0]

    var : test
    Probabilities:
    true: 0.9230769230769231
    false: 0.0769230769230769
    mean = [true = 0.923077,false = 0.0769231]

    Another run:
    var : disease
    Probabilities:
    true: 1.0000000000000000
    mean = [true = 1.0]

    var : test
    Probabilities:
    true: 0.9565217391304348
    false: 0.0434782608695652
    mean = [true = 0.956522,false = 0.0434783]

  * Observe disease == false
    var : disease
    Probabilities:
    false: 1.0000000000000000
    mean = [false = 1.0]

    var : test
    Probabilities:
    false: 0.9900995049752488
    true: 0.0099004950247512
    mean = [false = 0.9901,true = 0.0099005]

*/
go ?=>
  reset_store(),
  run_model(100_000,model,[show_probs_trunc,mean]),
  nl,
  % fail,
  nl.
go => true.

model() =>
  %  Suppose this disease is actually quite rare, occurring randomly in the general population
  %  in only one of every 10,000 people.
  ProbOfDisease = 1/10000,
  % ProbOfDisease = 1/1000,
    
  %  Probability that a person has a disease (and there's a reason to do a test)
  Disease = flip(ProbOfDisease),
    
  %  The test is quite accurate: It shows correct result (test is positive if disease) in 99%.
  %  However, in 1% of the case it shows incorrect result (positive even if there is no disease).
  Reliability = 0.99,  %  => disease:0.009803921568627446
  % Reliability 0.999, %  A more reliable test => disease:0.09083469721767588
  % Reliability 0.9999, %  An even more reliable test => disease:0.500000000000028
    
  Test = condt(Disease,
               flip(Reliability),
               flip(1-Reliability)),
    
  % observe(Test),
  observe(not Test),
  % observe(Disease),
  % observe(not Disease),     

  if observed_ok then
    add("disease",Disease),
    add("test",Test),    
  end.