7
$\begingroup$

Given solely the first $n$ moments $m_1,\dots,m_n$ of a random variables $X\in\mathbb{R}$, I was wondering whether there exists a direct methodology to approximate $X$ with a Gaussian Mixture ?

$\endgroup$
  • 1
    $\begingroup$ Are you (implicitly) assuming some fixed number of components $k \leq (n+1)/3$? Or, you wish to try to determine $k$ as well? $\endgroup$ – cardinal Mar 23 '15 at 13:26
  • $\begingroup$ Note that for a simple mixture of $k=2$ Gaussian, you have 5 parameters to estimate, so you would use the first five moments. The equations are a bit involved and you would need a iterative algorithm to find the solutions... EM seems more simple! $\endgroup$ – Elvis Mar 23 '15 at 14:04
  • $\begingroup$ cardinal, yes I also would like to determine the number of components $k$. $\endgroup$ – user149575 Mar 23 '15 at 17:37
  • $\begingroup$ Elvis, I would like to know whether a methodology exists if one aims to match more than the first 5 moments $\endgroup$ – user149575 Mar 23 '15 at 17:38
  • $\begingroup$ For fixed $k$, you need to use at least $n = 3k-1$ moments. If you use exactly this number of moments, you should find a unique solution. It is the usual way to use the method of moments. If you wish to use more moments, you need to define some criterion, like "minimize some distance between the observed and predicted moments" — I am not aware of any use of such method. $\endgroup$ – Elvis Mar 24 '15 at 9:26
3
$\begingroup$

The method of moments can always be used; I assume its properties for Gaussian mixture have been studied but I don’t know any references.

Let’s have a look on the mixture of two Gaussian $\mathcal N(\mu_1, \sigma_1^2)$ and $\mathcal N(\mu_2, \sigma_2^2)$ in proportions $p$, $1-p$. We have five parameters to estimate so we will use the first five moments.

The moment generating function of this mixture is $$p \exp\left(\mu_1 t + {1\over 2} \sigma_1^2 t^2\right) + (1-p) \exp\left(\mu_2 t + {1\over 2} \sigma_2^2 t^2\right)$$ which gives the first five moments as $$\begin{aligned} m_1 &= p \mu_1 + (1-p) \mu_2 \\ m_2 &= p (\mu_1^2 + \sigma_1^2) + (1-p)(\mu_2^2 + \sigma_2^2) \\ m_3 &= p (\mu_1^3 + 3 \mu_1 \sigma_1^2) + (1-p)(\mu_2^3 + 3 \mu_2 \sigma_2^2)\\ m_4 &= p (\mu_1^4 + 6 \mu_1^2 \sigma_1^2 + 3\sigma_1^2) + (1-p)(\mu_2^4 + 6 \mu_2^2 \sigma_2^2 + 3\sigma_2^4)\\ m_5 &= p (\mu_1^5 + 10 \mu_1^3 \sigma_1^2 + 15 \mu_1 \sigma_1^4) + (1-p)(\mu_2^5 + 10 \mu_2^3 \sigma_2^2 + 15 \mu_2 \sigma_2^4) \end{aligned}$$

The difficulty is to solve these equations in $p$, $\mu_1$ and $\mu_2$ for given moments $m_1$, $m_2$ and $m_3$. This is not easy! We need here a iterative algorithm. There is surely something clever to do here but I’ll use brute force, with a quasi-Newton to minimize the norm of the difference:

f <- function(m, p, mu1, s1, mu2, s2) {
  mm1 <- c(mu1, mu1**2 + s1, 3*mu1*s1 + mu1**3, 3*s1**2 + 6*s1*mu1**2 + mu1**4, 15*mu1*s1^2 + 10*s1*mu1^3 + mu1^5)
  mm2 <- c(mu2, mu2**2 + s2, 3*mu2*s2 + mu2**3, 3*s2**2 + 6*s2*mu2**2 + mu2**4, 15*mu2*s2^2 + 10*s2*mu2^3 + mu2^5)
  mm <- p*mm1 + (1-p)*mm2;
  sum( (m-mm)**2 )
}

set.seed(1)
x <- c( rnorm(100, 0, 1), rnorm(200, 4, 0.5) )
m <- c(mean(x), mean(x**2), mean(x**3), mean(x**4), mean(x**5) )

r <- optim(c(0.5,0,1,4,0.25), function(x) f(m, x[1], x[2], x[3], x[4], x[5]), method="BFGS")$par

Let’s see:

hist(x, freq=FALSE, breaks=20)
t <- seq(-3,6,length=501)
lines(t, r[1]*dnorm(t, mean=r[2], sd=sqrt(r[3])) + (1-r[1])*dnorm(t, mean=r[4], sd=sqrt(r[5])), col="red")

enter image description here

This does not look very good. I am pretty sure maximum likelihood behaves better. Moreover it is easy to implement with an EM. I don’t think this deserves more investigations.

$\endgroup$
  • $\begingroup$ Thank you for the interesting answer. It seems that the code you are providing is much faster than using EM, isn't it ? $\endgroup$ – user149575 Mar 24 '15 at 12:37
  • $\begingroup$ If the number of observations is large, yes, it can be much faster. Once the means of the $x_i^k$'s are computed, it does not depends on the sample size anymore, whereas each iteration of the EM is $O(n)$. Was your aim to reduce the time needed to estimate the parameters? $\endgroup$ – Elvis Mar 24 '15 at 13:05
  • 2
    $\begingroup$ I had in mind an application where the moments are available but not the sample. Given that the moments characterize the distribution, it should be that a large number of moments allow to derive a "nice" Gaussian Mixture approximation. There is a large literature on density approximation using polynomials and moments, I was wondering about equivalent techniques with GM. $\endgroup$ – user149575 Mar 24 '15 at 13:20
  • 1
    $\begingroup$ The question then is quite different. I might give it an other shot. $\endgroup$ – Elvis Mar 24 '15 at 13:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.