Calculation of Equation (7) in our sampling paper.

main.cpp

#include <iostream>
#define _USE_MATH_DEFINES
#include <cmath>
#include <boost/math/special_functions/beta.hpp>

double Q_k0( double f0, double k0, double k ) {
  double c = 1.0 / ( f0 * (k0+1.0) );

  double x = f0 / (1.0-f0);
  double a = k + 1.0;
  double b = k0 - k + 1.0;
  double ibeta = boost::math::ibeta( a, b, x );

  return c * ibeta;
}

double Gaussian( double sigma, double mu, double x ) {
  double c = 1.0 / std::sqrt( 2.0 * M_PI * sigma*sigma );
  double d = std::exp( - (x-mu)*(x-mu) / (2.0*sigma*sigma) );
  return c*d;
}

double lognormal( double mu, double sigma, double x ) {
  double c = 1.0 / (x*sigma * std::sqrt(2.0*M_PI) );
  double d = std::exp( - (std::log(x)-mu)*(std::log(x)-mu) / (2.0*sigma*sigma) );
  return c*d;
}

double Levy( double c, double mu, double x ) {
  if( x < mu ) { return 0.0; }

  double y1 = std::sqrt( c/(2.0*M_PI) ) / std::pow( x-mu, 1.5 );
  double y2 = exp( -c/(2.0*(x-mu)) );
  return y1 * y2;
}

// return original degree distribution
double P_0( double k ) {
  // return Gaussian( 150.0, 50.0, k );
  // return lognormal( std::log(200.0), std::log(2.0), k );
  return Levy( 150.0, 0.0, k );
}

double P( double k, double f0 ) {
  const double k0_max = 10000;

  double sum = 0.0;
  for( double i=k; i < k0_max; i+=1.0 ) {
    sum += P_0(i) * Q_k0( f0, i, k );
  }
  return sum;
}

int main() {
  /*
  for( long i=0; i < 500; i++) {
    std::cout << i << ' ' << P_0( static_cast<double>(i) ) << std::endl;
  }
  */
  const long max = 500;
  for( long i=0; i < max; i++) {
    std::cout << i << ' ' << P( static_cast<double>(i), 0.1 ) << ' ' << P_0( static_cast<double>(i) ) << std::endl;
    // std::cout << i << ' ' << P( static_cast<double>(i), 0.1 ) << ' ' << Q_k0( 0.1, 150.0, static_cast<double>(i)) << std::endl;
  }
  for( long i=max; i < 1000; i++) {
    std::cout << i << ' ' << 0.0 << ' ' << P_0( static_cast<double>(i) ) << std::endl;
  }


  return 0;
}