I have read around a bit and never asked a question on stackoverflow before. It seems there is a solution to my problem, but I don't understand the way it is written and computed.

Fit mixture of two gaussian/normal distributions to a histogram from one set of data, python

This seems to be the solution to my problem. Let me pose the problem.

At first I thought everything works well, but not quite well

I have a bimodal distribution that is probably somewhat gaussian in its parts, not very much off, and it would be easy, I presumed for a gaussian mixture model to finde a solution.

import numpy as np
from sklearn.mixture import GaussianMixture
from scipy.stats import norm

X = np.asanyarray(histodata)
X = X.reshape(-1, 1)
gm = GaussianMixture(covariance_type="diag",n_components=2, reg_covar=1e-12, tol= 0.1, init_params= "kmeans",random_state=0, max_iter=100).fit(X)
print(gm.means_)
print(gm.covariances_)
print(gm.weights_)

gm.get_params()

So this one is okay. But when I try for another one with two components, it does not work

# find useful parameters
mean = gm.fit(X).means_
covs  = gm.fit(X).covariances_
weights = gm.fit(X).weights_

# create necessary things to plot
x_axis = np.arange(-100, 470, 5)
y_axis0 = norm.pdf(x_axis, float(mean[0][0]), np.sqrt(float(covs[0])))*weights[0] # 1st gaussian
y_axis1 = norm.pdf(x_axis, float(mean[1][0]), np.sqrt(float(covs[1])))*weights[1] # 2nd gaussian

plt.plot(x_axis, y_axis0, lw=3, c='C0')
plt.plot(x_axis, y_axis1, lw=3, c='C1')
plt.plot(x_axis, y_axis0+y_axis1, lw=3, c='C2', ls='dashed')

plt.hist(histodata, density=True, stacked = True, bins = 470, color="black",)

But the peaks are not high enough, the widths not wide enough.