import numpy as np import scipy as scp from matplotlib import pyplot as plt from RubberBand import RubberBand NUM_BANDS = int(1E5) N = 100 length_dist = np.zeros((NUM_BANDS,)) for i in range(NUM_BANDS): band = RubberBand(N, a=1) length_dist[i] = band.length/band.a P_true = lambda N, L: scp.special.binom(N, (L+N)/2) / 2**N L_norm = np.linspace(-N, N, 100) # normalised (no Link length a) range plt.plot(L_norm, P_true(N, L_norm)) plt.hist(length_dist, range=(-N, N), bins=N, density=True) plt.xlabel("$L/a$") plt.ylabel("P($L/a$)") plt.show()