I am interesting in creating a mixture of gaussians with defined properties for mean, stdev, skewness, and kurtosis. My thought was to figure out the expressions for each of these descriptive statistics and then solve for the values of $\mu_1$, $\mu_2$, $\sigma_1$, $\sigma_2$, and $p$ ($p$ = proportion of density in first gaussian). This was my approach:
First define the moments for the gaussian mixture: $$ M1 = p \mu_1 + (1 - p) \mu_2\\ M2 = p (\mu_1^2 + \sigma_1^2) + (1 - p) (\mu_2^2 + \sigma_2^2)\\ M3 = p (\mu_1^3 + 3 \mu_1 \sigma_1^2) + (1 - p) (\mu_2^3 + 3 \mu_2 \sigma_2^2)\\ M4 = p (\mu_1^4 + 6 \mu_1^2 \sigma_1^2 + 3 \sigma_1^4) + (1 - p) (\mu_2^4 + 6 \mu_2^2 \sigma_2^2 + 3 \sigma_2^4) $$
Then define the central moments: $$ MM4 = M4 - 4 M3 M1 + 6 M2 M1^2 + 3 M1^4\\ MM3 = M3 - 3 M1 M2 + 2 M1^3\\ MM2 = M2 - M1^2\\ MM1 = M1\\ $$
(NOTE: these equations are simply derived using the equation $MMn=E[(X-E[X])^n]$ and then simplifying)
Then I define the statistics I'm interested in and set them to the values I want, for example:: $$ kurtosis = \frac{MM4 - 3 MM2^2}{MM2^2} = 5\\ skewness = \frac{MM3}{MM2^{\frac{3}{2}}} = 3\\ variance = MM2 = 2\\ mean = MM1 = 0\\ $$
and then solve for $\mu_1$, $\mu_2$, $\sigma_1$, and $\sigma_2$ (4 equations and 4 unknowns).
However, when I plug this into mathematica, it struggles to give me anything (ran for about an hour and then gave up).
I'm wondering if there is a better approach to solve this problem?