I am trying to find a numerical solution to the following nonlinear system of equations:
M_1  w_1 θ_1 α_1  w_2 θ_2 α_2 = 0
M_2  w_1 θ_1^2 (α_1+1) α_1  w_2 θ_2^2 (α_2+1) α_2 = 0
M_3  w_1 θ_1^3 (α_1+2) (α_1+1) α_1  w_2 θ_2^3 (α_2+2) (α_2+1) α_2 = 0
M_4  w_1 θ_1^4 (α_1+3) (α_1+2) (α_1+1) α_1  w_2 θ_2^4 (α_2+3) (α_2+2) (α_2+1) α_2 = 0
M_5  w_1 θ_1^5 (α_1+4) (α_1+3) (α_1+2) (α_1+1) α_1
~~~~~~  w_2 θ_2^5 (α_2+4) (α_2+3) (α_2+2) (α_2+1) α_2 = 0
M_6  w_1 θ_1^6 (α_1+5) (α_1+4) (α_1+3) (α_1+2) (α_1+1) α_1
~~~~~~  w_2 θ_2^6 (α_2+5) (α_2+4) (α_2+3) (α_2+2) (α_2+1) α_2 = 0
The unknowns are \alpha_1, \alpha_2, \theta_1, \theta_2, w_1 and w_2; the M_i~ (i = 1, ..., 6) are known.
In addition, it is known that all unknowns are greater than zero, and that w_1 + w_2 = 1.
(For context: The M_i are moments of a 2component mixture distribution, (\alpha_1, \theta_1) and (\alpha_2, \theta_2) are the parameters of the two mixture components, and w_1 and w_2 are their respective weights.)
I tried the following approach using NLsolve:
using NLsolve
function f!(F, x, M)
# x: [α1, α2, θ1, θ2, w1, w2]
F[1] = M[1]  x[5]*x[3]*x[1]  x[6]*x[4]*x[2]
F[2] = M[2]  x[5]*x[3]^2*(x[1]+1)*x[1]  x[6]*x[4]^2*(x[2]+1)*x[2]
F[3] = M[3]  x[5]*x[3]^3*(x[1]+2)*(x[1]+1)*x[1]  x[6]*x[4]^3*(x[2]+2)*(x[2]+1)*x[2]
F[4] = M[4]  x[5]*x[3]^4*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1]  x[6]*x[4]^4*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2]
F[5] = M[5]  x[5]*x[3]^5*(x[1]+4)*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1]  x[6]*x[4]^5*(x[2]+4)*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2]
F[6] = M[6]  x[5]*x[3]^6*(x[1]+5)*(x[1]+4)*(x[1]+3)*(x[1]+2)*(x[1]+1)*x[1]  x[6]*x[4]^6*(x[2]+5)*(x[2]+4)*(x[2]+3)*(x[2]+2)*(x[2]+1)*x[2]
end
As a test, I ran the following example, whose correct solution is (approximately):
\alpha_1 = 0.5, \alpha_2 = 2.0, \theta_1 = 1.0, \theta_2 = 1.5, w_1 = 0.3, w_2 = 0.7
M = [2.250, 9.675, 57.263, 427.219, 3836.109, 40234.852] # moments
initial_guess = [0.5; 0.5; 0.5; 0.5; 0.5; 0.5]
nlsolve((F, x) > f!(F, x, M), initial_guess, autodiff=:forward)
Output:
Results of Nonlinear Solver Algorithm
* Algorithm: Trustregion with dogleg and autoscaling
* Starting Point: [0.5, 0.5, 0.5, 0.5, 0.5, 0.5]
* Zero: [1282.030044474002, 1.3349431507591707, 0.013431556054807586, 1.5969473131487903, 0.00013193288212788804, 1.4550855337869975]
* Infnorm of residuals: 10.878979
* Iterations: 1000
* Convergence: false
* x  x' < 0.0e+00: false
* f(x) < 1.0e08: false
* Function Calls (f): 1001
* Jacobian Calls (df/dx): 573
Unfortunately this is not close to the correct solution, which leaves me with two questions:

Is there a way to incorporate the bounds/constraints on the unknowns (i.e., that they are all positive and that the two weights sum to 1)? This would at least rule out “impossible” solutions.

In general, are there alternatives to NLsolve that would be better suited to this problem and would give more accurate and reliable approximations to the solutions, at a similar or higher speed?
Thanks for any help!