daily pastebin goal
53%
SHARE
TWEET

Untitled

a guest Sep 25th, 2018 66 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. # Let X be multivariate hypergeometrically distributed with s
  2. # categories having assigned values x_1,...x_s having f_1,...f_s
  3. # instances in the urn.  m is the sum of all the f_i.
  4.  
  5. # Let Y be the sum of of n observations from X. Y is discrete random
  6. # variable of s' values y_1,...y_s' with frequency g_1,...g_s'.
  7.  
  8. # Y will range from n*x_0 to n*x_s, and s' is given by
  9. # | Outer Sum from 1..n of the set { x_0,...,x_s } |
  10.  
  11. from __future__ import division
  12. import sys
  13. # Factorial of N
  14. def factorial( N ):
  15.     if N == 0:
  16.         return 1.0
  17.     else:
  18.         r = 1.0
  19.         for i in range( 2, N+1 ):
  20.             r *= i
  21.         return r
  22.  
  23. # Number of combinations of r items taken from n
  24. def combinations( n, r ):
  25.     nCr = 1
  26.     if (r*2) > n:
  27.         r = n-r
  28.     for i in range( r-1, -1, -1 ):
  29.         nCr = ( nCr * ( n - i ) ) // ( r - i )
  30.     return nCr
  31.  
  32. # Enumerate probability mass pieces from the sum-of-multivariate-hypergeometrics distribution
  33. # Y, g, s, m, n, f, x are described at the top
  34. # i is the number of the category from X that we are processing
  35. # prob[Y] is the probability of observing sum Y
  36. # free is the number of observations (out of n) that remain to be used
  37. # maxfree[i] is the maximum allowed value of free at point i
  38. # k is the number of observations to use from category i
  39. def enumerate_mhyper_sum( i, prob, free, maxfree, Y, g, s, m, n, f, x  ):
  40.     if free == 0:
  41.         if Y not in prob: prob[Y] = 0
  42.         prob[Y] += g
  43.     else:
  44.         for k in range( max( 0, free-maxfree[i+1] ), min( f[i], free ) + 1 ):
  45.             enumerate_mhyper_sum( i+1, prob, free-k, maxfree, Y+(x[i]*k), g*combinations(f[i],k), s, m, n, f, x )
  46.  
  47. # Construct the sum-of-multivariate-hypergeometrics distribution
  48. # Return a discrete distribution such that values[i] has a frequency
  49. # of mass[i]
  50. def construct_mhyper_sum( x, f, n ):
  51.     s = len(x)
  52.     m = sum(f)
  53.     prob = dict()
  54.     maxfree = list( f )
  55.     for i in range( s-2, -1, -1 ):
  56.         maxfree[i] += maxfree[i+1]
  57.     maxfree += [ 0 ]
  58.     enumerate_mhyper_sum( 0, prob, n, maxfree, 0, 1, s, m, n, f, x )
  59.     values = []
  60.     mass = []
  61.     for k in sorted( prob.keys() ):
  62.         values.append( k )
  63.         mass.append( prob[k] / combinations(m,n) )
  64.     return ( values, mass )
  65.  
  66. n = int(sys.argv[4])
  67. x = [ -1, 0, 1 ]
  68. f = [ int(sys.argv[1]),int(sys.argv[2]) , int(sys.argv[3]) ]
  69. t = int(sys.argv[5])
  70.  
  71. #print f
  72. (values,mass) = construct_mhyper_sum( x, f, n )
  73.  
  74. #for (v,d) in zip(values,mass):
  75. #   print v, d
  76.  
  77. #print "Probability Mass Sums to: %f" % sum(mass)
  78. #print "Population Mean (expected sum of %d observations) is: %f" % (n, sum( v*d for (v,d) in zip(values,mass) ) )
  79. pval = 0
  80. #print t #debug
  81. for num in range(len(values)):
  82.     if t == values[num]:
  83.         if t >= 0:
  84.             for calc in range( num, len(values) ):
  85.             #       print mass[calc]        #debug
  86.                 pval = pval + mass[calc]
  87.         else:
  88.             for calc in range( 0, num+1 ):
  89.                 pval = pval + mass[calc]
  90.                 #print mass[calc]           #debug
  91.  
  92. print "1-tailed: ", pval
  93. print "2-tailed: ", pval*2
  94. #from matplotlib import pylab as p
  95. #p.plot( vals, density )
  96. #input()
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top