#include "TRandom.h"
void Bayes(long N=1000000){
   double Pm= 0.00002; // a priori
   //double Pm= 0.001; // a priori
   double Ppm=0.993;  // P(positivo|malato)
   double Pns=0.9999; // P(negativo|sano)
   //TRandom3 R(123457);
   TRandom3 R(1556457);

   long Pos=0;      // contatore esiti positivi 
   long Neg=0;      // contatore esiti positivi 
   long MVero=0;    // contatore malati (veri) 

   long FalsiNegativi=0;
   long VeriNegativi =0;
   long FalsiPositivi=0;
   long VeriPositivi =0;

   for (int i=0; i<N ; i++) {
    int     Sog=0;  // =0 -> sano ; =1 -> malato
    double  x=R.Uniform();
    if (x<Pm) {
      Sog=1;   // malato
      MVero++; // contiamo i veri malati  
    }
    x=R.Uniform();
    if (Sog==0)
      {
	// se soggetto sano
	if (x<Pns) {
	  // esito negativo per soggetto sano 
	  VeriNegativi++;
	  Neg++; // esisto negativo 
	}
	else
	  {
	    Pos++; // esisto positivo per soggetto sano 
	    FalsiPositivi++;
	  }
      }
    else
      {
	// se soggetto malato
	if (x<Ppm) {
	  // esito positivo per soggetto malato  
	  Pos++;  // esisto positivo
	  VeriPositivi++;
	}
	else
	  {
	    Neg++; // esisto negativo 
	    FalsiNegativi++;
	  }
      }
   }


   printf("\n-----------------------------------------\n");
   printf("Incidenza dell'epidemia:  %f \n",Pm);
   printf("Sensibilita' del test (P(pos|malato)) =  %f \n",Ppm);
   printf("Specificita' del test (P(neg|sano))   =  %f \n",Pns);
   printf("-----------------------------------------\n\n");
   printf("Nel campione analizzato i sani  erano %li e i malati erano  %li\n",N-MVero,MVero);
   printf("Al test risultano  %li positivi e %li N-positivi \n",Pos,N-Pos);
   printf("Al test risultano  %li positivi e %li   negativi \n\n",Pos,Neg);
   printf("Risultati:   %li VeriPositivi e %li FalsiPositivi \n", VeriPositivi ,FalsiPositivi);
   printf("Risultati:   %li VeriNegativi e %li FalsiNegativi \n", VeriNegativi ,FalsiNegativi);
   printf("Somma dei 4 casi: %li \n\n", VeriPositivi+FalsiPositivi+VeriNegativi+FalsiNegativi);

   printf("Risultati:   %e VeriPositivi/Positivi e %e FalsiPositivi/Positivi \n", double(VeriPositivi)/double(Pos) ,  double(FalsiPositivi)/double(Pos));
   printf("Risultati:   %e VeriNegativi/Negativi e %e FalsiNegativi/Negativi \n", double(VeriNegativi)/double(N-Pos) ,double(FalsiNegativi)/double(N-Pos));
   printf("-----------------------------------------\n");
   double Pps = (1.-Pns);
   double Pmp = Ppm*Pm      / (Ppm*Pm + Pps*(1-Pm));
   double Psp = Pps*(1.-Pm) / (Ppm*Pm + Pps*(1-Pm));

   // probablilita' di malato o sano se l'esito e' negativo
   double Pnm = (1.-Ppm);
   double Pmn = Pnm*Pm      / (Pnm*Pm + Pns*(1-Pm));
   double Psn = Pns*(1.-Pm) / (Pnm*Pm + Pns*(1-Pm));

   printf("Secondo Bayes:   %e P(VeriPositivi) %e P(FalsiPositivi) \n",Pmp,Psp);
   printf("Secondo Bayes:   %e P(VeriNegativi) %e P(FalsiNegativi) \n",Psn,Pmn);
   printf("-----------------------------------------\n\n");
 
   
   printf("Stima di Pm = Pos/N = %e \n",double(Pos)/double(N));

   
}