import random
from math import exp
from matplotlib import pyplot as plt
"""
Statistical results of random trade:
Start from the uniform situation where everyone has the same money,
repeat the cycle to exchange a constant dollars with a randomly selected person.
"""
#=============================
# 大域変数の定義
#=============================
# Number of persons
npersons = 50
# The initial money that the persons have, equal to the average money
value = 50
# The money traded at each cycle
vtrade = 1
# The number of maximum iteration
niter = 500
# The cycle to update graph
nplotinterval = 100
# The number of x mesh for distribution function
nfx = 21
# Treat argments
argv = sys.argv
if len(argv) == 1:
print("Usage: python randomtrade.py npersons value(average) vtrade n(maxiteration) n(plotinterval) n(distribution func)")
print(" ex: python randomtrade.py 200 50 1 10000 100 21")
exit()
if len(argv) >= 2:
npersons = int(argv[1])
if len(argv) >= 3:
value = int(argv[2])
if len(argv) >= 4:
vtrade = int(argv[3])
if len(argv) >= 5:
niter = int(argv[4])
if len(argv) >= 6:
nplotinterval = int(argv[5])
if len(argv) >= 7:
nfx = int(argv[6])
def main():
# 初期化
xindex = range(npersons) # list of indexes of the persons to be used for x axis
v = [] # list of money the persons have
for i in range(npersons):
v.append(value)
#=============================
# グラフの表示
#=============================
# theoretical distribution: f = A * exp(-bv)
# A and b are determined by the number of persons and the money average
print("")
print("Analyatical distribution:")
xmax = value * 10.0
nftx = 101
xstep = xmax / (nftx - 1)
b = 1.0 / value
A = npersons * b
print(" x range: 0 - %g" % (xmax))
print(" x step: %g" % (xstep))
print(" f(x) = %g * exp(-%g * x)" % (A, b))
ftx = []
fty = []
for i in range(nftx+1):
x = i * xstep
ftx.append(x)
fty.append(A * exp(-b * x))
# Prepare plot
fig = plt.figure()
ax1 = fig.add_subplot(3, 1, 1)
ax2 = fig.add_subplot(3, 1, 2)
ax3 = fig.add_subplot(3, 1, 3)
plt.subplots_adjust(wspace = 0.2, hspace = 0.3)
# trade
for i in range(niter):
for ip in range(npersons):
itarget = int(random.random() * npersons)
if v[ip] >= vtrade:
v[itarget] += vtrade
v[ip] -= vtrade
else:
v[itarget] += v[ip]
v[ip] = 0
vs = sorted(v)
if i == 0 or (i+1) % nplotinterval == 0 or i == niter - 1:
# calculate distribution function f(x)
vmax = max(v)
vstep = vmax / (nfx - 1)
fx = []
fy = []
for ifx in range(nfx):
vlim0 = ifx * vstep
vlim1 = (ifx+1) * vstep
fx.append(vlim0)
fy.append(0.0)
for ip1 in range(npersons):
if vlim0 < vs[ip1] <= vlim1:
fy[ifx] += 1.0 / vstep
ax1.cla()
ax2.cla()
ax3.cla()
ax1.bar(xindex, v, width = 1.5)
ax1.plot([0, npersons], [value, value], linestyle = 'dashed', color = 'red', linewidth = 0.5)
ax1.set_title('Random trade (i = %d/%d, n=%d, avg=%d)' % (i+1, niter, npersons, value))
ax1.set_xlabel("i")
ax1.set_ylabel("value")
ax1.set_xlim([0, npersons])
ax2.bar(xindex, vs, width = 1.5)
ax2.plot([0, npersons], [value, value], linestyle = 'dashed', color = 'red', linewidth = 0.5)
ax2.set_xlim([0, npersons])
ax2.set_xlabel("i")
ax2.set_ylabel("value (sorted)")
ax3.plot(fx, fy, label = 'distrib')
ax3.plot(ftx, fty, label = 'exponential', color = "red", linewidth = 0.5, linestyle = 'dashed')
ax3.set_xlim([0, vmax])
ax3.set_ylim([0, max(fy) + 1])
ax3.set_xlabel("value")
ax3.set_ylabel("frequency / value")
ax3.legend()
plt.pause(0.1)
print("Press ENTER to exit>>", end = '')
input()
if __name__ == '__main__':
main()
