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