I have a non-linear system of equations in Matlab and I would like to find the set of its solutions. The system is represented below.
Imagine for a moment that the Matlab function mvncdf can be used together with solve. My system of equations looks like this
clear
%SOME VALUES USED BELOW
U1true=1;
U2true=-1;
U3true=-2;
U4true=1;
rhotrue=-0.5;
mutrue=[0 0];
sigmatrue=[2*(1-rhotrue) 1-rhotrue; 1-rhotrue 2*(1-rhotrue)];
%SET VARIABLES
syms U1 U2 U3 U4 rho
%SOLVE SYSTEM WITH 6 EQUATIONS
S=solve(mvncdf([ U1 -U2+U1], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U1true -U2true+U1true], mutrue, sigmatrue) ,...
mvncdf([ U2 U2-U1], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U2true U2true-U1true], mutrue, sigmatrue),...
mvncdf([-U1 -U2], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([-U1true -U2true], mutrue, sigmatrue),...
mvncdf([ U3 -U4+U3], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U3true -U4true+U3true], mutrue, sigmatrue) ,...
mvncdf([ U4 U4-U3], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([ U4true U4true-U3true], mutrue, sigmatrue),...
mvncdf([-U3 -U4], [0 0], [2*(1-rho) 1-rho; 1-rho 2*(1-rho)])==mvncdf([-U3true -U4true], mutrue, sigmatrue));
%DISPLAY SOLUTIONS
M=[S.U1,S.U2,S.U3,S.U4, S.rho];
This approach clearly does not work for many reasons, most importantly because mvncdf cannot be used together with solve, as discussed here for example.
As an alternative method to solve my system, I thought about using slicesample as follows
This code implements my idea.
%MAIN
clear
rng default
Utrue=[1;-1;-2;1];
rhotrue=-0.5;
number = 100000;
initials=[Utrue;rhotrue]; %[U1;U2;U3;U4;rho]
rnd = slicesample(initials,number,'pdf',@dens);
%ANCILLARY FUNCTIONS TO SAVE
function density=dens(param)
T=10^(-6);
if param(end)>=1 || param(end)<=-1 %rho outside (-1,1)
density=0; %force density to be equal to zero
else
Utrue=[1;-1;-2;1];
rhotrue=-0.5;
mutrue=[0 0];
sigmatrue=[2*(1-rhotrue) 1-rhotrue; 1-rhotrue 2*(1-rhotrue)];
truth=mvncdf([ Utrue(1) -Utrue(2)+Utrue(1);...
Utrue(2) Utrue(2)-Utrue(1);...
-Utrue(1) -Utrue(2) ;...
Utrue(3) -Utrue(4)+Utrue(3);...
Utrue(4) Utrue(4)-Utrue(3);...
-Utrue(3) -Utrue(4)], mutrue, sigmatrue);
density=exp(-criterion(param,truth)/T);
end
end
function Q=criterion(param, truth)
U=param(1:end-1);
rho=param(end);
mu=[0,0];
sigma=[2*(1-rho) 1-rho; 1-rho 2*(1-rho)];
candidates=mvncdf([ U(1) -U(2)+U(1);...
U(2) U(2)-U(1);...
-U(1) -U(2) ;...
U(3) -U(4)+U(3);...
U(4) U(4)-U(3);...
-U(3) -U(4)], mu, sigma);
Q=(candidates(1)-truth(1))^2+(candidates(2)-truth(2))^2+(candidates(3)-truth(3))^2+...
(candidates(4)-truth(4))^2+(candidates(5)-truth(5))^2+(candidates(6)-truth(6))^2; %this is zero at a solution of the system
end
Question: using the approach that I propose, as I increase the variable number I get more and more solutions for different values of rho, which is in line with comment above that the system has one and only one solution with respect to U for every value of rho. However, to span the entire (-1,1) interval for rho it seems to me that I would need to set the variable number immensely high. Do you have (i) suggestions that can help me to improve the code above to solve this issue, or (ii) better approaches to provide the set of solutions of my system of equations?
