main.cpp

#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <fstream>
#include <iomanip>
#include <cmath>

#include "statusbar.h"

using namespace std;

//double seed=time(NULL);
double seed=33149495;
double p(double x);
int yearrand();

int main(int argc, char* argv[]){
if(argc>1) if(atoi(*(argv+1))==0) return 0;
	srand(seed);
	statusbar mybar;
	int n=500000;

//I experiment
	ofstream save;		  //open file
  save.open("montecarlo.csv", ios::out);
  while(!save.is_open()) cout <<"error creating file!"<<endl;	
int year[365]={0}, people=55;
mybar.maxevents=people;
int cntr=0;
 for(int i=2;i<=people;i++){//table per people
   cntr=0;
   for(int j=0;j<n;j++){//generates n events
     for(int l=0;l<365;l++) year[l]=0;	//reset year
     for(int k=0;k<i;k++) ++year[yearrand()%365]; //put i random birthday in year
     for(int l=0;l<365;l++) if(year[l]>=2) {cntr++;break;}//check if in this of n cases there are 2 people with same birthday
   }
	 save<<i<<";"<<(double)cntr/n<<endl;
   cout<<i<<" people\tp:  "<<setw(6)<<(double)cntr/n*100<<"%"<<"\t";
	 mybar.drawwpercent(i);cout<<endl;

}
	
	cout<<"\n\n***************************\n\n";
	//II experiment
	double average=0;
	mybar.setmaxevents(n);
		mybar.setr(1);
	for(int j=1; j<=n;j++){//n simulations
		for(int l=0;l<365;l++) year[l]=0; //reset
		cntr=0;
		
	while(true){
		cntr++;
		if(++year[yearrand()%365]==2) break;	 //count how many people before two have same birthday
	}
		average=((average*(j-1))+cntr)/j;
		mybar.draw(j);
}
	cout<<"\nWith "<<n<<" simulation, I need a average of "<<average<<" people to be sure two have same birthday\n\n";
  return 0;
}


/*
 function of probablity based on italian born in 2014
 data source: http://demo.istat.it/altridati/IscrittiNascita/2014/T2.8.xls
*/
 double p(double x){
	return (-3.27*pow(10.,-9)*x*x*x+1.86*pow(10.,-6)*x*x-2.63*pow(10.,-4)*x+0.0889)/15.643;
}

int yearrand(){
	int n=500000,i=0;
	double xi=0,yi=0;  
	int a=0,b=365;
	double fmax=p(9);
	while(i<n){  
		xi=rand()%(b-a+1)+a;
		yi=rand()/(double)RAND_MAX*fmax;
		if(p(xi)>yi) {return xi; i++;}
	}
	return -1;
}