import numpy.random as rd

def f(n) :
    L = [1]
    for k in range(n) :
        nb = 2*k+2
        b = L[-1]
        r = rd.randint(1,nb+1)
        if r <= b :
            b += 1
        L.append(b)
    return L
    
def approx(n) :
    E,E2 = 0,0
    N = 100
    for k in range(N) :
        b = f(n)[-1]
        E += b
        E2 += b**2
    return E/N,E2/N
    
N = [k for k in range(31)]
B = [approx(n) for n in N]
plt.plot(N,[B[k][0] for k in range(len(B))])
plt.plot(N,[B[k][1] for k in range(len(B))])
plt.show()