CoCalc Shared FilesSimulations.sagews
Author: Paul Zeitz
Views : 26
#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!
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


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

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)

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']
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
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)




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
#################################
#
#
###################################

#    First, create a "first head" random variable, which simulates flipping a coin until it is heads,
#    and records which flip is the 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):

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):
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.

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):
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.