SharedSimulations.sagewsOpen in CoCalc
Author: Paul Zeitz
Views : 6
#random number random()
0.4430272557901488
#random integer randint(9,13)
11
#make a table of random stuff for _ in range(10): print randint(1,6)
#make a list of random stuff myDecimalData = [random() for _ in range(100)] myDice = [randint(1,6) for _ in range(10000)]
myDice
#do stats with lists mean(myDice)
8833/2500
#get MAD! def MAD(myList): m = mean(myList) absDevs=[abs(x-m) for x in myList] return mean(absDevs)
MAD(myDice)
1767/1250
histogram(myDice, bins=[0..7],align='left')
[myDice.count(k) for k in [1..6]]
[1632, 1636, 1660, 1687, 1688, 1697]
#sums of dice sumDice = [randint(1,6)+randint(1,6) for _ in range(100000)]
histogram(sumDice,bins=[0..13],align='left')
sumsOfMany = [tossDice(10) for _ in range(10000)] histogram(sumsOfMany,bins=[0..61],align='left')
#sums of many dice def tossDice(nTosses): total = 0 for _ in range(nTosses): total += randint(1,6) #increment total return total
tossDice(500)
1756
#counting heads #simulate tossing a coin several times and counting heads #simulate a coin with a random number that is either 0 or 1 #heads = 1, tails = 0 def count_heads(nFlips): outcomes=[randint(0,1) for _ in range(nFlips)] return sum(outcomes)
count_heads(100000)
49890
#Now do many trials nFlips = 1000 nTrials = 10000 coinData = [count_heads(nFlips) for _ in range(nTrials)] histogram(coinData,bins=[0..nFlips+1],align='left',xmin=-1) print float(MAD(coinData)/nFlips)
0.01245004704
#check how data deviates from mean sum([coinData.count(k) for k in [4..6]])
Error in lines 2-2 Traceback (most recent call last): File "/projects/dc288336-f2cd-4333-8eaf-56382a3bcf40/.sagemathcloud/sage_server.py", line 879, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'coinData' is not defined
#examine coin statistics for streaks #count maximum streak in list def streak(myList): start = myList[0] streakLength = 1 maxStreakLength = 1 for a in myList[1:]: if a == start: streakLength += 1 if streakLength>maxStreakLength: maxStreakLength = streakLength else: start = a streakLength = 1 return maxStreakLength #fair coin nTosses = 100 nTrials = 1000 #1=H, 0=T #outcome = [randint(0,1) for _ in range(nTosses)] streaks = [streak([randint(0,1) for _ in range(nTosses)]) for _ in range(nTrials)] #outcomes = [[randint(0,1) for _ in range(nTosses)] for _ in range(nTrials)] #for r in range(nTrials): #rOutcome = outcomes[r] #print rOutcome, streak(rOutcome) meanS = float(mean(streaks)) minS = min(streaks) maxS = max(streaks) print('%s trials of flipping a coin %s times'%(nTrials,nTosses)) histogram(streaks,bins=range(minS-1,maxS+2),align='left',alpha=0.8,color='green') print('Minimum streak length is %s, maximum is %s, mean is %s'%(minS,maxS,meanS))
#simulate selection of card from deck #and check whether there are facecards or not vals=['A','2','3','4','5','6','7','8','9','10','J','Q','K'] suits = ['\color{red}\diamondsuit','\clubsuit','\color{red}\heartsuit','\spadesuit'] nVals = len(vals) nSuits = len(suits) nCards = nVals*nSuits nTrials = 10 nDealt = 3 #number of cards to deal nFacecards = 0 #number of trials that have at least one facecard for _ in range(nTrials): cards = sample(range(nCards),nDealt) facecards = false for c in cards: s = floor(c/nVals) v = c- s*nVals if v > 9: facecards = true html('%s$%s$&nbsp;&nbsp;&nbsp;'%(vals[v],suits[s]),hide=False) print facecards if facecards: nFacecards += 1 print 'Out of %s trials, %s had face cards.'%(nTrials,nFacecards)
2\color{red}\diamondsuit   3\clubsuit   J\color{red}\diamondsuit   
True
9\spadesuit   10\color{red}\heartsuit   10\color{red}\diamondsuit   
False
Q\color{red}\diamondsuit   10\color{red}\heartsuit   3\spadesuit   
True
J\clubsuit   6\color{red}\heartsuit   4\color{red}\heartsuit   
True
J\color{red}\diamondsuit   2\color{red}\heartsuit   K\color{red}\heartsuit   
True
4\color{red}\heartsuit   3\spadesuit   A\color{red}\heartsuit   
False
A\color{red}\heartsuit   7\clubsuit   7\spadesuit   
False
8\clubsuit   3\color{red}\diamondsuit   6\spadesuit   
False
J\clubsuit   9\color{red}\diamondsuit   9\spadesuit   
True
6\color{red}\diamondsuit   10\clubsuit   4\color{red}\heartsuit   
False Out of 10 trials, 5 had face cards.
#unfair coin def coin(pHeads): if random()<pHeads: return 1 else: return 0
#count maximum streak in list def streak(myList): start = myList[0] streakLength = 1 maxStreakLength = 1 for a in myList[1:]: if a == start: streakLength += 1 if streakLength>maxStreakLength: maxStreakLength = streakLength else: start = a streakLength = 1 return maxStreakLength @interact def _(nFlips=[10,20,100,1000],nTrials=[1000,2000]): streaks = [streak([randint(0,1) for _ in range(nFlips)]) for _ in range(nTrials)] meanS = float(mean(streaks)) minS = min(streaks) maxS = max(streaks) print('%s trials of flipping a coin %s times'%(nTrials,nFlips)) h=histogram(streaks,bins=range(minS-1,maxS+2),align='left',alpha=0.8,color='green') print('Minimum is %s, maximum is %s, mean is %s'%(minS,maxS,meanS)) show(h)
Interact: please open in CoCalc
def lottery(): prize = 1 pPrize = .25 pTicket = 0.25 spin = random() if spin < pPrize: return prize else: if spin <pPrize+pTicket: return lottery() else: return 0
lottery()
1
outcomes= [lottery() for _ in range(10000)]
float(mean(outcomes))
0.3341
################################# # # ST PETERSBURG PARADOX # ################################### # First, create a "first head" random variable, which simulates flipping a coin until it is heads, # and records which flip is the first head. def first_head(): coin=0 # 0 for tails, 1 for heads. Start by setting the coin to tails count=0 # Keep track of the number of flips before first head while coin == 0: # the two lines in this "while loop" will be executed as long as the coin value is 0 # but once it equals 1, the loop ends and we move to the final line of the function def. coin=randint(0,1) # set coin equal to 0 or 1 randomly count += 1 # increment the count return count # this is the count value when the coin first equals 1 (and thus, the while loop ends). # test a few outputs and compute the average payoff (which, theoretically, is INFINITE). sum=0 nTrials=10 for _ in range(nTrials): fh = first_head() print('First head at flip # %s. You win %s dollars!'%(fh, 2^fh)) sum += 2^fh print('Average payoff was %s.' %float(sum/nTrials))
First head at flip # 2. You win 4 dollars! First head at flip # 1. You win 2 dollars! First head at flip # 1. You win 2 dollars! First head at flip # 1. You win 2 dollars! First head at flip # 2. You win 4 dollars! First head at flip # 4. You win 16 dollars! First head at flip # 1. You win 2 dollars! First head at flip # 4. You win 16 dollars! First head at flip # 2. You win 4 dollars! First head at flip # 4. You win 16 dollars! Average payoff was 6.8.
################################################################################ # Here is one simulation of playing 1000 games, where the price to play each # game is 5 dollars. The result is a plot of the money on hand after each play, # using the list_plot function. ################################################################################ price = 15 #cost to play the game once nGames = 1000 accumulation = 0 #player's profit sequence = [] #keep track of the first_head() outputs daily = [] #accumlation on each day; to be built for _ in range(nGames): fh = first_head() sequence.append(fh) accumulation += 2^fh-price daily.append(accumulation) list_plot(daily,plotjoined = True)
######################################################################################### # This is an INTERACTIVE simulation which allows you to change the price per game # and the total number of games. It also prints out the maximum payout, and # when this happened, the final amount of money on hand, the minimum amount, and the maximum, # and the average payout (which, theoretically, is INFINITY). # # Note the use of the @interactive command, followed by a function definition (the # name of the function can be anything, but we choose the generic "_" since we will never use the # function's name. You can look up the use of Sage interactive online, or just study this example. # It's pretty flexible, and easier to use than Mathematica's Manipulate command, although less slick. ############################################################################################## # We include the already-defined first_head function, so that this is a stand-alone cell. def first_head(): coin=0 # 0 for tails, 1 for heads. Start by setting the coin to tails count=0 # Keep track of the number of flips before first head while coin == 0: # the two lines in this "while loop" will be executed as long as the coin value is 0 # but once it equals 1, the loop ends and we move to the final line of the function def. coin = randint(0,1) # set coin equal to 0 or 1 randomly count += 1 # increment the count return count # this is the count value when the coin first equals 1 (and thus, the while loop ends). ############################################################################################################ # interactive part starts here ############################################################################################################ @interact def _(price= [5,10,15,20,100], nGames = [100,1000,10000,100000,1000000]): # note that the selection boxes for price and nGames are defined in the function defn above. accumulation = 0 #player's profit sequence = [] #keep track of the first_head() outputs daily = [] #accumlation on each day; to be built for _ in range(nGames): fh = first_head() sequence.append(fh) accumulation += 2^fh-price daily.append(accumulation) # Compute average payout totalPayout = 0 for i in range(len(sequence)): totalPayout += 2^sequence[i] avg = float(totalPayout/nGames) # Compute longest run of tails, max min, etc. longestRun=max(sequence) longestRunIndex=sequence.index(longestRun) print('Luckiest: game # %s, first head at position %s, winning $%s.'%(longestRunIndex, longestRun, 2^longestRun)) print('Lowest: %s'%min(daily)) print('Highest: %s'%max(daily)) print('Final: %s'%daily[nGames-1]) print('Average payout: %s.'%avg) # Plot the money on hand. l = list_plot(daily,plotjoined = True,gridlines=True) show(l) #in interactive mode, the plot doesn't display with just the list_plot() command. #instead, we use the show() function.
Interact: please open in CoCalc