# Method of moments

For mixtures of three Gaussians

Consider three Gaussian random variables $X_1,X_2,X_3$ with means $\mu_1,\mu_2,\mu_3$ and variances $\sigma_1^2,\sigma_2^2,\sigma_3^2$. The density of $X_i$ is

$$\phi_i(x) = \frac{1}{\sqrt{2\pi}}\; e^{-\frac{(x-\mu_i)^2}{2\sigma_i^2}}.$$

A mixture of the three random variables is the random variable $U$ with density

$$\psi(x) = a_1\, \phi_1(x)\, + \,a_2 \,\phi_2(x)\, + \,a_3 \,\phi_3(x), \quad\text{ for } \quad a_1+a_2+a_3 =1.$$

The method of moments recovers $\psi$ from the moments

$$m_k = \int x^k \,\psi(x)\; \mathrm{d}x.$$

Since we have 8 unknowns, we expect to need at least 8 moments to recover $\psi$. Let us set up a system for this in Julia.

@polyvar a[1:3] μ[1:3] σ²[1:3]

m0 = a+a+a;
m1 = a*μ+a*μ+a*μ;
m2 = a*(μ^2+σ²)+a*(μ^2+σ²)+a*(μ^2+σ²);
m3 = a*(μ^3+3*σ²*μ)+a*(μ^3+3*σ²*μ)+a*(μ^3+3*σ²*μ);
m4 = a*(μ^4+6*σ²*μ^2+3*σ²^2)+a*(μ^4+6*σ²*μ^2+3*σ²^2)+
a*(μ^4+6*σ²*μ^2+3*σ²^2);
m5 = a*(μ^5+10*σ²*μ^3+15*μ*σ²^2)+
a*(μ^5+10*σ²*μ^3+15*μ*σ²^2)+
a*(μ^5+10*σ²*μ^3+15*μ*σ²^2);
m6 = a*(μ^6+15*σ²*μ^4+45*μ^2*σ²^2+15*σ²^3)+
a*(μ^6+15*σ²*μ^4+45*μ^2*σ²^2+15*σ²^3)+
a*(μ^6+15*σ²*μ^4+45*μ^2*σ²^2+15*σ²^3);
m7 = a*(μ^7+21*σ²*μ^5+105*μ^3*σ²^2+105*μ*σ²^3)+
a*(μ^7+21*σ²*μ^5+105*μ^3*σ²^2+105*μ*σ²^3)+
a*(μ^7+21*σ²*μ^5+105*μ^3*σ²^2+105*μ*σ²^3);
m8 = a*(μ^8+28*σ²*μ^6+210*μ^4*σ²^2+420*μ^2*σ²^3+105*σ²^4)+
a*(μ^8+28*σ²*μ^6+210*μ^4*σ²^2+420*μ^2*σ²^3+105*σ²^4)+
a*(μ^8+28*σ²*μ^6+210*μ^4*σ²^2+420*μ^2*σ²^3+105*σ²^4)

f = [m0, m1, m2, m3, m4, m5, m6, m7, m8]


Let us consider the following moments:

p₀ =  [1; -1; 3; -5.5; 22.45; -50.75; 243.325; -635.725; 3420.7375]


Solving for $f - p₀ = 0$ directly is difficult. The total number of paths is

julia> bezout_number(f - p₀)
362880


yet, Amendola, Faugere and Sturmfels showed that the number of (complex) solutions of the polynomial system $f = p₀$ is 1350 (much less than 362880!).

Instead of solving $f - p₀$ directly, we use monodromy for solving another instance $f - q₀$. Then we move the 1350 computed solutions from $f - q₀$ to $f - p₀$. First, we generate a start solution:

a₀ = randn(ComplexF64, 3)
a₀ = a₀ / sum(a₀)
μ₀ = rand(ComplexF64, 3)
σ²₀ = rand(ComplexF64, 3)
start_sol = [a₀; μ₀; σ²₀]

q₀ = [m([a; μ; σ²] => start_sol) for m in f]


which we can use as input data for the monodromy method for q₀ being the parameters. Here, it doesn’t matter that complex solutions have no interpretation as parameters. They only serve to solve one specific instance of the polynomial system.

julia> @polyvar p[1:9]
julia> R = monodromy_solve(f - p, start_sol, q₀, parameters=p; target_solutions_count = 1350)
MonodromyResult
==================================
• 1350 solutions (0 real)
• return code → success
• 11009 tracked paths


Now, we can track the solutions from $q₀$ to $p₀$.

R2 = solve(f - p, solutions(R), parameters = p, start_parameters=q₀, target_parameters = p₀)


There are

julia> nreal(R2)
36


real solutions. Let us keep only those for which the variances are positive:

true_real_solutions  = filter(r -> all(r[7:9] .> 0), real_solutions(R2))


There are 12 of those. Those 12 come in groups of 6, because the group that permutes ${1,2,3}$ acts on the solutions. We can filter this as follows

julia> S₃ = SymmetricGroup(3)
julia> relabeling(v) = map(p -> (v[1:3][p]..., v[4:6][p]..., v[7:9][p]...), S₃)
julia> Mults = multiplicities(true_real_solutions, group_action = relabeling)
[3, 4, 8, 9, 11, 12]
[1, 2, 5, 6, 7, 10]


Let us thus look at the first and at the third real solution

julia> true_real_solutions
[-0.182136, -1.74012, 2.92226, -1.60721, 0.184633, -0.33243, 0.0927855, 1.11301, 1.76594]
julia> true_real_solutions
[1.30343, -0.396056, 0.0926208, -0.346632, 0.755685, -2.68724, 1.98223, 0.553889, 0.392889]


Those are the parameters of the two mixtures Gaussian that give our moments q₀.

The group action can also be exploited in the monodromy computation:

julia> R_with_group_action = monodromy_solve(f - p, start_sol, q₀, parameters=p, group_action = relabeling; target_solutions_count = 225)
MonodromyResult
==================================
• 225 classes of solutions (modulo group action) (0 real)
• return code → heuristic_stop
• 8072 tracked paths


The full 1350 solutions can then be returned as

relabeling.(solutions(R_with_group_action))


Then, we can proceed with the parameter homotopy as above.

Cite this example:
@Misc{ method_of_moments2019 ,
author =  { Paul Breiding },
title = { Method of moments },
howpublished = { \url{ https://www.JuliaHomotopyContinuation.org/examples/method_of_moments/ } },
note = { Accessed: November 5, 2019 }
}