How to take random points from a moving Gaussian distribution (OU process)

I want to take some points from an OU process, i.e a Gaussian distribution that is moving until it converges to its stationary distribution (here the pdf from wikipedia).

def P(V,D,g,t):
    delta_vt_sqrd=(D/g)*(1-np.exp(-2*g*t))
    avg_vt=1*np.exp(-g*t)
    return (1/np.sqrt(2*np.pi*delta_vt_sqrd))*np.exp(-(V-avg_vt)**2/(2*delta_vt_sqrd))

and I want to take samples from this distribution. But using something like np.random.choice on samples from this and then calculating the mean and variance, doesn't seem to reproduce the original expected variance and mean.

This is also true if I try it with a simple Gaussian distribution, for example:

def P0(d,b): ##OU stationayry
   x=np.arange(-3.2,3.2,0.01)
   a=np.sqrt(b/(2*d*np.pi))
   return a*np.exp(-b*x**2/(2*d))

noise=[]
for i in range (0,1000):
   ou=np.random.choice(P0(1,1))
   noise.append(ou)

won't give me back variance 1 and mean 0 of this sample when I use np.mean(noise),np.var(noise)

So my question is, what am I doing wrong, or where does my misunderstanding lie?