538

riddler_jan15.R

library("Rcpp")
library("microbenchmark")

simulate_throws <- cppFunction('
  double simulate_throw(NumericVector start, int nthrows, int trials) {
      
      int trials_so_far = 0;
      int final_throw_successes_so_far = 0;
      int hits_so_far = 0;
      int nstart = start.size();
      
      double success_rate = 0;
      NumericVector out(nthrows);

      for (int i = 0; i < nstart; i++) {
          out[i] = start[i];
          hits_so_far += out[i];
      }

      int starting_hits_so_far = hits_so_far;

      for (int t = 0; t < trials; t++) {
          hits_so_far = starting_hits_so_far;

          double p_hit = 0;
          double trial = 0;

          for (int i = nstart; i < nthrows; i++) {
              p_hit = (double)hits_so_far / (double)i;
              trial = (double)rand() / RAND_MAX;

              if (trial < p_hit){
                  out[i] = 1;
                  hits_so_far++;
              } else {
                  out[i] = 0;
              }
          }

          if (out[nthrows - 2] == 1) {
              trials_so_far++;
              final_throw_successes_so_far += out[nthrows - 1];
              success_rate = (double)final_throw_successes_so_far / 
                  (double)trials_so_far;
              if (trials_so_far % 100000 == 0){
                  Rcpp::Rcout << success_rate << "\\n";
              }
          }
      }
      return success_rate;
  }
')