Is there an equivalent to Fourier decomposition using normal distributions instead of sin / cos? I need to create a desired light pattern (car front light) from a large number of normal distributions (small LEDs, maybe different size). One approach would be to use an evolutionary algorithm, but I was wondering if there is a smarter mathematical approach to this. 
Is there a way to decompose a desired pattern, say for a start a 1d box function, to a list of 1d normal distributions with different sigmas.
This might work similar to Fourier, maybe with the approach to use as less normal distributions as possible, so it will not use a infinite number of very small, narrow normal distributions and place them side by side.
Any ideas how to solve this?
 A: You might want to try out Gaussian Mixture models for your data. 
For example, to decompose a mixture of $\mathcal{N}(10, 5), \mathcal{N}(22, 3)$, using flexmix package
library(flexmix)
set.seed(42)

m1 <- 10
m2 <- 22

sd1 <- 5
sd2 <- 3

N1 <- 1000
N2 <- 5000

D <- c(rnorm(mean = m1, sd = sd1, n = N1), rnorm(mean = m2, sd = sd2, n = N2))

kde <- density(D)

mix1 <- FLXMRglm(family = "gaussian")
mix2 <- FLXMRglm(family = "gaussian")
fit <- flexmix(D ~ 1, data = as.data.frame(D), k = 2, model = list(mix1, mix2))

component1 <- parameters(fit, component=1)[[1]]
component2 <- parameters(fit, component=2)[[1]]

m1.estimated <- component1[1]
sd1.estimated <- component1[2]

m2.estimated <- component2[1]
sd2.estimated <- component2[2]

weights <- table(clusters(fit))


plot(kde)
lines(kde$x, (weights[1]/sum(weights)*dnorm(mean = m1.estimated, sd = 
sd1.estimated, x = kde$x)), col = "red", lwd = 2)
lines(kde$x, (weights[2]/sum(weights)*dnorm(mean = m2.estimated, sd = 
sd2.estimated, x = kde$x)), col = "blue", lwd = 2)


A: There are several important differences between sinusoidal functions are normal distributions. Sinusoidal functions provide an orthogonal basis for "nice" functions (I won't go into defining "nice" rigorously). So not only can every "nice" function be written as a linear combination, but the coefficients can be calculated independently as well, because functions with different frequencies are orthogonal. Normal functions are not orthogonal (they are never negative, so their product is never negative, so it's impossible for the dot product to have negative components, and thus cannot be zero).  Furthermore, if you vary only the sigmas, they can't possibly be a basis for any set that contains functions that are not symmetric about the mean. If you vary the mean and variance, then you have a two dimensional space of functions to search over. This set will span the set of "nice" functions, but will in some sense not be linearly independent, and coefficients can't be obtained by simply taking the dot product of the function with gaussian functions, as they can with sinusoidal functions. Thus, while there are packages that approximate a given function with gaussians, they are more complicated than fourier packages, are non-linear, and generally ask you to specify beforehand how many gaussians you wish to use to approximate the given function. (Note that in rightskewed's answer, they tell the package that they're looking to decompose the function into two gaussians with the parameter k=2).
