1#  _______________________________________________________________________
2#
3#  PECOS: Parallel Environment for Creation Of Stochastics
4#  Copyright (c) 2011, Sandia National Laboratories.
5#  This software is distributed under the GNU Lesser General Public License.
6#  For more information, see the README file in the top Pecos directory.
7#  _______________________________________________________________________
8
9#!/projects/vuq/anaconda/bin/python
10import math
11
12p = 0.10
13N = 20
14
15def fval(x, order):
16  global p, N
17
18  if order == 0:
19    return 1.0
20
21  elif order == 1:
22    return (p*N-x)/(p*N)
23
24  elif order == 2:
25    return (p*p*N*(1.0-N)+(1.0-2.0*p+2.0*p*N)*x-x*x)/(p*p*N*(1.0-N))
26
27  elif order > 2:
28    om1 = order-1
29    fm2 = fval(x, order-2)
30    fm1 = fval(x, order-1)
31    return ((p*(N-om1)+om1*(1.0-p)-x)*fm1 - om1*(1.0-p)*fm2)/(p*(N-om1))
32
33  else:
34    raise("Bad order")
35
36
37def nCrBad(n,r):
38  f = math.factorial
39  return f(n) / f(r) / f(n-r)
40
41
42def nCr(n,k):
43  combination = 1.
44  nmk = n - k;
45  if k <= nmk:
46    for i in range(k):
47      combination *= float(n-i)/float(k-i);
48  else:
49    for i in range(nmk):
50      combination *= float(n-i)/float(nmk-i);
51
52  return combination
53
54
55def poch(r,n):
56  if n == 0:
57    return 1.0
58  else:
59    poch = r
60    for i in range(1,n):
61      poch *= r+i;
62    return poch
63
64
65def exact(order):
66  global p, N
67  return math.pow(-1,order)*math.factorial(order)/poch(-N,order)*math.pow((1.0-p)/p,order)
68
69
70def doSum(n, o1, o2):
71  global p, N
72  sum = 0.0
73  for x in range(0,n):
74    sum = sum + nCr(N,x)*p**x*(1.0-p)**(N-x)*fval(x,o1)*fval(x,o2)
75  return sum
76
77
78
79def main():
80  global N
81  print fval(1.23, 5)
82  #print poch(-2.0,1)
83  #print nCr(6,3)
84  order = 3
85  for i in range(0,N+1):
86    sum = doSum(i,order,order)
87    print i, sum, (sum - exact(order))/exact(order)*100.0
88  print 'Exact = '+str(exact(order))
89
90
91main()
92