Fitting a mixture of two normal distributions for a data set?

EDIT: Raw data can be found here: https://gist.github.com/Kagaratsch/65a931d8d78fcdd81f7e346429a02afd

Consider the following binned example data:

hl={{-(153/400), 1}, {-(151/400), 0}, {-(149/400), 0}, {-(147/400), 0}, {-(29/80), 0}, {-(143/400), 0}, {-(141/400), 0}, {-(139/400), 0}, {-(137/400), 0}, {-(27/80), 0}, {-(133/400), 0}, {-(131/400), 0}, {-(129/400), 0}, {-(127/400), 0}, {-(5/16), 0}, {-(123/400), 0}, {-(121/400), 0}, {-(119/400), 0}, {-(117/400), 0}, {-(23/80), 0}, {-(113/400), 1}, {-(111/400), 0}, {-(109/400), 0}, {-(107/400), 0}, {-(21/80), 0}, {-(103/400), 0}, {-(101/400), 0}, {-(99/400), 0}, {-(97/400), 0}, {-(19/80), 0}, {-(93/400), 0}, {-(91/400), 0}, {-(89/400), 0}, {-(87/400), 0}, {-(17/80), 0}, {-(83/400), 3}, {-(81/400), 0}, {-(79/400), 0}, {-(77/400), 1}, {-(3/16), 0}, {-(73/400), 0}, {-(71/400), 1}, {-(69/400), 3}, {-(67/400), 4}, {-(13/80), 4}, {-(63/400), 5}, {-(61/400), 3}, {-(59/400), 2}, {-(57/400), 5}, {-(11/80), 8}, {-(53/400), 4}, {-(51/400), 8}, {-(49/400), 8}, {-(47/400), 11}, {-(9/80), 13}, {-(43/400), 10}, {-(41/400), 11}, {-(39/400), 18}, {-(37/400), 13}, {-(7/80), 21}, {-(33/400), 24}, {-(31/400), 28}, {-(29/400), 18}, {-(27/400), 35}, {-(1/16), 40}, {-(23/400), 39}, {-(21/400), 40}, {-(19/400), 41}, {-(17/400), 45}, {-(3/80), 58}, {-(13/400), 47}, {-(11/400), 59}, {-(9/400), 55}, {-(7/400), 71}, {-(1/80), 85}, {-(3/400), 70}, {-(1/400), 65}, {1/400, 83}, {3/400, 85}, {1/80, 83}, {7/400, 68}, {9/400, 73}, {11/400, 66}, {13/400, 61}, {3/80, 70}, {17/400, 60}, {19/400, 63}, {21/400, 48}, {23/400, 52}, {1/16, 46}, {27/400, 34}, {29/400, 43}, {31/400, 36}, {33/400, 27}, {7/80, 21}, {37/400, 23}, {39/400, 13}, {41/400, 17}, {43/400, 26}, {9/80, 9}, {47/400, 15}, {49/400, 6}, {51/400, 7}, {53/400, 5}, {11/80, 5}, {57/400, 8}, {59/400, 2}, {61/400, 2}, {63/400, 4}, {13/80, 2}, {67/400, 4}, {69/400, 3}, {71/400, 3}, {73/400, 5}, {3/16, 1}, {77/400, 3}, {79/400, 0}, {81/400, 3}, {83/400, 1}, {17/80, 1}, {87/400, 0}, {89/400, 1}, {91/400, 0}, {93/400, 5}, {19/80, 0}, {97/400, 1}, {99/400, 1}, {101/400, 0}, {103/400, 0}, {21/80, 1}, {107/400, 0}, {109/400, 0}, {111/400, 0}, {113/400, 0}, {23/80, 2}, {117/400, 0}, {119/400, 1}, {121/400, 0}, {123/400, 0}, {5/16, 0}, {127/400, 0}, {129/400, 0}, {131/400, 1}, {133/400, 0}, {27/80, 1}, {137/400, 0}, {139/400, 0}, {141/400, 0}, {143/400, 0}, {29/80, 0}, {147/400, 0}, {149/400, 0}, {151/400, 0}, {153/400, 0}, {31/80, 0}, {157/400, 0}, {159/400, 0}, {161/400, 0}, {163/400, 0}, {33/80, 0}, {167/400, 0}, {169/400, 0}, {171/400, 0}, {173/400, 0}, {7/16, 0}, {177/400, 0}, {179/400, 0}, {181/400, 1}, {183/400, 1}, {37/80, 0}, {187/400, 0}, {189/400, 0}, {191/400, 0}, {193/400, 0}, {39/80, 0}, {197/400, 0}, {199/400, 0}, {201/400, 0}, {203/400, 0}, {41/80, 1}};
ListLinePlot[hl]


