// Bootstrap para la correlación

void correlacion(vector<double> xdata, vector<double> ydata, double & corr_xy, double & ucorr_xy){

   double xsum = std::accumulate(xdata.begin(), xdata.end(), 0.0);
   double xmean = xsum / xdata.size();
   //cout<<"x mean: "<<xmean<<endl;
    
   double ysum = std::accumulate(ydata.begin(), ydata.end(), 0.0);
   double ymean = ysum / ydata.size();
   //cout<<"y mean: "<<ymean<<endl;
    
   double numerador=0;
   for (int i = 0; i < xdata.size(); ++i) numerador+=(xdata[i]-xmean)*(ydata[i]-ymean);
    
   double sumax=0;
   double sumay=0;
   for (int i = 0; i < xdata.size(); ++i) sumax+=pow(xdata[i]-xmean,2);
   for (int i = 0; i < xdata.size(); ++i) sumay+=pow(ydata[i]-ymean,2);
	double denominador=pow(sumax*sumay,0.5);
	
   corr_xy = numerador/denominador;
	ucorr_xy = (1-pow(corr_xy,2))/pow(xdata.size()-3,0.5); // Bajo la aproximacion de que x, y es normal bivariada
  
}

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

void bootstrap_correla(){
	
	const int population_size=19;
	const int sample_size=7;
	int replicas=1000;
     
	// Cargamos los datos de la población
   double tiempo[population_size]={45, 40, 35, 20, 60, 30, 60, 40, 35, 45, 35, 15, 45, 20, 30, 60, 25, 15, 45};
   double distancia[population_size]{10, 9, 6, 3.1, 12.7, 10, 9, 7.6, 5.8, 16, 4.9, 6, 6.2, 5, 5, 10, 7.5, 2, 5};
   int alturas[population_size]{536, 5477, 1175, 2229, 6751, 663, 2746, 5388, 4752, 2520, 4041, 1035, 2260, 2397, 492, 3969, 3161, 2552, 833};

   // Elegimos los índices de la población que determinan la muestra
   int indice[sample_size]={1, 5, 6, 7, 9, 12, 16};

   vector<double> xpoblacion;  
   vector<double> ypoblacion;  
   vector<double> xsample;  
   vector<double> ysample;  
   vector<double> bootstrap_sample; 

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

   TRandom3 X_Random_variable(0); 

   gROOT->SetStyle("Plain");
   TCanvas * canvas = new TCanvas("canvas","canvas",2000,500);
   canvas->Divide(3,1);

   canvas->cd(1);
   canvas->cd(1)->SetGridx();
   canvas->cd(1)->SetGridy();
   
   for (int i = 0; i < population_size; ++i) {
		xpoblacion.push_back(tiempo[i]);
		ypoblacion.push_back(distancia[i]);
   }
   
   TGraph * gr_population = new TGraph(xpoblacion.size(), &xpoblacion[0], &ypoblacion[0]);
   gr_population->SetTitle("Poblacion");
   gr_population->GetXaxis()->SetTitle("Tiempo (minutos)");
   gr_population->GetYaxis()->SetTitle("Distancia (km)"); 
   gr_population->SetMarkerStyle(21);

   gr_population->Draw("AP");
    
   double corr_xy;
   double ucorr_xy;
   correlacion(xpoblacion, ypoblacion, corr_xy, ucorr_xy);
   cout<<"Correlacion de la poblacion: "<<corr_xy<<"  +/- "<<ucorr_xy<<endl;
    
////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
	
   canvas->cd(2);
   canvas->cd(2)->SetGridx();
   canvas->cd(2)->SetGridy();
    
////////////////////////////////////////////////////////////////////////////////////////////////////
// Tomo la muestra sin reposicion

    for (int i = 0; i < sample_size; ++i) {
      xsample.push_back(xpoblacion[indice[i]]);
      ysample.push_back(ypoblacion[indice[i]]);
    }
   
    TGraph * gr_sample= new TGraph(xsample.size(), &xsample[0], &ysample[0]);
    gr_sample->SetTitle("Muestra");
    gr_sample->GetXaxis()->SetTitle("Tiempo (minutos)");
    gr_sample->GetYaxis()->SetTitle("Distancia (km)"); 
    gr_sample->SetMarkerStyle(21);

    gr_sample->Draw("AP");
    
    double corr_xy_sample;
    double ucorr_xy_sample;
    correlacion(xsample, ysample, corr_xy_sample, ucorr_xy_sample);
    cout<<"Correlacion de la muestra: "<<corr_xy_sample<<"  +/- "<<ucorr_xy_sample<<endl;

////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////
		
	TH1D * histo_bootstrap  =  new TH1D("histo_bootstrap","", 1000, 0, 1);
	histo_bootstrap -> Sumw2();
	
   canvas->cd(3);
   canvas->cd(3)->SetGridx();
   canvas->cd(3)->SetGridy();
    
   vector<double> b_xsample;  
   vector<double> b_ysample; 
   double corr_xy_replica;
	double ucorr_xy_replica;
    
   for (int i = 0; i < replicas; ++i) {
		double b_sample;
		double b_sample_suma;
		for (int j = 0; j < sample_size; ++j) {
			int indice=X_Random_variable.Uniform(1,sample_size);
			b_xsample.push_back(xsample[indice]);
			b_ysample.push_back(ysample[indice]);	
		}

		correlacion(b_xsample, b_ysample, corr_xy_replica, ucorr_xy_replica);
		
		double statistic=corr_xy_replica;
		bootstrap_sample.push_back(statistic);
		histo_bootstrap->Fill(statistic);
   }
    
   double ucorr_xy_bootstrap;
   double sum = std::accumulate(bootstrap_sample.begin(), bootstrap_sample.end(), 0.0);
	double mean = sum / replicas;
	
	double var=0;
	for (int i = 0; i < replicas; ++i) {
		var+=pow(bootstrap_sample[i]-mean,2)/(replicas-1);
	}
	ucorr_xy_bootstrap=pow(var,0.5);
	
   histo_bootstrap->SetTitle("Histograma del rho de las replicas");
   histo_bootstrap->GetXaxis()->SetTitle("correlacion de la replica");
	histo_bootstrap->GetYaxis()->SetTitle("#");
	histo_bootstrap->SetLineWidth(2);

	histo_bootstrap->Draw("HIST");
	
	double ucorr_xy_bootstrap2=histo_bootstrap->GetStdDev();
	cout<<"Incerteza en la correlación obtenida por bootstrap: "<<ucorr_xy_bootstrap<<endl;
	cout<<"Incerteza en la correlación obtenida por bootstrap del histograma: "<<ucorr_xy_bootstrap2<<endl;
	
}