// Este Scripts compara la performance de Bootstrap contra Toys Montecarlo
// Ayuda a visualizar que tamaño de muestra y número de réplicas son suficientes
// para obtener una estimación razonable de la incerteza en el promedio de una muestra Poissoniana

void bootstrap_vs_toy(){
	
	int sample_size=20;
	int replicas=200;
	
	double mu=50;
	double sigma = sqrt(50);
	double sigmam = sigma/sqrt(sample_size);
	int emin=mu-3*sigma;
	int emax=mu+3*sigma;
	int nsigmas=10;
	 	
	int bines = (emax-emin);
	
    vector<double> population;  
    vector<double> sample;  
    vector<double> bootstrap_sample; 
     
	TRandom3 X_Random_variable(0); 

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////

    gROOT->SetStyle("Plain");
    TCanvas * canvas = new TCanvas("canvas","canvas",2000,1000);
    canvas->Divide(2,2);

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////

    canvas->cd(1);
    canvas->cd(1)->SetGridx();
    canvas->cd(1)->SetGridy();
    
    TF1 *pdf = new TF1("pdf","TMath::Poisson(x,50)",emin,emax);
    pdf->SetLineColor(kRed);
    pdf->Draw();


////////////////////////////////////////////////////////////////////////////////////////////////////
// Aca empiezo los Toys
		
	TH1D * histo_toy  =  new TH1D("histo_toy","", bines, mu-nsigmas*sigmam, mu+nsigmas*sigmam);
	histo_toy -> Sumw2();
	
    canvas->cd(2);
    canvas->cd(2)->SetGridx();
    canvas->cd(2)->SetGridy();

////////////////////////////////////////////////////////////////////////////////////////////////////
// Hago los toys   
   
	double toys_mean;
    for (int i = 1; i < replicas+1; ++i) {
		double t_sample;
		double t_sample_suma;
		
		for (int j = 1; j < sample_size+1; ++j) {
			t_sample = X_Random_variable.Poisson(mu);
			t_sample_suma+= t_sample;
		}
		double mean=double(t_sample_suma)/sample_size;
		t_sample_suma=0;
		double statistic=mean;
		toys_mean+=mean;
		histo_toy->Fill(statistic);
    }
    toys_mean=toys_mean/replicas;
    cout<<"El promedio de los toys es:   "<<toys_mean<<endl;
    
    histo_toy->SetTitle("Toys del promedio");
    histo_toy->GetXaxis()->SetTitle("promedio de la replica");
	histo_toy->GetYaxis()->SetTitle("#");
	histo_toy->SetLineWidth(2);
	
	histo_toy->SetLineColor(kRed);
	histo_toy->Draw("E0 HIST");
	
////////////////////////////////////////////////////////////////////////////////////////////////////
// Aca empiezo el bootstrap
	
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
	
	TH1D * histo_sample  =  new TH1D("histo_sample","", bines, emin, emax);
	histo_sample -> Sumw2();
	
    canvas->cd(3);
    canvas->cd(3)->SetGridx();
    canvas->cd(3)->SetGridy();

    
////////////////////////////////////////////////////////////////////////////////////////////////////
// Tomo la muestra 
	
	double x_sample_mean;
    for (int i = 1; i < sample_size+1; ++i) {
		double x_sample;
		//x_sample=X_Random_variable.Uniform(emin,emax);
		x_sample=X_Random_variable.Poisson(mu);
		x_sample_mean+=x_sample;
		sample.push_back(x_sample);
		histo_sample->Fill(x_sample);
    }
    
    x_sample_mean=double(x_sample_mean)/sample_size;
    cout<<"El promedio de la muestra es: "<<x_sample_mean<<endl;
   
    histo_sample->SetTitle("Muestra");
    histo_sample->GetXaxis()->SetTitle("realizacion");
	histo_sample->GetYaxis()->SetTitle("#");
	histo_sample->SetLineWidth(2);

	histo_sample->SetLineColor(kBlue);
	histo_sample->Draw("E0 HIST");
	
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
		
	TH1D * histo_bootstrap  =  new TH1D("histo_bootstrap","", bines, mu-nsigmas*sigmam, mu+nsigmas*sigmam);
	histo_bootstrap -> Sumw2();
	
    canvas->cd(4);
    canvas->cd(4)->SetGridx();
    canvas->cd(4)->SetGridy();

////////////////////////////////////////////////////////////////////////////////////////////////////
// Hago las replicas, estas son con reposicion

	double bootstrap_mean;
    for (int i = 1; i < replicas+1; ++i) {
		double b_sample;
		double b_sample_suma;
		for (int j = 1; j < sample_size+1; ++j) {
			int indice=X_Random_variable.Uniform(1,sample_size);
			b_sample=sample[indice];
			b_sample_suma+=b_sample;
		}
		double mean=b_sample_suma/sample_size;
		b_sample_suma=0;
		double statistic=mean;
		bootstrap_mean+=mean;
		//cout<<"Estadistico: "<<statistic<<endl;
		bootstrap_sample.push_back(statistic);
		histo_bootstrap->Fill(statistic);
    }
    bootstrap_mean=bootstrap_mean/replicas;
    cout<<"El promedio de las replicas es:   "<<bootstrap_mean<<endl;
    double bias=abs(bootstrap_mean-x_sample_mean);
    cout<<"por lo tanto el bias es:          "<<bias<<endl;
    
    histo_bootstrap->SetTitle("Bootstrap para el promedio");
    histo_bootstrap->GetXaxis()->SetTitle("promedio de la replica");
	histo_bootstrap->GetYaxis()->SetTitle("#");
	histo_bootstrap->SetLineWidth(2);

	histo_bootstrap->SetLineColor(kBlue);
	histo_bootstrap->Draw("E0 HIST");
	

	
}