JTank.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
 
#include "TRandom3.h"
 
#include "JTools/JStats.hh"
#include "JMath/JMathSupportkit.hh"
 
#include "Jeep/JParser.hh"
 
 
/**
 * \file
 * Test application for estimate of maximam number from limited sampling.
 * This estimate is known as the German tank problem.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  size_t    N;
  size_t    k;
  size_t    start;
  size_t    M;
  Int_t     seed;
  int       debug;
 
  try { 
 
    JParser<> zap;
 
    zap['N'] = make_field(N);
    zap['k'] = make_field(k);
    zap['s'] = make_field(start)          = 1;
    zap['M'] = make_field(M)              = 1;
    zap['S'] = make_field(seed)           = 0;
    zap['d'] = make_field(debug)          = 3;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  gRandom->SetSeed(seed);
 
 
  /*
  // https://en.wikipedia.org/wiki/German_tank_problem
  cout << getMinimumVarianceUnbiasedEstimator(60, 4) << endl;
  cout << getBayesianMedian(60, 4) << endl;
  cout << getBayesianMean(60, 4) << endl;
  */
 
  vector<size_t> buffer;
 
  for (size_t i = 0; i != N; ++i) {
    buffer.push_back(start + i);                 // serial numbers
  }
 
  JStats Q[] = {
    JStats("frequentist"),
    JStats("Bayesian median"),
    JStats("Bayesian mean")
  };
 
  for (size_t counter = 0; counter != M; ++counter) {
 
    size_t m = 0;
 
    for (size_t i = 0; i != k; ++i) {
 
      const size_t l = gRandom->Integer(N - i);  // remaining serial numbers
 
      if (buffer[l] > m) {
        m = buffer[l];
      }
 
      swap(buffer[l], buffer[N - i - 1]);        // remove observation
    }
 
    Q[0].put(getMinimumVarianceUnbiasedEstimator(m, k, start));
    Q[1].put(getBayesianMedian(m, k, start));
    Q[2].put(getBayesianMean(m, k, start));
  }
 
  if (debug >= debug_t) {
 
    for (size_t i = 0; i != sizeof(Q)/sizeof(Q[0]); ++i) {
 
      cout << setw(24) << left << Q[i].getTitle() << right << FIXED(7,3) << Q[i].getMean();
 
      if (Q[i].getCount() > 1) {
        cout << ' ' << FIXED(7,3) << Q[i].getSTDev();
      }
 
      cout << endl;
    }
  }
}