histo-1p.awk

#!/usr/bin/awk -f

function frac(a) { return a-int(a+1e-10) }
function mod(a,b) { return b*frac(a/b) }
function round(a,b) { return b*int(int(a/b+0.5)+1) }
function bound(x,a,b) { return x<a ? a : x>b?b:x }

#function getbin(x,max,n) { return int(bound(x*n/max,0,n)+offs) }
function getbin(x,max,n) { return x<0 ? -1 : x>max ? n+1 : int(x*n/max+offs) }

function putdat(x, max,n, bins, weight) {
  i=getbin(x,max,n)
#  print i,x,max,n
  bins["n"]+=weight
  bins[i,0]+=weight
}

function pline(i) { return i<19 ? $(i+5) : $(i+6) }
function bparty(i) { return pline(i-1+19) }

function fact(n) {
  if(lastn && n==lastn) { return lastfact }
  res = n<=0.5?1:n*fact(n-1)
  lastn=n
  lastfact=res
  return res
}

function pow(a,b) { return a>0?exp(b*log(a)):0 }

BEGIN {
  OFS="\t"; FS="\t"
  binn=200
  itern=2
  fail_p=0.017

  nd_target=3
  print "#", "itern=" itern, "nd_target=" nd_target, "fail_p=" fail_p, "binn=" binn
}

/^[^#U]/ {
  all=pline(3)+pline(4)+pline(5)
  normp=pline(1)+0
  normb=pline(2)+0
  er=bparty(6)
  sr=bparty(1)
  #if(normp==0) normp=1
  #turn=100*all/(normp+rand()-0.5)
  if(all==0) all=1
    for(iter=0; iter<itern; iter++) {
  erp=100*(er+rand()-0.5)/(all+rand()-0.5)

#  putdat(erp, 100,binn, bin, 1)
  if(pline(9)==nd_target) putdat(erp, 100,binn, bin, 1)
  putdat(erp, 100,binn, binmod, pow(fail_p*all,nd_target)/fact(nd_target)*exp(-fail_p*all))
} }

END {
  print "# total=" bin["n"] " modtotal=" binmod["n"]
  for(i=-1; i<=binn+1; i++) {
    print (i+0.5-offs)*100/binn, bin[i,0]/itern+0,  binmod[i,0]/itern+0
  }
}