Calculation of $\pi$ with Monte Carlo algorithm using direct sampling and the calculation of rms deviation.
%pylab inline
import random, math
def direct_pi(N):
n_hits = 0
for i in range(N):
x, y = random.uniform(-1.0, 1.0), random.uniform(-1.0, 1.0)
if x ** 2 + y ** 2 < 1.0:
n_hits += 1
return n_hits
n_runs = 10
n_trials = 4000
for run in range(n_runs):
print 4.0 * direct_pi(n_trials) / float(n_trials)
#Calculation of rms deviation for 500 runs:
n_runs = 100
n_trials_list = []
sigmasqs = []
for poweroftwo in range(4, 13):
n_trials = 2 ** poweroftwo
sigmasq = 0.0
for run in range(n_runs):
pi_est = 4.0 * direct_pi(n_trials) / float(n_trials)
sigmasq += (pi_est - math.pi) ** 2
sigmasqs.append(math.sqrt(sigmasq/(n_runs)))
n_trials_list.append(n_trials)
pylab.plot(n_trials_list, sigmasqs, 'o')
pylab.plot([10.0, 10000.0], [1.642 / math.sqrt(10.0), 1.642 / math.sqrt(10000.0)])
pylab.xscale('log')
pylab.yscale('log')
pylab.xlabel('number of trials')
pylab.ylabel('root mean square deviation')
pylab.title('Direct sampling of pi: root mean square deviation vs. n_trials')
pylab.savefig('direct_sampling_rms_deviation.png')
pylab.show()
Populating the interactive namespace from numpy and matplotlib 3.138 3.153 3.139 3.152 3.168 3.115 3.16 3.176 3.185 3.161
Bunching method to calculate the error without knowing the value of $\pi$.
import random, math
n_trials = 400000
n_hits = 0
var = 0.0
var_est = 0.0
Obs_exp = 0.0
Obs2_exp = 0.0
for iter in range(n_trials):
x, y = random.uniform(-1.0, 1.0), random.uniform(-1.0, 1.0)
Obs = 0.0
if x**2 + y**2 < 1.0:
n_hits += 1
Obs = 4.0
Obs_exp+=Obs/n_trials #expectation value of the observable
Obs2_exp+=(Obs)**2/n_trials #expectation value of the square of the observable
var+=(Obs-math.pi)**2/n_trials #calculation of the variance if the value of pi had been known beforehand
var_est = Obs2_exp - Obs_exp**2 #best estimate to the variance
print 4.0 * n_hits / float(n_trials), math.sqrt(var), math.sqrt(var_est)
3.14122 1.64244240401 1.64244236172
Calculation of $\pi$ with Monte Carlo algorithm using Markov chain sampling, multiple runs.
import random
def markov_pi(N, delta):
x, y = 1.0, 1.0
n_hits = 0
n_accept = 0
for i in range(N):
del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)
if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:
n_accept += 1
x, y = x + del_x, y + del_y
if x**2 + y**2 < 1.0: n_hits += 1
return n_hits
n_runs = 10
n_trials = 4000
delta = 0.1
for run in range(n_runs):
print 4.0 * markov_pi(n_trials, delta) / float(n_trials)
3.109 3.178 2.955 3.107 3.225 3.161 2.85 3.397 2.841 2.97
import random
def markov_pi_acceptance(N, delta):
x, y = 1.0, 1.0
n_hits = 0
n_accept = 0
for i in range(N):
del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)
if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:
n_accept += 1
x, y = x + del_x, y + del_y
if x**2 + y**2 < 1.0: n_hits += 1
return n_accept
n_runs = 1
n_trials = 2**12
for delta in [0.062,0.125,0.25,0.5,1.0,2.0,4.0]:
print markov_pi(n_trials, delta)/ float(n_trials)
0.970458984375 0.96240234375 0.873046875 0.747314453125 0.548583984375 0.24609375 0.068359375
%pylab inline
import random, math, pylab
def markov_pi(N, delta):
x, y = 1.0, 1.0
n_hits = 0
for i in range(N):
del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)
if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:
x, y = x + del_x, y + del_y
if x**2 + y**2 < 1.0: n_hits += 1
return n_hits
n_runs = 500
for delta in [0.062, 0.125, 0.25, 0.5, 1.0, 2.0, 4.0]:
n_trials_list = []
sigmas = []
for poweroftwo in range(4, 13):
n_trials = 2 ** poweroftwo
sigma = 0.0
for run in range(n_runs):
pi_est = 4.0 * markov_pi(n_trials, delta) / float(n_trials)
sigma += (pi_est - math.pi) ** 2
sigmas.append(math.sqrt(sigma/(n_runs)))
n_trials_list.append(n_trials)
pylab.plot(n_trials_list, sigmas, 'o', ms = 8, label = '$\delta = $' + str(delta))
pylab.xscale('log')
pylab.yscale('log')
pylab.xlabel('number of trials')
pylab.ylabel('root mean square deviation')
pylab.plot([10,10000],[1.642 / math.sqrt(10.0), 1.642 / math.sqrt(10000.0)], label = 'direct')
pylab.title('Markov-chain sampling of pi: root mean square deviation vs. n_trials')
pylab.legend(loc='upper right')
pylab.savefig('markov_sampling_rms_deviation.png')
pylab.show()
Populating the interactive namespace from numpy and matplotlib
Bunching method for Markov-chain sampling. By this method, Markov-chain data are analyzed as direct-sampling ones and the error is underestimated at the beginning since the data is correlated, but that at later stages of the bunching they become more and more uncorrelated, just as direct-sampling data and the saturation (plateau) comes from the fact that the data are almost unrelated. For uncorrelated data, one has a plateau.
%pylab inline
import random, pylab, math
def markov_pi_all_data(N, delta):
x, y = 1.0, 1.0
data = []
for i in range(N):
del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)
if abs(x + del_x) < 1.0 and abs(y + del_y) < 1.0:
x, y = x + del_x, y + del_y
if x ** 2 + y ** 2 < 1.0:
data.append(4.0)
else:
data.append(0.0)
return data
poweroftwo = 20
n_trials = 2 ** poweroftwo
delta = 0.1
data = markov_pi_all_data(n_trials, delta)
errors = []
bunches = []
for i in range(poweroftwo):
new_data = []
mean = 0.0
mean_sq = 0.0
N = len(data)
while data != []:
x = data.pop()
y = data.pop()
mean += x + y
mean_sq += x ** 2 + y ** 2
new_data.append((x + y) / 2.0 )
errors.append(math.sqrt(mean_sq / N - (mean / N) ** 2) / math.sqrt(N))
bunches.append(i)
data = new_data[:]
print mean / float(N), 'mean value, estimate of pi'
print mean / float(N) - pi, 'mean value, estimate of pi'
pylab.plot(bunches, errors, 'o')
pylab.xlabel('iteration')
pylab.ylabel('apparent error')
pylab.title('Bunching: naive error vs iteration number')
pylab.savefig('apparent_error_bunching.png')
pylab.show()
Populating the interactive namespace from numpy and matplotlib 3.15369415283 mean value, estimate of pi 0.0121014992422 mean value, estimate of pi
Simulation of the 3$\times$3 pebble game with animation.
%pylab qt
import random, pylab
sigma = 0.4 # sigma and s_map are needed for the graphical output
s_map = [(1.0, 1.0), (2.0, 1.0), (3.0, 1.0),
(1.0, 2.0), (2.0, 2.0), (3.0, 2.0),
(1.0, 3.0), (2.0, 3.0), (3.0, 3.0)]
neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],
[4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],
[7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]
site = 8
N_runs = 50
for run in range(N_runs):
if run < 10: number_string = '0'+str(run)
else: number_string = str(run)
# Begin of graphical output
pylab.clf() #clears the previous figure
cir = pylab.Circle(s_map[site], radius=sigma, fc='r')
pylab.gca().add_patch(cir)
pylab.plot([0.5, 3.5], [1.5, 1.5], 'b')
pylab.plot([0.5, 3.5], [2.5, 2.5], 'b')
pylab.plot([1.5, 1.5], [0.5, 3.5], 'b')
pylab.plot([2.5, 2.5], [0.5, 3.5], 'b')
pylab.title('t = '+ number_string)
pylab.axis('scaled')
pylab.axis([0.5, 3.5, 0.5, 3.5])
pylab.xticks([])
pylab.yticks([])
pylab.pause(0.05) #Pause to generate a real time animation
pylab.show()
# End of graphical output
site = neighbor[site][ random.randint(0, 3)]
Populating the interactive namespace from numpy and matplotlib
Basic pebble throw game with homogeneous steady state probabilities. Multiple runs. The probability of landing on each grid element after a specified number of runs is visualised by a two dimensional histogram.
import random
#define the grid in terms of transitions:[0=right,1=up,2=left,3=down]
neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],
[4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],
[7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]
t_max = 4 #maximum number of pebble throws at each run
N_runs = 10 #number of runs
for run in range(N_runs): #run many times
site = 8 #initial site
t = 0 #initialise time
while t < t_max: #a pebble game of t_max throws
t += 1
site = neighbor[site][random.randint(0, 3)] #Pick a random transition from the initial site
print site #prints the final site after t_max throws in each run
7 7 7 7 4 2 7 6 4 6
%matplotlib qt
import random
import numpy as np
import matplotlib.pyplot as plt
xvec = {1:3, 2:2, 3:1, 4:3, 5:2, 6:1, 7:3, 8:2, 9:1}
yvec = {1:1, 2:1, 3:1, 4:2, 5:2, 6:2, 7:3, 8:3, 9:3}
neighbor = {1 : [2, 4, 1, 1], 2 : [3, 5, 1, 2], 3 : [3, 6, 2, 3],
4 : [5, 7, 4, 1], 5 : [6, 8, 4, 2], 6 : [6, 9, 5, 3],
7 : [8, 7, 7, 4], 8 : [9, 8, 7, 5], 9 : [9, 9, 8, 6]}
N_runs = 10
for run in range(N_runs):
list_vec = []
if run < 10: run_str= '0'+str(run)
else: run_str = str(run)
for n_runs in range(100000):
pos = 9
for iter in range(run):
pos = neighbor[pos][ random.randint(0, 3)]
list_vec.append(pos)
x = [xvec[k] for k in list_vec]
y = [yvec[k] for k in list_vec]
plt.xticks([])
plt.yticks([])
H, xedges, yedges = np.histogram2d(x, y, bins=(3, 3),
range=[[1,3],[1,3]], normed=True)
print H
H /= np.sum(H)
print H
extent = [yedges[0], yedges[-1], xedges[-1], xedges[0]]
plt.clf() #clears the previous figure
histo = plt.imshow(H, extent=extent, interpolation='nearest', vmin=0, vmax=1.00)
histo.set_cmap('hot')
plt.colorbar()
plt.title('t = '+str(run),fontsize=22)
plt.savefig('3x3_pebble_run_'+run_str+'.png')
plt.pause(0.05)
plt.show()
[[0. 0. 2.25] [0. 0. 0. ] [0. 0. 0. ]] [[0. 0. 1.] [0. 0. 0.] [0. 0. 0.]] [[0. 0.5585625 1.1299725] [0. 0. 0.561465 ] [0. 0. 0. ]] [[0. 0.24825 0.50221] [0. 0. 0.24954] [0. 0. 0. ]] [[0.137655 0.4276125 0.841815 ] [0. 0.2833425 0.419985 ] [0. 0. 0.13959 ]] [[0.06118 0.19005 0.37414] [0. 0.12593 0.18666] [0. 0. 0.06204]] [[0.1775925 0.4267125 0.6258825] [0.10782 0.212625 0.4200525] [0. 0.1052325 0.1740825]] [[0.07893 0.18965 0.27817] [0.04792 0.0945 0.18669] [0. 0.04677 0.07737]] [[0.2203875 0.358515 0.5262525] [0.1221075 0.26874 0.360585 ] [0.052425 0.1224675 0.21852 ]] [[0.09795 0.15934 0.23389] [0.05427 0.11944 0.16026] [0.0233 0.05443 0.09712]] [[0.2317275 0.3440025 0.44424 ] [0.1649925 0.2417625 0.33723 ] [0.0882 0.1676925 0.2301525]] [[0.10299 0.15289 0.19744] [0.07333 0.10745 0.14988] [0.0392 0.07453 0.10229]] [[0.2418075 0.311715 0.39024 ] [0.18441 0.2544525 0.31446 ] [0.12519 0.18315 0.244575 ]] [[0.10747 0.13854 0.17344] [0.08196 0.11309 0.13976] [0.05564 0.0814 0.1087 ]] [[0.2434725 0.3036825 0.356175 ] [0.19872 0.2446875 0.304065 ] [0.152775 0.20106 0.2453625]] [[0.10821 0.13497 0.1583 ] [0.08832 0.10875 0.13514] [0.0679 0.08936 0.10905]] [[0.2458575 0.287865 0.3297825] [0.20898 0.248625 0.2894625] [0.18018 0.2109375 0.24831 ]] [[0.10927 0.12794 0.14657] [0.09288 0.1105 0.12865] [0.08008 0.09375 0.11036]] [[0.2522925 0.2791125 0.2991825] [0.2209725 0.24426 0.28458 ] [0.1943325 0.2231775 0.25209 ]] [[0.11213 0.12405 0.13297] [0.09821 0.10856 0.12648] [0.08637 0.09919 0.11204]]
Exact simulation of the 3$\times$3 pebble game can be done using the transfer matrix for 50 throws. The evolution of the probability vector indicates that at large times the effect of the multiplication by the transfer matrix on the probability vector is negligible compared to the changes at small times. That means that the equilibrium probability vector is an eigenvector of the transfer matrix with eigenvalue 1.
import numpy
neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],
[4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],
[7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]
transfer = numpy.zeros((9, 9)) #initialise the transfer matrix
#the transfer matrix is constructed by adding the probabilities for all of the the terms
#that contribute to a specific transition to determine each element in the matrix
for k in range(9):
for neigh in range(4): transfer[neighbor[k][neigh], k] += 0.25
#initial probability vector: pebble is at position 8 with certainty
position = numpy.zeros(9)
position[8] = 1.0
for t in range(50):
print t,' ',["%0.5f" % i for i in position]
position = numpy.dot(transfer, position) #probability vector at time t
0 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '1.00000'] 1 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.25000', '0.00000', '0.25000', '0.50000'] 2 ['0.00000', '0.00000', '0.06250', '0.00000', '0.12500', '0.18750', '0.06250', '0.18750', '0.37500'] 3 ['0.00000', '0.04688', '0.07812', '0.04688', '0.09375', '0.18750', '0.07812', '0.18750', '0.28125'] 4 ['0.02344', '0.05469', '0.09766', '0.05469', '0.11719', '0.16016', '0.09766', '0.16016', '0.23438'] 5 ['0.03906', '0.07324', '0.10254', '0.07324', '0.10742', '0.15234', '0.10254', '0.15234', '0.19727'] 6 ['0.05615', '0.08057', '0.10767', '0.08057', '0.11279', '0.13989', '0.10767', '0.13989', '0.17480'] 7 ['0.06836', '0.08929', '0.10895', '0.08929', '0.11023', '0.13379', '0.10895', '0.13379', '0.15735'] 8 ['0.07883', '0.09421', '0.11024', '0.09421', '0.11154', '0.12758', '0.11024', '0.12758', '0.14557'] 9 ['0.08652', '0.09871', '0.11057', '0.09871', '0.11089', '0.12373', '0.11057', '0.12373', '0.13657'] 10 ['0.09261', '0.10167', '0.11089', '0.10167', '0.11122', '0.12044', '0.11089', '0.12044', '0.13015'] 11 ['0.09714', '0.10410', '0.11098', '0.10410', '0.11106', '0.11818', '0.11098', '0.11818', '0.12530'] 12 ['0.10062', '0.10582', '0.11106', '0.10582', '0.11114', '0.11638', '0.11106', '0.11638', '0.12174'] 13 ['0.10322', '0.10716', '0.11108', '0.10716', '0.11110', '0.11508', '0.11108', '0.11508', '0.11906'] 14 ['0.10519', '0.10814', '0.11110', '0.10814', '0.11112', '0.11408', '0.11110', '0.11408', '0.11707'] 15 ['0.10666', '0.10889', '0.11110', '0.10889', '0.11111', '0.11334', '0.11110', '0.11334', '0.11557'] 16 ['0.10777', '0.10944', '0.11111', '0.10944', '0.11111', '0.11278', '0.11111', '0.11278', '0.11446'] 17 ['0.10861', '0.10986', '0.11111', '0.10986', '0.11111', '0.11236', '0.11111', '0.11236', '0.11362'] 18 ['0.10923', '0.11017', '0.11111', '0.11017', '0.11111', '0.11205', '0.11111', '0.11205', '0.11299'] 19 ['0.10970', '0.11041', '0.11111', '0.11041', '0.11111', '0.11182', '0.11111', '0.11182', '0.11252'] 20 ['0.11005', '0.11058', '0.11111', '0.11058', '0.11111', '0.11164', '0.11111', '0.11164', '0.11217'] 21 ['0.11032', '0.11071', '0.11111', '0.11071', '0.11111', '0.11151', '0.11111', '0.11151', '0.11190'] 22 ['0.11052', '0.11081', '0.11111', '0.11081', '0.11111', '0.11141', '0.11111', '0.11141', '0.11171'] 23 ['0.11067', '0.11089', '0.11111', '0.11089', '0.11111', '0.11133', '0.11111', '0.11133', '0.11156'] 24 ['0.11078', '0.11094', '0.11111', '0.11094', '0.11111', '0.11128', '0.11111', '0.11128', '0.11145'] 25 ['0.11086', '0.11099', '0.11111', '0.11099', '0.11111', '0.11124', '0.11111', '0.11124', '0.11136'] 26 ['0.11092', '0.11102', '0.11111', '0.11102', '0.11111', '0.11121', '0.11111', '0.11121', '0.11130'] 27 ['0.11097', '0.11104', '0.11111', '0.11104', '0.11111', '0.11118', '0.11111', '0.11118', '0.11125'] 28 ['0.11101', '0.11106', '0.11111', '0.11106', '0.11111', '0.11116', '0.11111', '0.11116', '0.11122'] 29 ['0.11103', '0.11107', '0.11111', '0.11107', '0.11111', '0.11115', '0.11111', '0.11115', '0.11119'] 30 ['0.11105', '0.11108', '0.11111', '0.11108', '0.11111', '0.11114', '0.11111', '0.11114', '0.11117'] 31 ['0.11107', '0.11109', '0.11111', '0.11109', '0.11111', '0.11113', '0.11111', '0.11113', '0.11116'] 32 ['0.11108', '0.11109', '0.11111', '0.11109', '0.11111', '0.11113', '0.11111', '0.11113', '0.11114'] 33 ['0.11109', '0.11110', '0.11111', '0.11110', '0.11111', '0.11112', '0.11111', '0.11112', '0.11114'] 34 ['0.11109', '0.11110', '0.11111', '0.11110', '0.11111', '0.11112', '0.11111', '0.11112', '0.11113'] 35 ['0.11110', '0.11110', '0.11111', '0.11110', '0.11111', '0.11112', '0.11111', '0.11112', '0.11113'] 36 ['0.11110', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112', '0.11111', '0.11112', '0.11112'] 37 ['0.11110', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112', '0.11111', '0.11112', '0.11112'] 38 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112'] 39 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11112'] 40 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 41 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 42 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 43 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 44 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 45 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 46 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 47 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 48 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111'] 49 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111']
Find all of the eigenvalues of the transfer matrix. What do the eigenvalues mean, apart from 1, which is associated with equilibrium?
np.set_printoptions(suppress=True) #suppress scientific format
eigenvalues, eigenvectors = numpy.linalg.eig(transfer) #eigenvalues of the transfer matrix
idx = np.argsort(-eigenvalues) #sort the eigenvalues in descending order
eigenvalues = eigenvalues[idx]
print eigenvalues #print the sorted eigenvalues
# you may print the eigenvectors by uncommenting the following lines...
#for iter in range(9):
# print eigenvalues[iter]
# for i in range(9):
# print eigenvectors[i][iter]
[ 1. 0.75 0.75 0.5 0.25 0.25 -0. -0. -0.5 ]
To obtain the deviation from the equilibrim value at each time, subtract the equilibrium probabilities from each probability vector and take the absolute value. Specifically, the evolution of the deviation in the site 0 is plotted in semi-log scale. It is found that the resulting line is approximately parallel to $0.75^t$ in semi-log scale.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy
neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],
[4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],
[7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]
transfer = numpy.zeros((9, 9))
site_0=numpy.zeros(50)
fit=numpy.zeros(50)
for k in range(9):
for neigh in range(4): transfer[neighbor[k][neigh], k] += 0.25
position = numpy.zeros(9)
position[8] = 1.0
for t in range(50):
print t,' ',["%0.5f" % abs(i- 1.0 / 9.0) for i in position]
position = numpy.dot(transfer, position)
fit[t]=0.75**t
site_0[t]=abs(position[0]- 1.0 / 9.0)
plt.clf()
plt.semilogy(site_0, label="$\pi_0(t)-\pi_0^{eq}$")
plt.semilogy(fit, label="$0.75^t$")
plt.ylabel('Deviation from 1/9')
plt.xlabel('t')
plt.grid(True)
plt.legend(shadow=False, fancybox=True)
plt.show()
0 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.88889'] 1 ['0.11111', '0.11111', '0.11111', '0.11111', '0.11111', '0.13889', '0.11111', '0.13889', '0.38889'] 2 ['0.11111', '0.11111', '0.04861', '0.11111', '0.01389', '0.07639', '0.04861', '0.07639', '0.26389'] 3 ['0.11111', '0.06424', '0.03299', '0.06424', '0.01736', '0.07639', '0.03299', '0.07639', '0.17014'] 4 ['0.08767', '0.05642', '0.01345', '0.05642', '0.00608', '0.04905', '0.01345', '0.04905', '0.12326'] 5 ['0.07205', '0.03787', '0.00857', '0.03787', '0.00369', '0.04123', '0.00857', '0.04123', '0.08615'] 6 ['0.05496', '0.03054', '0.00345', '0.03054', '0.00168', '0.02878', '0.00345', '0.02878', '0.06369'] 7 ['0.04275', '0.02182', '0.00216', '0.02182', '0.00088', '0.02268', '0.00216', '0.02268', '0.04624'] 8 ['0.03228', '0.01690', '0.00087', '0.01690', '0.00043', '0.01647', '0.00087', '0.01647', '0.03446'] 9 ['0.02459', '0.01241', '0.00054', '0.01241', '0.00022', '0.01262', '0.00054', '0.01262', '0.02546'] 10 ['0.01850', '0.00944', '0.00022', '0.00944', '0.00011', '0.00933', '0.00022', '0.00933', '0.01904'] 11 ['0.01397', '0.00701', '0.00014', '0.00701', '0.00005', '0.00707', '0.00014', '0.00707', '0.01419'] 12 ['0.01049', '0.00529', '0.00005', '0.00529', '0.00003', '0.00527', '0.00005', '0.00527', '0.01063'] 13 ['0.00789', '0.00395', '0.00003', '0.00395', '0.00001', '0.00397', '0.00003', '0.00397', '0.00795'] 14 ['0.00592', '0.00297', '0.00001', '0.00297', '0.00001', '0.00297', '0.00001', '0.00297', '0.00596'] 15 ['0.00445', '0.00223', '0.00001', '0.00223', '0.00000', '0.00223', '0.00001', '0.00223', '0.00446'] 16 ['0.00334', '0.00167', '0.00000', '0.00167', '0.00000', '0.00167', '0.00000', '0.00167', '0.00335'] 17 ['0.00250', '0.00125', '0.00000', '0.00125', '0.00000', '0.00125', '0.00000', '0.00125', '0.00251'] 18 ['0.00188', '0.00094', '0.00000', '0.00094', '0.00000', '0.00094', '0.00000', '0.00094', '0.00188'] 19 ['0.00141', '0.00070', '0.00000', '0.00070', '0.00000', '0.00070', '0.00000', '0.00070', '0.00141'] 20 ['0.00106', '0.00053', '0.00000', '0.00053', '0.00000', '0.00053', '0.00000', '0.00053', '0.00106'] 21 ['0.00079', '0.00040', '0.00000', '0.00040', '0.00000', '0.00040', '0.00000', '0.00040', '0.00079'] 22 ['0.00059', '0.00030', '0.00000', '0.00030', '0.00000', '0.00030', '0.00000', '0.00030', '0.00059'] 23 ['0.00045', '0.00022', '0.00000', '0.00022', '0.00000', '0.00022', '0.00000', '0.00022', '0.00045'] 24 ['0.00033', '0.00017', '0.00000', '0.00017', '0.00000', '0.00017', '0.00000', '0.00017', '0.00033'] 25 ['0.00025', '0.00013', '0.00000', '0.00013', '0.00000', '0.00013', '0.00000', '0.00013', '0.00025'] 26 ['0.00019', '0.00009', '0.00000', '0.00009', '0.00000', '0.00009', '0.00000', '0.00009', '0.00019'] 27 ['0.00014', '0.00007', '0.00000', '0.00007', '0.00000', '0.00007', '0.00000', '0.00007', '0.00014'] 28 ['0.00011', '0.00005', '0.00000', '0.00005', '0.00000', '0.00005', '0.00000', '0.00005', '0.00011'] 29 ['0.00008', '0.00004', '0.00000', '0.00004', '0.00000', '0.00004', '0.00000', '0.00004', '0.00008'] 30 ['0.00006', '0.00003', '0.00000', '0.00003', '0.00000', '0.00003', '0.00000', '0.00003', '0.00006'] 31 ['0.00004', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00004'] 32 ['0.00003', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00000', '0.00002', '0.00003'] 33 ['0.00003', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00003'] 34 ['0.00002', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00002'] 35 ['0.00001', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00001'] 36 ['0.00001', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00000', '0.00001', '0.00001'] 37 ['0.00001', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00001'] 38 ['0.00001', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00001'] 39 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 40 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 41 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 42 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 43 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 44 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 45 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 46 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 47 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 48 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000'] 49 ['0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000', '0.00000']
Dual pebble game with a finite probability of switching between the games to ensure irreducability of the transfer matrix.
%pylab qt
import random, pylab
random.seed('1234')
sigma = 0.4
epsilon = 0.4 # probability to switch from red to blue pebble, and vice versa
pylab.figure()
s_map_red = [(1.0, 1.0), (2.0, 1.0), (3.0, 1.0),
(1.0, 2.0), (2.0, 2.0), (3.0, 2.0),
(1.0, 3.0), (2.0, 3.0), (3.0, 3.0)]
offset = 3.0
s_map_blue = [(x+offset,y-offset) for (x,y) in s_map_red]
neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],
[4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],
[7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]
color = 'red' #chose 'red' or 'blue'
site = 8
tmax = 240
for iter in range(tmax):
period = 4
if (iter%period) == 0:
# Begin of graphical output
pylab.clf()
maxlength = len(str(tmax-1))
number_string = str(iter).zfill(maxlength)
if color == 'red': cir = pylab.Circle(s_map_red[site], radius=sigma, fc='r')
if color == 'blue': cir = pylab.Circle(s_map_blue[site], radius=sigma, fc='b')
pylab.figure()
pylab.gca().add_patch(cir)
pylab.plot([0.5, 3.5], [0.5, 0.5], 'r')
pylab.plot([0.5, 3.5], [1.5, 1.5], 'r')
pylab.plot([0.5, 3.5], [2.5, 2.5], 'r')
pylab.plot([1.5, 1.5], [0.5, 3.5], 'r')
pylab.plot([2.5, 2.5], [0.5, 3.5], 'r')
pylab.plot([3.5, 3.5], [0.5, 3.5], 'r')
pylab.plot([0.5+offset, 3.5+offset], [1.5-offset, 1.5-offset], 'b')
pylab.plot([0.5+offset, 3.5+offset], [2.5-offset, 2.5-offset], 'b')
pylab.plot([0.5+offset, 3.5+offset], [3.5-offset, 3.5-offset], 'b')
pylab.plot([0.5+offset, 0.5+offset], [0.5-offset, 3.5-offset], 'b')
pylab.plot([1.5+offset, 1.5+offset], [0.5-offset, 3.5-offset], 'b')
pylab.plot([2.5+offset, 2.5+offset], [0.5-offset, 3.5-offset], 'b')
pylab.title('t = '+ number_string)
pylab.axis('scaled')
pylab.axis([0.5, 6.5, -2.5, 3.5])
pylab.xticks([])
pylab.yticks([])
number_string_filename = str(iter/period).zfill(3)
pylab.pause(0.05)
pylab.close()
# End of graphical output
newsite = neighbor[site][ random.randint(0, 3)]
newcolor = color
if (color == 'red') and (site == 2) and (newsite == 2):
if random.random() < epsilon:
newcolor = 'blue'
newsite = 6
print "transition red->blue at time = ", iter
if (color == 'blue') and (site == 6) and (newsite == 6):
if random.random() < epsilon:
newcolor = 'red'
newsite = 2
print "transition blue->red at time = ", iter
site = newsite
color = newcolor
Populating the interactive namespace from numpy and matplotlib transition red->blue at time = 176 transition blue->red at time = 197 transition red->blue at time = 198 transition blue->red at time = 199 transition red->blue at time = 211 transition blue->red at time = 213
Recurrent pebble game with $2\times 2$ grid with small aperiodicity.
%pylab qt
import math, random, pylab
sigma = 0.4
epsilon = 0.1
pylab.figure()
s_map = [(1.0, 1.0), (2.0, 1.0)]
neighbor = [[1], [0]]
pos = 0
tmax = 20
for iter in range(tmax):
# Begin of the graphics output
pylab.figure()
number_string = str(iter).zfill(len(str(tmax)))
cir = pylab.Circle(s_map[pos], radius=sigma, fc='r')
pylab.gca().add_patch(cir)
pylab.plot([1.5, 1.5], [0.5, 1.5], 'b')
pylab.title('t = '+ number_string)
pylab.axis('scaled')
pylab.axis([0.5, 2.5, 0.5, 1.5])
pylab.xticks([])
pylab.yticks([])
#pylab.savefig('2x1pebble_epsilon'+number_string+'.png', transparent=True)
pylab.pause(0.05)
pylab.clf()
pylab.close()
# End of the graphics output
newpos = neighbor[pos][0]
if random.random() < epsilon:
newpos = pos
pos = newpos
Populating the interactive namespace from numpy and matplotlib
Basic pebble throw game with inhomogeneous steady state probabilities.
import random
#define the grid in terms of transitions:[right,up,left,down]
histo = [0, 0, 0, 0, 0, 0, 0, 0, 0]
neighbor = [[1, 3, 0, 0], [2, 4, 0, 1], [2, 5, 1, 2],
[4, 6, 3, 0], [5, 7, 3, 1], [5, 8, 4, 2],
[7, 6, 6, 3], [8, 7, 6, 4], [8, 8, 7, 5]]
weight = [3.0, 0.5, 1.0, 0.5, 1.0, 0.5, 2.0, 0.5, 1.0]
pos = 8
n_iter = 1000000
for iter in range(n_iter):
new_pos = neighbor[pos][random.randint(0, 3)]
if random.random() < weight[new_pos] / weight[pos]:
pos = new_pos
histo[pos] += 1
norm = sum(weight)
print 'comparison: weight, histogram'
for k in range(9):
print 'site: ', k,' weight: ', weight[k], ' histo: ', norm * histo[k] / float(n_iter)
comparison: weight, histogram site: 0 weight: 3.0 histo: 3.00848 site: 1 weight: 0.5 histo: 0.49825 site: 2 weight: 1.0 histo: 0.99388 site: 3 weight: 0.5 histo: 0.50321 site: 4 weight: 1.0 histo: 1.00448 site: 5 weight: 0.5 histo: 0.49447 site: 6 weight: 2.0 histo: 2.01049 site: 7 weight: 0.5 histo: 0.50085 site: 8 weight: 1.0 histo: 0.98589