Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #! /usr/bin/python
- # -*- coding: utf-8 -*-
- import copy as CO
- import bisect as BI
- import math as MA
- import matplotlib.pyplot as PLT
- '''
- dice game
- copyright 2018 Marc Dechico
- this one add matplotlib for visualizing the distri
- '''
- class CountingValueError(ValueError):
- pass
- def arangements_ie( n, l ):
- """ arrangements_ie ( n,( k1, k2, k3,...) )
- Return the number of possible arrangements of n elements
- with k1, k2,.... elements being indistinguishable.
- """
- a= MA.factorial( n )
- K=0
- for k in l:
- K += k
- if K > n:
- raise CountingValueError ("Can't have k indistinguishable objects taken from n objects if k > n")
- a= a / MA.factorial( k )
- return a
- def result(result):
- print "gain\t\tresult"
- for r in result:
- pass
- def storegain(trialgain,trialproba):
- #print "### trial gain:",trialgain,"|proba:",trialproba*100,"%"
- global tp
- global gains
- gains[0].append(trialgain)
- gains[1].append(trialproba)
- tp+=trialproba
- i = BI.bisect_right(totwin,trialgain) -1
- totproba [i] += trialproba
- def trial(t):
- global g_no
- global g_N
- trialproba = 0.0
- trialgain = 0.0
- for i in range (g_no):
- k=t[i]
- if k!=0:
- p = gamewin[i][0]
- gain = gamewin[i][1]
- trialgain += gain * k
- proba=MA.pow(p,k)
- if trialproba == 0:
- trialproba= proba
- else:
- trialproba*= proba
- #print "-- outcome:",i,"|k:",k,"|p:",p,"|totgain:", gain * k
- trialproba*=arangements_ie(g_N,t)
- storegain(trialgain, trialproba)
- def callback ( t):
- #print t
- #print "."
- global cpt
- cpt+=1
- trial(t)
- def rec_gen ( co , t ,s ):
- '''
- N: number of experiment still available
- co: current outcome (position in t )
- t : list of occurence number for each outcome
- '''
- for i in xrange ( 0,g_N-s+1 ):
- t[ co ] = i
- ss=s+i
- if ss == g_N :
- callback ( t )
- continue
- if co > 0 :
- t2 = CO.copy ( t )
- rec_gen ( co-1 , t2,ss )
- del(t)
- return
- def gen ( no , N ):
- '''
- no: number outcome
- N : Number of subsequent experiment (trial of N experiment)
- t : contain the number of occurence for each outcomes
- '''
- global g_no # number of outcome
- g_no = no
- global g_N # number of experiment
- g_N = N
- print "===", g_no,g_N
- t = [ 0 ] * no
- rec_gen ( no-1 , t,0 )
- def plot1(totwin,totproba):
- PLT.title("distri")
- PLT.xlabel("win")
- PLT.ylabel("proba")
- PLT.plot(totwin,totproba)
- PLT.grid(True)
- PLT.show()
- if __name__ == "__main__":
- # number of expiriment (how many time I throw the dices)
- N = 100
- #list of outcomes (proba, gain)
- gamewin = [ [1.0/6, 200], [2.0/6, 10], [3.0/6, 0] ]
- # game expectation E(X)
- gameexpect= 0
- for g in gamewin:
- gameexpect += g[ 0 ] * g[ 1 ]
- print "Game expectation = ", gameexpect
- # game variance var(X)=E(X²)-E²(X)
- gamevar = 0
- for g in gamewin:
- gamevar += MA.pow( g[ 1 ], 2 ) * g[ 0 ]
- gamevar -= MA.pow( gameexpect , 2 )
- print "Game variance = ", gamevar
- print "Game Standard Error =", MA.sqrt(gamevar)
- # maxwin after N experiments
- maxwin = 200 * N
- print "Maxwin after N throws:", maxwin
- print "Expectation after N throws:", gameexpect * N
- # to have the total gains we get after N experiments well distributed in
- # the result list [(gains,proba), .... ]
- totwin = range( 0, int(maxwin), int(maxwin/500))
- totproba = [0.0] * len( totwin )
- cpt = 0
- tp=0
- gains=[[],[]]
- gen( len(gamewin),N )
- print "cpt",cpt,"tp",tp
- result = zip( totproba, totwin )
- # verifying that the calculus are corrects:
- # total probability for all final gains = 1
- print "verify the sum of all probabilty = 1:",sum(totproba)
- # game expectation E(Xn) after N throws
- totgameexpect= 0
- for g in result:
- totgameexpect += g[ 0 ] * g[ 1 ]
- print "Game expectation after N throws= ", totgameexpect
- # game variance after N throws var(Xn)=E(Xn²)-E²(Xn)
- totgamevar = 0
- for g in result:
- totgamevar += MA.pow( g[ 1 ], 2 ) * g[ 0 ]
- totgamevar -= MA.pow( totgameexpect , 2 )
- print "Game variance N throws = ", totgamevar
- print "Game Standard Error N throws =", MA.sqrt(totgamevar)
- #plot1( totwin, totproba )
- plot1( gains[0], gains[1] )
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement