This is a question related to the statistical model behind independent component analysis (ICA). 


We assume that $Z \sim N(0,1)$. Our goal is to construct a random variable $X$ that follows a Gaussian mixture distribution and satisfies 
$$
E[X^j] = E[Z^j]
$$ 
for $j \leq k-1$, while 
$$
E[X^k] - E[Z^k] = \epsilon.
$$ 
Here $k \geq 2$ is an integer.


My question is: Is there a principled way to construct the Gaussian mixture random variable $X$?


For $k = 4$, one possible construction is $X\sim 1/2 \cdot N(0,1+\varepsilon^{1/2}/\sqrt{3}) + 1/2 \cdot N(0,1-\varepsilon^{1/2}/\sqrt{3})$. However, is there a construction for general $k$?



Moment matching: construction of a mixture of Gaussian distribution with lower moments identical to Gaussian
If non-zero means are allowed, the answer is yes (and it is based on quite different considerations, and so, is presented separately from the previous answer). 
Indeed, let $Z$ be a standard normal random variable (r.v.). 
Let $P\_k$ denote the set of all probability measures $\nu$ on $\mathbb R$ such that $\int\_{\mathbb R}x^j\nu(dx)=EZ^j$ for all $j=0,\dots,k$. Let $M\_k:=\sup\{\int\_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in P\_k\}$ and 
$m\_k:=\inf\{\int\_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in P\_k\}$. Later it will be shown that 
$$m\_k<EZ^{k+1}<M\_k,\tag{-1}
$$ 
with an emphasis on the strictness of the inequalities. 
Take any 
$$c\in(m\_k,M\_k). \tag{0} $$ 
Let us then show that there is a mixture $\mu\_c$ of normal distributions such that $\int\_{\mathbb R}x^{k+1}\mu\_c(dx)=c$ and $\mu\_c\in P\_k$, that is, $\int\_{\mathbb R}x^j\mu\_c(dx)=EZ^j$ for all $j=0,\dots,k$. 


