 # Untitled

Sep 6th, 2018
100
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #! /usr/bin/python
2. # -*- coding: utf-8 -*-
3.
4. import copy as CO
5. import bisect as BI
6. import math as MA
7. import matplotlib.pyplot as PLT
8. '''
9. dice game
11. this one add matplotlib for visualizing the distri
12. '''
13.
14. class CountingValueError(ValueError):
15.     pass
16.
17. def arangements_ie( n, l ):
18.     """ arrangements_ie ( n,( k1, k2, k3,...) )
19.    Return the number of possible arrangements of n elements
20.    with k1, k2,.... elements being indistinguishable.
21.    """
22.     a= MA.factorial( n )
23.     K=0
24.     for k in l:
25.         K += k
26.         if K > n:
27.             raise CountingValueError ("Can't have k indistinguishable objects taken from n objects if k > n")
28.         a= a / MA.factorial( k )
29.     return a
30. def result(result):
31.     print "gain\t\tresult"
32.     for r in result:
33.         pass
34.
35.
36. def storegain(trialgain,trialproba):
37.     #print "### trial gain:",trialgain,"|proba:",trialproba*100,"%"
38.     global tp
39.     global gains
40.     gains.append(trialgain)
41.     gains.append(trialproba)
42.     tp+=trialproba
43.     i = BI.bisect_right(totwin,trialgain) -1
44.     totproba [i] += trialproba
45.
46. def trial(t):
47.     global g_no
48.     global g_N
49.     trialproba = 0.0
50.     trialgain = 0.0
51.     for i  in range (g_no):
52.         k=t[i]
53.         if k!=0:
54.             p = gamewin[i]
55.             gain = gamewin[i]
56.             trialgain += gain * k
57.             proba=MA.pow(p,k)
58.             if trialproba == 0:
59.                 trialproba= proba
60.             else:
61.                trialproba*= proba
62.             #print "-- outcome:",i,"|k:",k,"|p:",p,"|totgain:", gain * k
63.     trialproba*=arangements_ie(g_N,t)
64.
65.     storegain(trialgain, trialproba)
66.
67.
68. def callback ( t):
69.
70.     #print t
71.     #print "."
72.     global cpt
73.     cpt+=1
74.     trial(t)
75.
76. def rec_gen ( co , t ,s ):
77.     '''
78.     N: number of experiment still  available
79.     co: current outcome (position in t )
80.     t : list of occurence number for each outcome
81.    '''
82.
83.     for i in xrange ( 0,g_N-s+1 ):
84.         t[ co ] = i
85.         ss=s+i
86.         if ss == g_N :
87.              callback ( t )
88.              continue
89.
90.         if  co > 0 :
91.             t2 = CO.copy ( t )
92.             rec_gen ( co-1 , t2,ss )
93.
94.     del(t)
95.     return
96.
97. def gen ( no , N ):
98.     '''
99.    no: number outcome
100.    N : Number of subsequent experiment (trial of N experiment)
101.    t : contain the number of occurence for each outcomes
102.    '''
103.     global g_no  # number of outcome
104.     g_no = no
105.     global g_N  # number of experiment
106.     g_N = N
107.     print "===", g_no,g_N
108.     t = [ 0 ] * no
109.     rec_gen ( no-1 , t,0 )
110.
111. def plot1(totwin,totproba):
112.
113.     PLT.title("distri")
114.     PLT.xlabel("win")
115.     PLT.ylabel("proba")
116.
117.     PLT.plot(totwin,totproba)
118.
119.     PLT.grid(True)
120.     PLT.show()
121.
122.
123.
124.
125.
126.
127. if __name__ == "__main__":
128.
129.     # number of expiriment (how many time I throw the dices)
130.     N = 100
131.
132.     #list of outcomes  (proba, gain)
133.     gamewin = [ [1.0/6, 200], [2.0/6, 10], [3.0/6, 0] ]
134.
135.     # game expectation E(X)
136.     gameexpect= 0
137.     for g in gamewin:
138.         gameexpect += g[ 0 ] * g[ 1 ]
139.     print "Game expectation = ", gameexpect
140.     # game variance var(X)=E(X²)-E²(X)
141.     gamevar = 0
142.     for g in gamewin:
143.         gamevar +=  MA.pow( g[ 1 ], 2 ) * g[ 0 ]
144.     gamevar -= MA.pow( gameexpect , 2 )
145.     print "Game variance = ", gamevar
146.     print "Game Standard Error =", MA.sqrt(gamevar)
147.
148.     # maxwin after N experiments
149.     maxwin = 200 * N
150.
151.     print "Maxwin after N throws:", maxwin
152.     print "Expectation after  N throws:", gameexpect * N
153.
154.     # to have the total gains  we get after N experiments  well distributed in
155.     #  the result  list [(gains,proba), .... ]
156.     totwin = range( 0, int(maxwin), int(maxwin/500))
157.     totproba = [0.0] * len( totwin )
158.
159.     cpt = 0
160.     tp=0
161.     gains=[[],[]]
162.
163.     gen( len(gamewin),N )
164.     print "cpt",cpt,"tp",tp
165.     result = zip( totproba, totwin )
166.
167.     # verifying that the calculus are corrects:
168.     # total probability for all final gains = 1
169.     print "verify the sum of all probabilty = 1:",sum(totproba)
170.
171.     # game expectation E(Xn) after N throws
172.     totgameexpect= 0
173.     for g in result:
174.         totgameexpect += g[ 0 ] * g[ 1 ]
175.     print "Game expectation after N throws= ", totgameexpect
176.     # game variance after N throws var(Xn)=E(Xn²)-E²(Xn)
177.     totgamevar = 0
178.     for g in result:
179.         totgamevar +=  MA.pow( g[ 1 ], 2 ) * g[ 0 ]
180.     totgamevar -= MA.pow( totgameexpect , 2 )
181.     print "Game variance N throws = ", totgamevar
182.     print "Game Standard Error N throws =", MA.sqrt(totgamevar)
183.
184.
185.     #plot1( totwin, totproba )
186.     plot1( gains, gains )