• API
• FAQ
• Tools
• Archive
daily pastebin goal
13%
SHARE
TWEET

# Untitled

a guest Sep 25th, 2018 72 Never
ENDING IN00days00hours00mins00secs

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.

Top