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 2-component 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: Trust-region 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]
* Inf-norm of residuals: 10.878979
* Iterations: 1000
* Convergence: false
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: 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!