Histogram of a Sample with Overlay of Population Density To familiarize myself with histograms and probability density functions, I decided to sample various distributions, plot samples' histograms and their corresponding probability distribution functions.
I started with Beta:
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta, norm

rng = np.random.default_rng()

# Generate data
a = 2.; b = 6 
s = rng.beta(a,b,10000)

# Plot histogram
fig, ax = plt.subplots()
ax.hist(s, 50, density=True, label=r'Sample counts: $\alpha$=2, $\beta$=6')

# Plot pdf
x = np.linspace(beta.ppf(0.001,a,b), beta.ppf(0.999,a,b), 100)
ax.plot(x, beta.pdf(x, a, b),'-', lw=2, color='red', alpha=0.8, label='Beta proba dist function')

ax.set_xlabel('x', fontsize=12, fontweight='bold')
ax.set_ylabel('Beta pdf', fontsize=12, fontweight='bold')
plt.legend(loc='upper right')

plt.show()


Thankfully the result is as-expected.
The trouble started with the normal distribution (standard or not):
mu = -2.
sigma = 2.

# Generate data
data = rng.normal(mu,sigma,(10000,1))

# Plot histogram
fig, ax = plt.subplots() 
count, bins, ignored = ax.hist(data, bins=100, color = (0.,0.,1,0.6))

# Plot pdf
x = bins
y = np.exp(- (bins - mu)**2 / (2 * sigma**2) ) / ( sigma * np.sqrt(2 * np.pi) )
# or
x = np.linspace(norm.ppf(0.001,loc=mu,scale=sigma), norm.ppf(0.999,loc=mu,scale=sigma), 1000)
y = norm.pdf(x,loc=mu,scale=sigma)
ax.plot(x, y, color=(1,0,0,0.8), lw=2, label='normal proba dist. function')

plt.axvline(data.mean(), 
            color='r', linestyle='dotted', linewidth=2, 
            label='Distribution mean' + ' (' + str(round(data.mean(),1)) + ')')
ax.set_facecolor((0.4,0.4,0.1,0.3))
ax.set_ylim(0,400)

plt.legend(fontsize=10, loc='upper left', bbox_to_anchor=(0, 1), ncol=1)

plt.show()


Obviously I expected to get a red bell-shaped curved located at x=-2 but, no dice, there is what is probably a trivial vertical scaling problem with the generated normal curve's pdf. I am missing something basic and don't know what it is.
Any pointer appreciated.

Thanks for editing your Question so it could be re-opened.
When showing simulated data or illustrating an analysis that involves simulation, it is always
a good idea to set the seed of the pseudorandom
number generator so others can replicate your work.
I believe @rxFt20 (+1) has the right idea; you need to
plot your histogram of normal data on a 'density'
scale if you want to overlay a density function on the histogram. That is, the sum of the areas of the bars must be unity. Note the label 'Density' on the vertical axis.
I will illustrate the same procedure in the base of R, where
a density histogram is invoked by the parameter
prob=T. (More sophisticated graphics in R would use different code, but my intention is to illustrate the
idea of a density histogram, not the R code.)
set.seed(2022)
x = rnorm(10^6, -2, 2)
hist(x, prob=T, br=50, col="skyblue2")
curve(dnorm(x, -2, 2), add=T, 
 col="brown", lwd=2)


In R, with an additional line of code, you can also show a kernel density estimate of the population density function, based on the data (and independent of the histogram binning). With 100,000 observations, there is not much difference between the actual population density and
the KDE (dotted green curve).
...
lines(density(x), col = "green3", 
 lwd=3, lty = "dotted")


