jacknifetest.cpp

Go to the documentation of this file.

#include "MTStation.h"
#include <boost/random/lagged_fibonacci.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/generator_iterator.hpp>
#include <vector>
#include <iostream>
#include <ctime>
#include <complex>
#include <algorithm>
#include <numeric>
#include <cmath>
 #include <fstream>
 #include "Jacknife.h"
 #include "Bootstrap.h"
 #include "MTTensor.h"
 #include <time.h>

void wait ( int seconds )
{
  clock_t endwait;
  endwait = time (0) + seconds ;
  while (time(0) < endwait) {}
}

 double CalcRhoa(const std::complex<double> Z, const double frequency)
{
        return 1./(2 * PI * frequency) * mu * (pow(Z.real(),2)+pow(Z.imag(),2)) * pow(1000.0,2);
    //return Z.real();
}
 
 
int main()
{
   int nsamples =1000;
  
   std::cout<< "Number of samples: " ;
   std::cin >> nsamples;
   
   std::vector<double> RhoValues(nsamples);
   boost::lagged_fibonacci607 generator(static_cast<unsigned int>(std::time(0)));
  
   CMTStation Station;
   
   Station.GetData("bel.mtt");
   CJacknife Jacknife(nsamples);
   wait(2);
   CBootstrap Bootstrap(nsamples);
   double JackMean, JackVar;
   double BootMean, BootVar;
   Jacknife.CalcErrors(&MTTensor::GetRhoxy,Station.GetMTData().at(0),JackMean,JackVar);
   Bootstrap.CalcErrors(&MTTensor::GetRhoxy,Station.GetMTData().at(0),BootMean,BootVar);
   const double Zerr = Station.GetMTData().at(0).GetdZxy();
   const double frequency = Station.GetMTData().at(0).GetFrequency();
   const std::complex<double> Zmean(Station.GetMTData().at(0).GetZxy()); 
   
   
  /* 
   boost::normal_distribution<> real_dist(Station.GetMTData().at(0).GetZxx().real(),Station.GetMTData().at(0).GetdZxx())
        , imag_dist(Station.GetMTData().at(0).GetZxx().imag(),Station.GetMTData().at(0).GetdZxx());
   //boost::uniform_real<> real_dist(Zmean.real()-Zerr,Zmean.real()+Zerr), imag_dist(Zmean.imag()-Zerr,Zmean.imag()+Zerr);
   boost::uniform_int<> int_dist(0,nsamples-1);
   boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > randomreal(generator, real_dist);
   boost::variate_generator<boost::lagged_fibonacci607&, boost::normal_distribution<> > randomimag(generator, imag_dist);
   boost::variate_generator<boost::lagged_fibonacci607&, boost::uniform_int<> > randomindex(generator, int_dist);
   std::ofstream zvalues("zvalues");  
   std::cout << "Rhoxx: " << Station.GetMTData().at(0).GetRhoxx() << std::endl;
   for(int i = 0; i < nsamples; i++)
   {
           std::complex<double> number(randomreal(),randomimag());
       RhoValues.at(i) = CalcRhoa(number,frequency);
       zvalues << number.real() << std::endl;
       //zvalues << " " << number.imag() << std::endl;
   }
   std::ofstream rhovalues("rhovalues");
  
   for(int i = 0; i < nsamples; i++)
           rhovalues << RhoValues.at(i) << std::endl; 
        
        std::vector<double> Pseudo(RhoValues);
        std::ofstream pseudovalues("pseudovalues");
        double initmean = std::accumulate(RhoValues.begin(),RhoValues.end(),0.)/RhoValues.size();
        for (int i = 0; i < nsamples; ++i)
        {
                
                std::vector<double> Sample(RhoValues);
                std::vector<double>::iterator it(Sample.begin());
                std::advance(it,i);
                Sample.erase(it);
                Pseudo.at(i) = nsamples * initmean - std::accumulate(Sample.begin(),Sample.end(),0.);
                pseudovalues << Pseudo.at(i) << std::endl;
        }
        double jackmean = std::accumulate(Pseudo.begin(),Pseudo.end(),0.)/nsamples;
        double sumdev = 0;
        for (int i = 0; i < nsamples; ++i)
           sumdev += pow(Pseudo.at(i) - jackmean,2);
        double jackvar =sqrt( 1./double(nsamples -1) * sumdev);
        std::cout << "Jackmean: " << jackmean << " JackVar: " << jackvar << std::endl;
        */
        double drhoa = 1./(PI * frequency) * mu * abs(Zmean) * Zerr * pow(1000.0,2);
        double rhoa = CalcRhoa(Zmean,frequency);
        std::cout << "Propmean: " << rhoa << " Properr: " << drhoa << std::endl;
        std::cout << "Jackmean: " << JackMean << " Jackerr: " << JackVar << std::endl;
        std::cout << "Bootmean: " << BootMean << " Booterr: " << BootVar << std::endl;
}

---