AICc is picking overly complex models - something stricter? I'd like to know if there are stricter alternatives to automated model selection than AICc / AIC / BIC.
We have approximately ten thousand curves, and for each we'd like to find the most parsimonious Gaussian mixture model. Our maximum complexity is a mixture of 5 Gaussians. The procedure so far is to fit these 5 models:
gauss1 = a1*exp(-((x-b1)/c1)^2)
gauss2 = a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)
...
gauss5 = a1*exp(-((x-b1)/c1)^2) + ... + a5*exp(-((x-b5)/c5)^2),

and then calculate AICc for each and choose the model with the lowest AICc value. We repeat this for each of the ten thousand curves.
To verify whether it's working, I generated some simulation data which looks a lot like our real data. The simulated data is generated by models gauss1, gauss2, ..., gauss5 + noise. To include noise, I added a random number to each point (between 0 and 10% of the maximum value of the simulated data, uniform distribution). This way I know what the model complexity should be. 
Unfortunately it looks like AICc is choosing overly complex models. An example fit for simulated data is below. The complexity of the generating model is 2 Gaussians (gauss2) and the best fitted model, as chosen by AICc, was gauss5.

Using 50 simulated curves, here is the model complexity breakdown. The simulated curves were chosen to be biased toward a simpler model. However the complexity of the fitted models, as chosen by AICc, is biased toward the more complex. For example from these 50 curves, AICc never chose gauss1, even though 19 out of 50 were generated by gauss1. Similar results are seen for AICc, AIC, and BIC.

Choosing between the 5 models with AICc (or AIC or BIC) does not look like it's working for us. Are there other accepted ways to pick between a small number of models? I don't think I can use LASSO / ridge regression since I believe that's only for linear regression. I also don't want to build something ad hoc.
 A: You could try cross-validation, fitting all five models for each fold and picking the model that has the lowest MSE (or whatever other error measure you are interested in) on the holdout folds. That is, for each curve, label all dots at random with labels "1", "2", ..., "10" (for ten-fold cross validation). Fit a one-Gaussian model to all points not labeled "1", and using this fitted model, predict the $y$ component of the holdout data (i.e., the points labeled "1"). Record the Mean Squared Error. Do the same for two-Gaussian, three-Gaussian etc. models. Then repeat the process for points labeled "2". Finally, summarize the MSEs per mixture model over all labels, and you will get grand MSEs per mixture model. Pick the model with the lowest overall MSE.
That said, maybe it's not all that problematic if you have "too complicated" models? The problem with too complex models is usually that they are too variable, but looking at your example with two Gaussians in the true DGP and five in the AICc-fitted model, it's not like the five-Gaussian fit is obviously problematic. This will depend on what you actually want to do with your data. In the particular example you show, it looks to me like (say) an interpolation based on five Gaussians wouldn't be overly dramatically different from one based on two Gaussians.
Edit: Alternatively, you could do a Bayesian approach. Assign prior probabilities to the number of mixture components, which are higher for one or two than for five Gaussians. Next, assign sufficiently uninformative priors to the actual parameters of these $k$ Gaussians. Derive posterior distributions and pick the number of Gaussians with the highest posterior probability, which would be a kind of a maximum a posteriori (MAP) approach. (Or use the full posterior distribution, which would be a kind of a "mixture of mixtures", since it would mix one-Gaussian, two-Gaussian and so forth mixtures.) This approach has the advantage that you can fine-tune the (non-)complexity you want, by changing the prior probabilities you assign to the numbers of mixture components.
A: K fold cross validation using the one standard deviation rule. With the sd rule, the model with the lowest error is not selected. Rather, the least complex model that is within one sd of the model with the lowest error is select as the best model. Please see the r package bestglm.
