Fitting 4-moment distribution with mixture gaussian I know that Mclust does the fit on its own but I am trying to implement an optimization with the aim to generate a mixture of 2 gaussians with the combine moments as closed as possible to the moment of my returns' distribution. 
The objective is to 
Min Abs((Mean Ret - MeanFit)/Mean Fit) + Abs((Std Ret -Stdev Fit)/Stdev) + Abs((Sk Ret-Sk fit)/Sk Fit) + Abs((Kurt Ret- Kurt Fit)) 
Taking into account that I fix the weight between the two gaussians at (0.2;0.8) I implement the below code in R: 
distance <-function(parameter,x) { 
 u=mean(x) 
 s=sd(x) 
 sk=skewness(x) 
 kurt=kurtosis(x) 
 d1=dnorm(x,parameter[1],parameter[2]) 
 d2=dnorm(x,parameter[3],parameter[4]) 
 dfit=0.2d1+0.8d2 
 ufit=mean(dfit) 
 sdfit=sd(dfit) 
 skfit=skewness(fit) 
 kurtfit=kurtosis(fit) 


 abs((u-ufit)/ufit)+abs(s-sdfit)/sdfit)+abs((sk-skfit)/skfit)+abs((kurt-kurtfit)/kurtfit)) 
} 
Parameter<-c(0,0.01,0,0.01)  # starting point of the optimization 
opp<-optim(parameter,distance,x=conv) 



*

*could anybody tell me whether it is the right approach ? 

*should I add some constraint like 
ufit=0.2*mean(d1)+0.8*mean(d2)...
thank you very much in advance for your time and help. 
Sam
 A: I am not exactly sure what your code is trying to do, but it seems like you should be using $rnorm()$ (with a large $n$) instead of $dnorm()$, since the functions $mean()$, $var()$, etc. are designed to be used on samples, not densities.
In any case, using the method-of-moments is easy to do algebraically.
Let's say your sample is $X=[x_1, \cdots, x_n]$.  Your sample moments are:
$$m_1(X) = (1/n) \sum_{i} x_i$$
$$m_2(X) = (1/n) \sum_{i} x_i^2$$
$$m_3(X) = (1/n) \sum_{i} x_i^3$$
$$m_4(X) = (1/n) \sum_{i} x_i^4$$
The moments of $Z \sim N(\mu, \sigma^2)$ are:
$$m_1(Z) = \mu$$
$$m_2(Z) = \mu^2 + \sigma^2$$
$$m_3(Z) = \mu^3 + 3\mu\sigma^2$$
$$m_4(Z) = \mu^3 + 6\mu^2\sigma^2 + 3(\sigma^2)^2$$
Furthermore, if your random variable $Y$ is a mixture of $Z_1, Z_2, \cdots, Z_m$, with $P[X=Z_i]=p_i$, then
$$m_1(Y) = p_1 m_1(Z_1) + p_2 m_1 (Z_2) + \cdots + p_m m_1 (Z_m)$$
$$m_2(Y) = p_1 m_2(Z_1) + p_2 m_2 (Z_2) + \cdots + p_m m_1 (Z_m)$$
$$\cdots$$
If you use enough gaussians you should be able to match a finite number of moments exactly.  Under constraint you may need to define a distance metric on your moments vector $(m_1, m_2, m_3, m_4)$.  The euclidean distance should work fine but you will want to differentially weight the moments, since the higher-order moments could be very large or very small depending on what distribution you are working with.  Also, the first and second moments are generally more important.
