[or-cvs] r15815: Verify we actually have a proper probability distribution vi (torflow/branches/gsoc2008/tools/BTAnalysis)

[mikeperry at seul.org]

Author: mikeperry
Date: 2008-07-10 03:01:35 -0400 (Thu, 10 Jul 2008)
New Revision: 15815

Modified:
   torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py
Log:

Verify we actually have a proper probability distribution via
some numerical integration. Yep, we do. GNUplot is just
insane.



Modified: torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py
===================================================================
--- torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py	2008-07-10 01:11:26 UTC (rev 15814)
+++ torflow/branches/gsoc2008/tools/BTAnalysis/shufflebt.py	2008-07-10 07:01:35 UTC (rev 15815)
@@ -10,6 +10,8 @@
 import getopt,sys,os
 import popen2
 import math,copy
+from scipy.integrate import *
+from numpy import trapz
 
 class Stats:
   def __init__(self,file):
@@ -137,10 +139,16 @@
 def pareto(k,Xm,N,fname):
   # gnuplot string for shifted, normalized exponential PDF
   # g(x,k,B) = (N * k*(Xm**k)/x**(k+1)))
-  ps = fname+'(x)=(x<='+str(Xm)+') ? 0 : (('+str((N*k)*(Xm**k))+')/(x-'+str(Xm)+ ')**('+str(k+1)+'))\n'
+  ps = fname+'(x)=(x<='+str(Xm)+') ? 0 : (('+str((N*k)*(Xm**k))+')/((x)**('+str(k+1)+')))\n'
   #ps = fname+'(x)='+str(N*k*(Xm**k))+'/x**('+str(k+1)+')\n'
   return ps
 
+def pypareto(x, k,Xm,N):
+  # gnuplot string for shifted, normalized exponential PDF
+  # g(x,k,B) = (N * k*(Xm**k)/x**(k+1)))
+  if x<Xm: return 0
+  else: return ((((N*k)*(Xm**k)))/((x)**((k+1))))
+
 def exp(mean,shift,N,fname):
   # gnuplot string for normalized exponential PDF
   # g(x,k,B) = N * l*exp(-l*(x-shift))
@@ -286,13 +294,17 @@
     parK = s.paretoK(mode)
     modeN = s.modeN(mode)
     modeMean = s.modeMean(mode)
+    # verify sanity by integrating scaled distribution:
+    modeNint = trapz(map(lambda x: pypareto(x, parK, mode, modeN),
+                     xrange(1,200000)))
+
     print 'Resolution of histogram:',res,'ms'
     print 'Mean: '+str(mean)+', mode: '+str(mode)
     print 'ParK: '+str(parK)
-    print 'ModeN',modeN 
+    print 'ModeN: '+str(modeN)+" vs integrated: "+str(modeNint)
     print '#successful runs:',len(s.values)
     # get stats
-  
+   
     if graph:
       # create gnuplot file
 #      if not sort and not truncate:
@@ -345,8 +357,8 @@
       f.close()
     
       # plot the file
-      #p = popen2.Popen4('gnuplot ' + plotname)
-      p = popen2.Popen4('gp4.2 ' + plotname) #peculiarity of fallon's system
+      p = popen2.Popen4('gnuplot ' + plotname)
+      #p = popen2.Popen4('gp4.2 ' + plotname) #peculiarity of fallon's system
     
       p.wait()
       for err in p.fromchild: