Advertisement
dechicom

Untitled

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