Advertisement
dechicom

Untitled

Sep 6th, 2018
150
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.61 KB | None | 0 0
  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
  10. copyright 2018 Marc Dechico
  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,trialproba])
  41.     tp+=trialproba
  42.     i = BI.bisect_right(totwin,trialgain) -1
  43.     totproba [i] += trialproba
  44.  
  45. def trial(t):
  46.     global g_no
  47.     global g_N
  48.     trialproba = 0.0
  49.     trialgain = 0.0
  50.     for i  in range (g_no):
  51.         k=t[i]
  52.         if k!=0:
  53.             p = gamewin[i][0]
  54.             gain = gamewin[i][1]
  55.             trialgain += gain * k
  56.             proba=MA.pow(p,k)
  57.             if trialproba == 0:
  58.                 trialproba= proba
  59.             else:
  60.                trialproba*= proba
  61.             #print "-- outcome:",i,"|k:",k,"|p:",p,"|totgain:", gain * k
  62.     trialproba*=arangements_ie(g_N,t)
  63.                
  64.     storegain(trialgain, trialproba)
  65.  
  66.  
  67. def callback ( t):
  68.  
  69.     #print t
  70.     #print "."
  71.     global cpt
  72.     cpt+=1
  73.     trial(t)
  74.  
  75. def rec_gen ( co , t ,s ):
  76.     '''
  77.     N: number of experiment still  available
  78.     co: current outcome (position in t )
  79.     t : list of occurence number for each outcome
  80.    '''
  81.    
  82.     for i in xrange ( 0,g_N-s+1 ):
  83.         t[ co ] = i
  84.         ss=s+i
  85.         if ss == g_N :
  86.              callback ( t )
  87.              continue
  88.  
  89.         if  co > 0 :
  90.             t2 = CO.copy ( t )
  91.             rec_gen ( co-1 , t2,ss )
  92.  
  93.     del(t)
  94.     return
  95.  
  96. def gen ( no , N ):
  97.     '''
  98.    no: number outcome  
  99.    N : Number of subsequent experiment (trial of N experiment)
  100.    t : contain the number of occurence for each outcomes
  101.    '''
  102.     global g_no  # number of outcome
  103.     g_no = no
  104.     global g_N  # number of experiment
  105.     g_N = N
  106.     print "===", g_no,g_N
  107.     t = [ 0 ] * no
  108.     rec_gen ( no-1 , t,0 )
  109.  
  110. def plot1(totwin,totproba):
  111.  
  112.     PLT.title("distri")
  113.     PLT.xlabel("win")
  114.     PLT.ylabel("proba")  
  115.    
  116.     PLT.plot(totwin,totproba)
  117.  
  118.     PLT.grid(True)
  119.     PLT.show()
  120.  
  121.  
  122.  
  123.  
  124.  
  125.    
  126. if __name__ == "__main__":
  127.  
  128.     # number of expiriment (how many time I throw the dices)
  129.     N = 100
  130.  
  131.     #list of outcomes  (proba, gain)
  132.     gamewin = [ [1.0/6, 200], [2.0/6, 10], [3.0/6, 0] ]
  133.  
  134.     # game expectation E(X)
  135.     gameexpect= 0
  136.     for g in gamewin:
  137.         gameexpect += g[ 0 ] * g[ 1 ]
  138.     print "Game expectation = ", gameexpect
  139.     # game variance var(X)=E(X²)-E²(X)
  140.     gamevar = 0
  141.     for g in gamewin:
  142.         gamevar +=  MA.pow( g[ 1 ], 2 ) * g[ 0 ]
  143.     gamevar -= MA.pow( gameexpect , 2 )
  144.     print "Game variance = ", gamevar
  145.     print "Game Standard Error =", MA.sqrt(gamevar)
  146.  
  147.     # maxwin after N experiments
  148.     maxwin = 200 * N
  149.    
  150.     print "Maxwin after N throws:", maxwin
  151.     print "Expectation after  N throws:", gameexpect * N
  152.  
  153.     # to have the total gains  we get after N experiments  well distributed in
  154.     #  the result  list [(gains,proba), .... ]  
  155.     totwin = range( 0, int(maxwin), int(maxwin/500))
  156.     totproba = [0.0] * len( totwin )
  157.  
  158.     cpt = 0
  159.     tp=0
  160.     gains=[]
  161.     gen( len(gamewin),N )
  162.     print "cpt",cpt,"tp",tp
  163.     result = zip( totproba, totwin )
  164.  
  165.     # verifying that the calculus are corrects:
  166.     # total probability for all final gains = 1
  167.     print "verify the sum of all probabilty = 1:",sum(totproba)
  168.  
  169.     # game expectation E(Xn) after N throws
  170.     totgameexpect= 0
  171.     for g in result:
  172.         totgameexpect += g[ 0 ] * g[ 1 ]
  173.     print "Game expectation after N throws= ", totgameexpect
  174.     # game variance after N throws var(Xn)=E(Xn²)-E²(Xn)
  175.     totgamevar = 0
  176.     for g in result:
  177.         totgamevar +=  MA.pow( g[ 1 ], 2 ) * g[ 0 ]
  178.     totgamevar -= MA.pow( totgameexpect , 2 )
  179.     print "Game variance N throws = ", totgamevar
  180.     print "Game Standard Error N throws =", MA.sqrt(totgamevar)
  181.  
  182.     # have                  
  183.     plot1(totwin,totproba)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement