Inverse CDF technique : analytic vs numerical

I am trying to draw samples from a distribution p(M) ~ M^(-1.5) in the interval 10^5 to 10^9. I am comparing my numerical approach to the analytic one. Here is what I have :

import numpy as np
import mpmath as mp
from scipy.integrate import quad
from scipy.optimize import fsolve
from scipy.optimize import brentq
import matplotlib.pyplot as plt

norm= 1/quad(lambda x : (x)**(-1.5),10**5,10**9)[0]
def p(M):
      return norm*(M**(-1.5))

def cdf(M):
      return quad(p,10**5,M)[0]

#analytic
sample_GM_th= lambda: (10.**((-0.5)*5)+(10.**((-0.5)*9)-10.**((-0.5)*5))*np.random.uniform())**(1./(-0.5))


#numerical 
def icdf(M):
     return quad(p,10**(5),M)[0] - np.random.uniform()

sample_GM_num= lambda: brentq(icdf,10**5,10**9)

Now, if I drew say 1000 samples from both of them :

   analytic=[]
   numerical=[]

   for i in range(1000):
         analytic.append(sample_GM_th())
         numerical.append(sample_GM_num())

I expected that the resulting distributions will coincide, but I found that the one from numerical approach is quite different from the analytic one.

Can anybody explain why this is happening?