#! /usr/bin/perl -w
use strict;

for my $n (50, 100, 200) {
  open(GNUPLOT, "| gnuplot > returns-$n.png")
    or die "Can't run gnuplot: $!";

  my @formulas = (
    "n = $n",
    "log_ret(m, x) = (m*log(1+x) + (n-m)*log(1-x))/n",
    "rate(x) = 100*(exp(x) - 1)",
    "0 lt -1 title ''",
  );
  
  for my $percentile (50, 25, 10, 5) {
    my $m = percentile_outcome($n, 0.6, $percentile);
    push @formulas
      , "rate(log_ret($m, x)) lw 2 title '${percentile}th Percentile'";
  }

  my $formula_string = join ", ", @formulas;

  print GNUPLOT qq{
    set terminal png
    set title "Average return after $n bets"
    set ylabel "Average % Return/Bet"
    set xlabel "Fraction of Stake Bet"
    plot [0:0.3] [-3:3] $formula_string
  };
}

sub log_factorial {
  my $n = shift;
  my $ans = 0;
  for my $i (1..$n) {
    $ans += log($i);
  }

  return $ans;
}

sub prob_outcome {
  my ($n, $m, $p) = @_;
  return exp(
      log_factorial($n)
    - log_factorial($m)
    + $m*log($p)
    - log_factorial($n-$m)
    + ($n-$m)*log(1-$p)
  );
}

sub percentile_outcome {
  my ($n, $p, $percentile) = @_;
  my $total_p = 0;
  for my $m (0..$n) {
    $total_p += prob_outcome($n, $m, $p);
    return $m if $percentile <= 100*$total_p;
  }
  return $n;
}