I would like to fit a sum of two normal distributions into this data, so I try

mod = NonlinearModelFit[hl, A1 Exp[-A2 (x - A3)^2] + B1 Exp[-B2 (x - B3)^2], {A1, A2, A3, B1, B2, B3}, x] // Normal;


Mathematica complains that there are convergence issues, and sure enough a plot of the result is very unsatisfactory:

Show[ListLinePlot[hl, PlotRange -> All], Plot[mod, {x, -0.3, 0.3}, PlotStyle -> Red]]


What is the proper way to do this fit in Mathematica, so that it actually converges to a sensible approximation?

EDIT

Interestingly, comparing the (normalized) naive fit to the mixture and smooth kernel distributions from the answer by JimB we see that the fit deviates from the distributions quite a bit

Show[Plot[PDF[mixture /. sol, z], {z, -0.4, 0.4}],
Plot[mod, {x, -0.4, 0.4}, PlotStyle -> Red],
Plot[SKD, {x, -0.4, 0.4}, PlotStyle -> Green]]


• Do you have frequency counts or does the data consist of pairs of measurements? If the former, then NonlinearModelFit is inappropriate. If the latter note that the model (a mixture of two curves with a similar shape as a normal distribution) assumes equal variability across all values which the data does not exhibit. There's much less variability in the tails than in the middle. – JimB Jun 22 '19 at 16:14
• @JimB Those are frequency counts. Right, so my question is - how to fit a sum of two Gaussians into a distorted bell curve? I don't have any strong attachment to the NonlinearModelFit function. Please, let me know if there is a better function for the job? – Kagaratsch Jun 22 '19 at 16:19
• @Kagaratsch: Writing respectively A2^2 and B2^2 you will get what you want. – TeM Jun 22 '19 at 16:22
• @TeM Amazing, you are right! That is very curious... – Kagaratsch Jun 22 '19 at 16:25
• And to be picky: you have a "mixture" of normal densities (which is a weighted sum of the densities) rather than a "sum" of two normal random variables. You might want to change "sum" in the title to "mixture". – JimB Jun 22 '19 at 18:41

Statistics is more than mathematics. One needs to account for how the data was collected rather than just starting with the data and applying some analysis procedure.

What you have is a random sample from a distribution that you've hypothesized to be a mixture of two normal distributions. (The initial attempt at using regression is a common misconception that seems to be prevalent in this forum. I have to believe that this approach must be (inappropriately) used in subject matter textbooks because it seems to occur so often.)

Using the data you provided it is relatively simple in Mathematica to fit a mixture of normal distributions:

mixture = MixtureDistribution[{w1, 1 - w1},
{NormalDistribution[μ1, σ1], NormalDistribution[μ2, σ2]}]

sol = FindDistributionParameters[data, mixture]
(* {w1 -> 0.964246, μ1 -> 0.00764751, σ1 -> 0.0853816, μ2 -> 0.208146, σ2 -> 0.189363} *)
Plot[PDF[mixture /. sol, z], {z, Min[data], Max[data]}]


Unfortunately FindDistributionParameters does not supply standard errors or covariance among the parameter estimators. But that is not too difficult either.

(* Log of the likelihood *)
logL = LogLikelihood[mixture, data];

(* Parameter covariance matrix *)
cov = -Inverse[(D[logL, {{w1, μ1, σ1, μ2, σ2}, 2}]) /. sol];

(* Standard errors *)
se = Thread[{sew1, seμ1, seσ1, seμ2, seσ2} -> Diagonal[cov]^0.5]
(* {sew1 -> 0.013437142118899128,seμ1 -> 0.0021502023883548864,
seσ1 -> 0.0018001069575776648,seμ2 -> 0.05745078807898059,
seσ2 -> 0.022206958940369257} *)


Plot[{PDF[NormalDistribution[Mean[data], StandardDeviation[data]], z],

• Let's just say that I knew of the hammer called NonlinearModelFit and so the problem looked a lot like a nail to me. :) – Kagaratsch Jun 22 '19 at 18:46
• Good way of putting it. (We all have analogous hammers.) You are not alone concerning NonlinearModelFit`. – JimB Jun 22 '19 at 18:51