bootstraptest.cpp

Go to the documentation of this file.

#include "CMTStation.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 <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <iostream>
#include <ctime>
#include <complex>
#include <algorithm>
#include <cmath>
 #include <fstream>
 #include "VecMat.h"
 
 namespace ublas = boost::numeric::ublas;
 
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;
   int nbootstrap = 1000;
   std::cout<< "Number of samples: " ;
   std::cin >> nsamples;
   std::cout<< "Number of bootstraps: " ;
   std::cin >> nbootstrap;
   
   gplib::rvec RhoValues(nsamples);
   boost::lagged_fibonacci607 generator(static_cast<unsigned int>(std::time(0)));
   const double Zerr = 0.1;
   const double frequency = 0.1;
   const std::complex<double> Zmean(1.,1.);
   CMTStation Station;
   boost::normal_distribution<> real_dist(Zmean.real(),Zerr), imag_dist(Zmean.imag(),Zerr);
   //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");  
   for(int i = 0; i < nsamples; i++)
   {
           std::complex<double> number(randomreal(),randomimag());
       RhoValues(i) = CalcRhoa(number,frequency);
       zvalues << number.real() << std::endl;
       //zvalues << " " << number.imag() << std::endl;
   }
   std::ofstream rhovalues("rhovalues");
   gplib::rvec Means(nbootstrap);
   for(int i = 0; i < nsamples; i++)
           rhovalues << RhoValues(i) << std::endl; 
        std::ofstream meanvalues("meanvalues");
        for (int i = 0; i < nbootstrap; ++i)
        {
                gplib::rvec Sample(RhoValues);
                for(int j = 0; j < nsamples; j++)
                {
                        int index = randomindex();
                        //std::cout << "j: " << j << " Sample(j): " << Sample(j) << " Index: "<< index <<" Sample(index): " << RhoValues(index) << std::endl;
                    Sample(j) = RhoValues(index);
                }
                //std::cout << Sample << std::endl;
                Means(i) = ublas::sum(Sample)/double(nsamples);
                meanvalues << Means(i) << std::endl;
                //std::cout << Means(i) << std::endl;
        }
        double bootmean = sum(Means)/double(nbootstrap);
        double sumdev = 0;
        for (int i = 0; i < nbootstrap; ++i)
           sumdev += pow(Means(i) - bootmean,2);
        double bootvar =sqrt( 1./double(nbootstrap -1) * sumdev);
        double rhomean = ublas::sum(RhoValues)/double(nsamples);
        double rhovar = 0;
        for (int i = 0; i < nsamples; ++i)
           rhovar += pow(RhoValues(i) - rhomean,2);
        rhovar =sqrt( 1./double(nsamples -1) * rhovar);
        std::cout << "Bootmean: " << bootmean << " BootVar: " << bootvar << std::endl;
        std::cout << "Rhomean: " << rhomean << " RhoVar: " << rhovar << std::endl;
        std::sort(Means.begin(),Means.end());
        int lowerbound = floor(nbootstrap * 0.16);
        int upperbound = ceil(nbootstrap * 0.84);
        std::cout << "Bootmedian: " << Means(nbootstrap/2) << " BootLower: " << Means(lowerbound) <<
           " BootUpper: " << Means(upperbound)  << " Booterr: " <<  Means(upperbound) - Means(lowerbound)<< std::endl;
        double drhoa = 1./(PI * 0.1) * mu * abs(Zmean) * Zerr * pow(1000.0,2);
        double rhoa = CalcRhoa(Zmean,frequency);
        std::cout << "Propmean: " << rhoa << " Properr: " << drhoa << std::endl;
}

---