/* 

  Bayes (Resampling Stats) in Picat.

  From Statistics101 (Resampling Stats)
  Program bayes.txt
  """
  Suppose 2 percent of the population has a particular disease.
  There is a test for the disease, but it is not perfectly accurate,
  as 3.2 percent of the population tests positive for the disease.
  There is a 75 percent chance that a person with the disease will
  test positive. Calculate the probability that a person who does 
  test positive really does have the disease.
  Solution: Out of 1000 people, 32 (3.2 percent) will test positive;
  but only 20 (2 percent) actually have the disease.

  NAME hasDisease noDisease
  NAME truePositive falsePositive
  COPY 2#hasDisease 98#noDisease population
  COPY 20#truePositive 12#falsePositive testPopulation

  REPEAT 10000
    SAMPLE 1 testPopulation testResult
    IF testResult = truePositive
      SAMPLE 1 population actualHealth
      IF actualHealth = hasDisease
         SCORE hasDisease results
      ELSE
         SCORE noDisease results
      END
    END
  END
  COUNT results = hasDisease truePositivesCount
  COUNT results = noDisease falsePositiveCount
  ADD truePositivesCount falsePositiveCount totalPositives
  DIVIDE truePositivesCount totalPositives probability
  PRINT probability

  -> probability: 0.019688790092092727
  """

  Cf my Gamble model gamble_bayes.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 : test tesult
  Probabilities:
  true_positive: 0.62109999999999999 (6211 / 10000)
  false_positive: 0.37890000000000001 (3789 / 10000)
  mean = [true_positive = 0.6211,false_positive = 0.3789]

  var : actual health
  Probabilities:
  no_disease: 0.98099999999999998 (981 / 1000)
  has_disease: 0.01900000000000000 (19 / 1000)
  mean = [no_disease = 0.981,has_disease = 0.019]


*/
go ?=>
  reset_store,
  % 2% of the population has the disease
  Population = ones(2,has_disease) ++ ones(98,no_disease),
  TestPopulation = ones(20,true_positive) ++ ones(12,false_positive),
  run_model(10_000,$model(Population,TestPopulation),[show_probs_rat_trunc,mean]),
  nl,
  % fail,
  nl.
go => true.
  
  
model(Population,TestPopulation) =>
  TestResult = uniform_draw(TestPopulation),
  ActualHealth = uniform_draw(Population),

  add("test tesult",TestResult),
  add("actual health",ActualHealth).