For $t\in(0,1)$, let 
$$c\_t:=c-(1- t^{(k+1)/2})EZ^{k+1}. 
$$
Since $c\_t\to c$ as $t\uparrow1$, one can choose $t$ close enough to $1$ so that, by $(0)$, 
$$c\_t\in(t^{(k+1)/2}m\_k,t^{(k+1)/2}M\_k).$$
So, by the convexity of the set $P\_k$ and my answer to the question at 
**[[mathoverflow.net/questions/229646](https://mathoverflow.net/questions/229646/gaussian-and-the-convex-hull-of-moment-curves)]**
there is a r.v. $X\_t$ taking only finitely many values such that 
$$EX\_t^j=E(\sqrt tZ)^j\text{ for }j=0,\dots,k\text{ and }EX\_t^{k+1}=c\_t. \tag{1}
$$ Without loss of generality, $X\_t$ is independent of $Z$. Let $Y\_t:=X\_t+\sqrt{1-t}Z$. Then the distribution of $Y\_t$ is a mixture of normal distributions (each of those distributions with variance $1-t$). 
Moreover, for $j=0,1,\dots$ 
$$EY\_t^j=\sum\_{i=0}^j\binom ji EX\_t^i\, E(\sqrt{1-t}Z)^{j-i}. \tag{2} 
$$
On the other hand, $Z$ is equal in distribution to $\sqrt t Z\_1+\sqrt{1-t}Z\_2$, where $Z\_1$ and $Z\_2$ are independent standard normal r.v.'s. So, 
$EZ^j=\sum\_{i=0}^j\binom ji E(\sqrt t Z)^i\, E(\sqrt{1-t}Z)^{j-i}$ 
for all $j=0,1,\dots$. 
Comparing this with $(2)$ and recalling $(1)$, we see that $EY\_t^j=EZ^j$ for $j=0,\dots,k$, whereas 
$$EY\_t^{k+1}=EX\_t^{k+1}+\sum\_{i=0}^k\binom{k+1}i EX\_t^i\, E(\sqrt{1-t}Z)^{k+1-i}$$
$$=EX\_t^{k+1}+\sum\_{i=0}^k\binom{k+1}i E(\sqrt t Z)^i\, E(\sqrt{1-t}Z)^{k+1-i}$$
$$
=c\_t+EZ^{k+1}-E(\sqrt t Z)^{k+1}=c, 
$$
as desired.


**Addendum.** 
As promised, I am adding the proof of inequalities $(-1)$. 


**Lemma.** *Let $\mu$ be any finite measure on $\mathbb R$ with finite first $k+1$ moments and with support set $S\_\mu$ of cardinality at least $k+2$, where $k$ is a natural number. For any interval $J\subseteq\mathbb R$, let 
$N\_k(\mu,J)$ denote the set of all measures $\nu$ on $\mathbb R$ such that $\int\_J x^j\nu(dx)=\int\_J x^j\mu(dx)$ for all $j=0,\dots,k$, with $N\_k(\mu):=N\_k(\mu,\mathbb R)$. Let $M\_k(\mu):=\sup\{\int\_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in N\_k(\mu)\}$ and 
$m\_k(\mu):=\inf\{\int\_{\mathbb R}x^{k+1}\nu(dx)\colon\nu\in N\_k(\mu)\}$. Then 
$$m\_k(\mu)<\int\_{\mathbb R}x^{k+1}\mu(dx)<M\_k(\mu),\tag{-1+}
$$ 
with an emphasis on the strictness of the inequalities.*


*Proof.* To obtain a contradiction, suppose that $\int\_{\mathbb R}x^{k+1}\mu(dx)=M\_k(\mu)$. 
Since $S\_\mu$ is of cardinality at least $k+2$, there are pairwise disjoint closed intervals $J\_0,\dots,J\_{k+1}$ of nonzero $\mu$-measure. Again by 
my answer to the question at 
**[[mathoverflow.net/questions/229646](https://mathoverflow.net/questions/229646/gaussian-and-the-convex-hull-of-moment-curves)]**, 
for each $i$ there is a measure $\nu\_i\in N\_{k+1}(\mu,J\_i)$ with a finite support set contained in $J\_i$, so that $p\_i:=\nu\_i(\{x\_i\})>0$ for some $x\_i\in J\_i$. It also follows that the points $x\_i$ are pairwise distinct. 
Let now $\nu(A):=\mu(A\setminus J\_0\setminus\dots\setminus J\_{k+1})+\sum\_{i=0}^{k+1}\nu\_i(A)$ for all Borel sets $A\subseteq\mathbb R$. 
Then $\nu\in N\_{k+1}(\mu)$ and, in particular, $\int\_{\mathbb R}x^{k+1}\nu(dx)=M\_k(\mu)$. 


For the pairwise distinct points $x\_i$ mentioned above, consider the system of equations 
$$\sum\_{i=0}^{k+1}a\_i x\_i^j=I\{j=k+1\}\quad\text{for } j=0,\dots,k+1,
$$
for the unknowns $a\_0,\dots,a\_{k+1}$, where $I\{\cdot\}$ is the indicator function. The determinant of this linear system is a nonzero Vandermonde determinant, and so, the system has a solution $(a\_0,\dots,a\_{k+1})$. For that solution and real $t>0$, let $\nu\_t:=\nu+t\sum\_{i=0}^{k+1}a\_i\delta\_{x\_i}$, where $\delta\_x$ is the Dirac measure at $x$. 
If $t$ is small enough (so that $p\_i+ta\_i\ge0$ for all $i$), then $\nu\_t\in N\_k(\mu)$ and $\int\_{\mathbb R}x^{k+1}\nu\_t(dx)=\int\_{\mathbb R}x^{k+1}\nu(dx)+t=M\_k(\mu)+t>M\_k(\mu)$, which is a contradiction. 
Similarly, the assumption $\int\_{\mathbb R}x^{k+1}\mu(dx)=m\_k(\mu)$ leads to a contradiction. QED


